ANABIOPLUS: A NEW PACKAGE FOR PARAMETER ESTIMATION AND SIMULATION OF BIOPROCESSES

– A very important step in the study of bioprocesses is the modeling and simulation of bioreactors. This paper presents the AnaBioPlus software package, which was developed for teaching and research purposes concerning the analysis of bioreactors. The package consists of two programs: OptimusFerm and SimulaFerm. OptimusFerm allows the analysis of experimental batch cultivation data sets in order to estimate parameters of kinetic models by means of a differential evolution method. SimulaFerm performs simulations obtaining values of cell growth, product formation, and limiting substrate consumption for batch, fed-batch and continuous cultivations. Thirty-two kinetic models of cell growth are available in OptimusFerm and SimulaFerm. Two case studies are described in order to illustrate the use of AnaBioPlus. The results show that AnaBioPlus is a very useful tool for the analysis of bioprocesses. Furthermore, it presents a user-friendly graphical interface and is freely available to users, in both English and Portuguese versions.


INTRODUCTION
Bioprocessesare present in many industrial processes, including the production of foods, antibiotics, amino acids, organic acids, and ethanol. They are also used in wastewater treatment. Bioprocess reactions are often carried out in bioreactors, involving consumption of substrate, formation of product, and growth of the microorganism responsible for the bioprocess. Therefore, modeling of bioreactors is an important stage in the evaluation of biochemical processes.
In bioprocess modeling, mathematical equations are used to describe the specific growth rate of biomass, which is used in mass balance equations. Many mathematical models are available in the literature to describe cell growth rates, with and without substrate, as well as product and/or cell inhibition (Monod, 1942;Moser, 1958;Contois, 1959;Aiba et al., 1968;Andrews, 1968;Levenspiel, 1980;Hoppe and Hansford, 1982;Lee et al., 1983;Wu et al., 1988), and can be used in simulation modeling of bioprocesses.
Several software packages have been developed for simulation of biochemical processes. OptFerm is a computational platform for the simulation and optimization of bioprocesses,where the objective is to calculate the best concentration of nutrients to initiate the fed-batch bioprocess. Moreover, the user can perform an estimation of kinetic and yield parameters in order to obtain the best possible model setup (Rocha et al., 2009). BioMASS is a simulator that allows users to specify the feed, bioreactor type, operating conditions, and cell growth kinetic model. The program then predicts substrate, product, and biomass concentrations over a given period of bioreactor operation (Phan-Thien, 2011). KINSOLVER is a simulator that uses a representation of the chemical reactions as the input file, instead of a system of differential equations. Concentrations of participating species are specified at the initial time, and the system is solved for concentrations at a final time. Kinetic parameters and initial conditions are all part of the input file (Aleman-Meza, 2011). KINSIM is software for simulation of the progress of kinetic curves, where the user inserts a representation of a chemical reaction scheme using conventional chemical format, after which differential equations representing the mechanism are generated and a direct numerical integration is performed. KINSIM software can be used to simulate any number of reaction mechanisms (Barshop et al., 1983). Commercial software can also be used for the simulation of biochemical processes, employing packages such as Aspen Plus ® orSuperPro Designer ® . EMSO (Environment for Modeling, Simulation, and Optimization) software is an equation-oriented process simulator for modeling dynamic or steady-state processes, distributed at no cost for academic purposes. The software has its own modeling language, and new models can be developed by the users (Soares and Secchi, 2003). The simulator platform provides a tool for parameter estimation. However, the license costs of the Aspen Plus and SuperPro Designer software packages and the need for previous knowledge in using the EMSO software are important considerations. When the focus is on estimating the parameters of kinetic models, programs that make use of numerical integration of rate equations are available, such as FITSIM (Zimmerle and Frieden, 1989), besides commercial software packages that have tools for solving optimization problems, such as Matlab ® , Scilab ® , MathCad ® , Berkeley Madonna ® , and Scientist ® .
The Global Kinetic Explorer (Johnson, 2009), COPASI (Hoops, 2006), KINFITSIM (Svir et al., 2002), and GEPASI (Mendes, 1983)tools simulate biochemical reactions and estimate reaction kinetics parameters. The commercial Global Kinetic Explorer program is available as a professional version and as a free student version that has limited resources. COPASI is a successor to GEPASI and is free for non-commercial use. The KINFITSIM commercial software is available at a discount for academic purposes. Several other important software tools are being developed for use in cell biology, such as CELL (Tomita et al., 1999), Virtual Cell (Loew and Schaff, 2001), E-MCell (Stiles and Bartol, 2001), StochSim (Novère and Shimizu, 2001), STOCKS (Kierzek, 2002), BioSPICE (Garvey, 2003), and ETA (Murphy and Young, 2013).
As can be seen, the literature offers various software options with different characteristics. However, no non-commercial software package is available that estimates kinetic parameters of cell growth models, while at the same time simulating different bioprocesses by means of a userfriendly graphical interface. The present work therefore presents AnaBioPlus, a software package for analysis of bioreactors. AnaBioPlus contains two software modules: OptimusFerm and SimulaFerm. These free and flexible tools can be easily used by persons with minimal training. OptimusFerm is software for estimation of the parameters of cell growth kinetic models, and allows the analysis of more than one batch cultivation experimental data set. Once the kinetic model has been chosen, its parameters are estimated by means of the differential evolution method. The values of the estimated parameters are then displayed, for a confidence interval of 95%, together with a graphical representation of the experimental points and the fitted curve.
SimulaFerm is a simulator of bioprocesses in which the user can specify the cultivation mode (batch, fed-batch, continuous without and with recycle, and with tanks in series), the initial and operating conditions, and the cell growth kinetic model. Simulation is then performed of the behavior of the substrate, cell, and product concentrations as a function of time under transient cultivation conditions or at steady state in continuous cultivation. The software provides help topics concerning its operation, bioprocess kinetics, and types of cultivation, including the corresponding equations. The user can simulate the same process under different cultivation conditions and overlap the curves predicting substrate, biomass, and product concentrations. Experimental data can also be added. AnaBioPlus installation is simple, and both software modules have thirty-two kinetic models of cell growth available.

PRESENTATION OF ANABIOPLUS
AnaBioPlus version 1.0 was developed with the objective of helping researchers and students in the analysis of bioreactors. It is a non-commercial software package consisting of two programs: OptimusFerm and SimulaFerm. OptimusFerm is software for estimating kinetic parameters of cell growth kinetic models, and SimulaFerm is a simulator for bioprocesses. The package has a user-friendly interface, which provides the user with easy access to program resources. AnaBioPlus has some routines (executable, *.dll) in Fortran ® . The software was developed using Microsoft Visual Basic ® and is available in Portuguese and English versions.

OptimusFerm
OptimusFerm version 1.0 is free software for estimation of the parameters of cell growth kinetic models. Parameter estimation consists of finding an optimal set of parameters that maximize (or minimize) one or more objective functions, with or without constraints. In this work, the problem to be solved is a nonlinear optimization, where the objective function (Φ) is minimized (Eq. 1). The constraints (Eq. 2) are the upper and lower bounds of the parameters to be estimated, b i , which are defined by the user. The experimental and calculated concentrations of substrate, biomass, and product are normalized in the objective function (Φ) in order to avoid excessiveweightings in the data set. The experimental and calculated values of each dependent variable are scaled using the maximum value found for that variable in the data set. When multiple batch runs are used, the maximum value of each variable is searched among the data set of each batch.
Many approaches can be used for parameter estimation, ranging from classical deterministic methods, such as Levenberg-Marquardt (LM), to newer methods based on ecological and biological behavior, such as evolutionary algorithms (EAs). The latter have many advantages over conventional nonlinear programming techniques, including no requirement for the gradients of the cost or constraint functions, simple implementation, and less chance of becoming trapped in a local minimum (Nelles, 2001;Long et al., 2013).
The evolutionary algorithm class contains several families of methods, which have different advantages and disadvantages. In developing OptimusFerm, a differential evolution (DE) algorithm was used. Differential evolution is a heuristic method for minimizing possible nonlinear and non-differentiable continuous space functions. Evolutionary algorithms have proven effective in solving many engineering optimization problems and have the advantage of being less susceptible to local minima. The DE algorithm has many useful characteristics, such aseasy implementation, requirement for only a few control variables, and robustness (Nelles, 2001). Furthermore, the DE algorithm does not require parameter initialization, in contrast to derivative methods (such as the LM method). In derivative methods, if an initial set of bad parameter values is set, the method could diverge. Moreover, most complex real-world applications can be solved with this method (Storn and Price, 1997).
OptimusFerm contains 32 non-segregated and unstructured kinetic models of cell growth (µ). These models are divided into different categories: no inhibition, substrate inhibition, product inhibition, cellular inhibition, and hybrid inhibition. Table 1 shows some of the models present in the software. (1) (2)

Andrews (1968)
The 4 th order Runge-Kutta method is used to solve the set of ordinary differential equations that describe the mass balance for biomass, substrate, and product in the operationalmode of the bioreactor. Five different types of cultivation are available in the SimulaFerm software (batch, fed-batch, continuous without cell recycle, continuous with external cell recycle, and continuous with internal cell recycle). As an example, a set of mass balances for batch cultivationis described by Eqs. 12-14.
In Eqs.12-14, the terms r X , r D , r S , r P , and r DP , are the rates of cell growth, cell death, substrate consumption, product formation, and product degradation, respectively (g.L -1 .h -1 ). These reaction rates are given in Eqs. 15-19. Substrate consumption is represented by three terms that account for cell growth, product formation, and cell maintenance (Eq. 17).Product formation can be not associated with growth, partially associated with growth,or both(Eq. 18). C X , C S , and C P are the cell, substrate, and product concentrations, respectively (g.L -1 ).
The rates of substrate consumption for biomass production or for cell growth (r SX ), the rate of consumption of substrate for product formation (r SP ), and the rate of consumption of substrate for cell maintenance (r Sm ) are described by Eqs. 20-22.

Wu et al. (1988)
Product Inhibition Hoppe and Hansford (1982) In Eqs.20 and 21, Y G and Y P are, respectively, thestoichiometric yield coefficient from substrate to cells, also known as the true cell yield coefficient, and the stoichiometric yield coefficient from substrate to product, also known as the true product yield coefficient. The term m s (Eq. 22) is the maintenance energy coefficient (h -1 ).
In practice, it is difficult to determine the exact amount of substrate consumed exclusively for cell growth (Y G ) or product formation (Y P ). Therefore, in mass balances it is usual to express all the substrate consumption by global yield coefficients (Y X/S and Y P/S ), which are easier to determine experimentally.
The SimulaFerm software allows simulationsto be performed using these two approaches. If the global yield coefficients are used, the mass balances for cells, substrate, and product in batch mode operation are described by Eqs. 23-25. If the global yield coefficients are used, the program calculates the value ofthe α coefficient. However, the β parameter should be provided by the user. However, the β parameter should be provided by the user. If the true cell coefficients are used (Y G and Y P ), the user needs to provide the values of m s and the α coefficient. The values of the k d and k dp coefficients are also provided by the user.
It is possible to input more than one experimental data and cultivation conditions set in OptimusFerm. A very useful tool in any program, which is also present in OptimusFerm and SimulaFerm, is the option to save and open files. A convenient resource for beginners in the use of OptimusFerm is the help file, which contains topics concerning the theory of the differential evolution method and the software resources, as well as instructions on how to use the software for the first time.
Estimation of kinetic parameters is performed in batch mode, which is the simplest cultivation method. The confidence interval for each estimated parameter is calculated after the convergence of the DE algorithm, taking into account the model chosen, the degree of freedom of the system and the Jacobian matrix, calculated numerically using finite differences. Good descriptions of the theory of parameter estimation for non-linear model scan be found inHimmelblau (1970), Bard (1974), and Nelles (2001).However, simulations can be performed in SimulaFerm with different cultivation conditions, using the cell growth kinetic model parameters previously estimated in OptimusFerm. Fig. 1 shows the OptimusFerm splash screen. Fig. 2 illustrates the flow of information throughout the OptimusFerm software.  (Silva et al., 2003(Silva et al., , 2005. The version presented in this paper has improvements in layout and source code and adds new graphical features, and for these reasons the name of the software has been changed to SimulaFerm. The package was named AnaBioPlus to differentiate it from the older AnaBio simulator. SimulaFerm performs bioreactor simulations, obtaining curves forcell growth, product formation, and limiting substrate consumption. These are presented as a function of time (for cultures operated under transient conditions), and as a function of dilution rate or residence time (for continuous cultivations operated at steady state).
Users can choose among seven different types of bioreactor configuration, such as batch, fed-batch, continuous without recycle, continuous with external cell recycle, continuous with internal cell recycle, multistage continuous without cell recycle, and multistage continuous with cell recycle. Cultures in batch and fed-batch modes are analyzed in the transient state (graphs of concentration versus time). Continuous cultures with or without cell recycle may be analyzed in the transient state (graphs of concentration versus time) or in steady state (graphs of concentration versus dilution rate or residence time), and the user can view the operating points in steady state. Multistage cultivations with or without recycle allow the user to view the operating points in steady state for three tanks connected in series. The 4 th order Runge-Kutta method is used to solve the set of ordinary differential equations, and the Newton-Raphson method is used to solve algebraic equations. Table 2 shows mass balances for substrate, biomass, and product in different bioreactor configurations provided in SimulaFerm. Fig. 3 shows the SimulaFerm splash screen. Fig. 4 illustrates the flow of information throughout the SimulaFerm software.
Thirty-two cell growth kinetic models are available, some of which are shown in Table 1. The user can choose among product formation types: product formation associated with cell growth; product formation not associated with cell growth; and product formation partially associated with cell growth. The program has operating conditions available with or without cell maintenance, cell death, and product degradation.
SimulaFerm has the option of including experimental points along with the simulated profiles for batch and fed-batch configurations. An additional feature is the analysis of cell growth models, in which the program allows animation of the curve of μ versus C S when kinetic parameters of the cell growth model vary (such as μ max and K S for the Monod model).This tool is of great importance, especially for teaching purposes.

CASE STUDIES
To illustrate the use of OptimusFerm and SimulaFerm, two examples are presentedthat involve the estimation of kinetic parameters and simulation of bioprocesses.

Parameter estimation using OptimusFerm (Case Study 1)
This case study concerns the bioprocess described by Baptista-Neto et al. (2000) for clavulanic acid production. Table 3 presents the experimental values of cultivation conditions and kinetic parameters estimated by the authors using a non-linear regression method (Marquardt algorithm) with 95% confidence level. The batch experiment was conducted using Streptomyces clavuligerus NRRL 3585, with glycerol as the carbon and energy source. The product Fed batch Continuous with external cell recycle

Cultivation Equation
Multistage continuous without recycle Multistage continuous with cell recycle    (Table 1), in which m max is the maximum specific growth rate (h -1 ) and K S is the saturation constant of the Monod model (g.L -1 ).
When starting the program, the window "Insert experiment" will open and the user can add the experimental data of the batch cultivation. Fig. 5 shows a window with the experimental substrate, biomass, and product concentration data, as a function of time, as presented by Baptista-Neto et al. (2000). In this example, only one experimental data set was used, but it is possible to input additional data sets. OptimusFerm has a variety of graphical resources available, including color for the markers of the experimental data, line style for the fitted curve, copy and save the graph, and remove or insert gridlines and legend.
Each experiment added must be saved, or it is not used in the estimation. The option "visible graph" may be selected or not, but if it is not selected, the graphs of the experiment will not be exhibited.
The graphs can be enlarged by double clicking on them or by right clicking and selecting Zoom. Fig. 6 shows the window that then opens with the enlarged graph, where it is possible to insert a title, add or remove the legend and the grid lines, delimit the time interval, copy, save, and print the graph.
The next step consists of inserting the cultivation conditions of each experiment. Fig. 7 shows a screen capture of the "Batch Cultivation" window for the problem studied. To choose the parameters to be estimated, it is necessary to select them in the field "Kinetic parameters", add values to start the algorithm, and fill in their respective upper and lower bounds. In this example, m max and K S were the estimated parameters of the Monod model (Eq. 3, Table 1). The experimental data and cultivation conditions can be saved in files with extension "*.abr", using the menus "File/ Save/Experiments" and "File/Save/Cultivation conditions/ Batch", respectively. They can then be opened whenever necessary.
In the field "Choose" (Fig. 7), the user can select the option "Estimate", which estimates the parameters of the kinetic model, or select the option "Simulate", if the user wants to simulate the process with manually inputted parameters. The sum of squared deviations and the estimated parameters are presented in the "Results" field ( Fig. 7).It is possible to stop the estimationby clicking the red button below the Results field. When the estimation is complete, the user can generate a report with cultivation conditions, estimated parameters, and graphs with the experimental data and fitted curves. It is also possible to export the fitted substrate, biomass, and product concentration data using the menu "File/Export/Data (Clipboard)". This enables graphs to be constructed using external software.
Visual inspection of the fitted curves for the substrate, cell, and product concentrations using the differential evolution method (Fig.7) showed that the parameter estimation could be considered of good quality. This was corroborated by the low values of the sums of squared deviations. The estimated parameter values and confidence intervals were m max =0.2056±0.0006 h -1 and K S =3.49±0.29 g.L -1 . The estimated parameters were within the range of values found in the literature (Baptista-Neto et al., 2000). The user should pay attention to the standard deviation of the estimated parameters (at 95% confidence interval), because it can be so wide that it includes values outside the limits informed by the user, or be in a range in Table 3. Bioprocess conditions and estimated parameters obtained by Baptista-Neto et al. (2000) Critical substrate concentration at which product formation begins Cs 2 * 11.14 g.L -1 Constant for non-growth-associated product formation β 0.0035 g P .g X -1 .h -1 The differential evolution method has control parameters that affect the performance of the algorithm. Some of these parameters, such as the population size and crossover probability, can be modified by the user in order to try to improve the estimation. In OptimusFerm, there are default values for these parameters. To change them, the user must access the menu "Fitting method" or use the shortcut Ctrl+M (Fig. 8). However, alteration of these values by the user should be performed with caution and with prior knowledge. There are some studies in the literature where self-adaptation of the parameters of the differential evolution algorithm is described (Ghosh et al., 2011;Pan et al., 2011), but as this is beyond the scope of the present work, the ranges of these values were only suggested and default values were assigned in accordance with the literature. For the estimation of kinetic parameters presented in this case study, the control parameters were not changed.
A convenient resource for beginners of OptimusFerm is the help file, which contains topics concerning the theory of the differential evolution method and the software resources, as well as instructions on how to use the software. The help file is available from the "Help/Help topics"menu. It is important to mention that the software does not handle the problem of ill-conditioned parameters.
The user must be aware of this problem.

Simulation using SimulaFerm (Case Study 2)
This second case study concerns the bioprocessdescribed by Badino and Hokka (1999) for ethanol production. Table  4 presents the cultivation conditions for this problem. The cell and product yield coefficients (Y X/S and Y P/S , respectively) must be calculated separately by independent equations (Eqs. 53and 54) and are part of the data supplied by the user.

Cs Cs
The batch experiment was conducted using Saccharomyces cerevisiae, with glucose as substrate. The product (ethanol) was associated with cell growth. Estimation of the kinetic parameters (m max and K S ) for this problem employed OptimusFerm with the Monod kinetic model (Eq. 3), and the values obtained were within the ranges found in the literature (Hoppe and Hansford, 1982;Badino and Hokka, 1999).
The user must first choose the type of cultivation (batch, fed-batch, continuous without recycle, continuous (53) (54) Figure 8. Window used to modify the control parameters of the differential evolution method (Case Study 1) with recycle, and continuous in tanks connected in series without or with recycle), the type of analysis (graphical analysis for transient operations, or graphical analysis and operations at steady state), the kinetic model of cell growth, and the production type (associated, partially associated, or not associated with cell growth). After selection of the model and production kinetics, the cultivation conditions are inserted as shown in Fig. 9. In this window, the user can view the profiles of the substrate, cell, and product concentrations over time ("View graph"), and simulate the process using different conditions ("Overlay graphs"). It is also possible to overlap the simulated curves in different colors, add experimental data ("Experimental points"), delete inserted data ("Clear data"), return to the step of choosing the cell growth model and the production type ("Back"), and exit the cultivation window ("Exit"). These options are also available from the toolbar. All cultivation condition windows provide an option to view the cell growth model that is being used.
The experimental data of substrate, cell, and product concentrations, as a function of time, can be added using the option "Experimental points" (Fig.  10), when selecting batch and fed batch cultivation modes. To illustrate this step, experimental data were added in Case Study 2. The user can choose which experimental points will be added, the exhibition mode, and the color and symbol used for the experimental datamarkers. After adding the experimental points, the cultivation conditions window is reopened (Fig. 11), displaying the graphs containing the experimental data and the simulation curves. The quality of the simulation can be seen from Fig. 11, since the experimental data of Badino and Hokka (1999) lie very close to the simulation values. In addition, the software allows the user to modify the kinetic parameters in order to improve the visual fitting of the kinetic model to the experimental data, enabling observation of the sensitivity of the model with respect to the kinetic parameters.
SimulaFerm has options available to copy and save the graph, remove or insert gridlines, and enlarge the graph. The simulation data can be saved in files with extensions "*.abr" or "*.dat", using the menu "File/Save simulated data as". This enables the data to be used in other programs (such as Microsoft Excel ® or Origin ® ). The input data (cultivation conditions) can also be saved for future use.
The program includes a help file containing information about its resources and the theory of bioprocess kinetics and bioreactor analysis. The animation of cell growth models is another option that can aid students and researchers, because it is possible to choose one of the available cell growth models and vary its parameter values, at the same time observing the results of the changes in the graph of μ versus C S . The user can compare the results of the chosen model with those of the classical Monod model.
The AnaBioPlusprogram is available for download at the following website: http://www.deq.ufscar.br/ downloads/anabioplus.zip

CONCLUSIONS
The AnaBioPlus software package has great potential for teaching and research purposes involving the analysis of bioreactors. It is a non-commercial software package consisting of two programs: OptimusFerm and SimulaFerm. OptimusFerm is software for estimating the kinetic parameters of cell growth kinetic models, while SimulaFerm is a simulator for bioprocesses in different types of bioreactor configuration. The use of OptimusFerm and SimulaFerm was illustrated here by considering two example problems involving the estimation of kinetic parameters and simulation of a bioprocess. Good results were obtained for the estimated parameters and the simulation of the batch mode bioprocess, indicating that AnaBioPlus is a very useful tool for both purposes. Furthermore, the package has a user-friendly interface, which allows the user easy access to program resources. upper bounds of the parameters (-) c: fraction of outlet stream leaving the bioreactor without passing through the cell separation unit (-) C P : product concentration (g.L -1 ) C P0 : initial product concentration or inlet product concentration (g.L -1 ) C P * : critical value of product concentration to product inhibition (g.L -1 ) C P1 : product concentration at bioreactor outlet stream in continuous cultivation with external recycle (g.L -1 ) C Pn : product concentration in tank n (g.L -1 ) C Pn-1 : product concentration in tank n-1 (g.L -1 ) C S : substrate concentration (g.L -1 ) C S0 : initial substrate concentration or inlet substrate concentration (g.L -1 ) C S1 * : critical substrate concentration at which cell death begins (g.L -1 ) C S2 * : critical substrate concentration at which product formation begins (g.L -1 ) C Sf : substrate concentration in the feed (g.L -1 ) C Sn : substrate concentration in tank n (g.L -1 ) C Sn-1 : substrate concentration in tank n-1 (g.L -1 ) C X * : critical value of cell concentration to cell inhibition (g.L -1 ) C X : cell concentration (g.L -1 ) C X0 : initial cell concentration or inlet cell concentration (g.L -1 ) C Xn : cell concentration in tank n (g.L -1 ) C Xn-1 : cell concentration in tank n-1 (g.L -1 ) DE: differential evolution EAs: evolutionary algorithms F: culture medium feed volumetric flow rate (L.h -1 ) h: decline factor of cellular concentration to pass through the separation unit