In-silico Optimization of a Batch Bioreactor for mAbs Production in Relationship to the Net Evolution of the Hybridoma Cell Culture

Bioreactor optimization is a common engineering problem difficult to be solved due to the large number of influential variables of high variability. Production of monoclonal antibody is a well-known method to synthesize a large number of identical antibodies (that is of uniform characteristics, also called monoclonal antibodies, mAb). Due to such reasons intense efforts have been invested to maximize the production of mAb by using hybridoma cell culture. Based on an adequate kinetic model from literature (experimentally checked) this paper focus on pointing-out the major role of the net evolution of the viable biomass (growth, and decay) in the location of the optimal operating setpoint (SP) of a three-phase mechanically agitated batch bioreactor (TPMAB) with immobilized hybridoma culture. This in-silico analysis opens the possibility I) to optimize the bioreactor performances by placing the SP in the most favourable location, by adjusting the substrate and biomass initial load in the bioreactor according to the preliminary determined characteristics of a modified / improved biomass; ii) to optimize the batch-to-batch operation mode (not approached here) according to the time-varying characteristics of the biomass culture, or iii) to determine the optimal operation of the bioreactor in a fed-batch operating mode (not approached here).

Over the last decades, there is a continuous trend to develop more and more effective bioreactors [1,2] to implement various biosyntheses for producing finechemicals or organic compounds in the food, pharmaceutical, or detergent industry, by using freesuspended or immobilized cell cultures in batch, semibatch (fed-batch), continuously operated fixed-bed, or continuously three-phase fluidized-bed (and/or mechanically agitated) reactors, aiming to replace complex chemical processes, energetically intensive and generating toxic wastes [3,4]. It is to be mentioned, among these, the fermentative processes for production of organic acids, alcohols, vinegar, amino-acids/proteins, yeast, hydrogen, food products and additives, recombinant proteins/ antibodies, etc., by using bioreactors with microbial / cell cultures [2,3], or bioreactors for waste removal in water treatment plants [5], with various cooling and aeration systems in simple or sophisticate constructive alternatives [2], by integrating genetic and engineering methods [6,7].
One application of tremendous importance in medicine is the industrial production of monoclonal antibodies (that is of uniform characteristics), by the so-called Hybridoma technology, [8][9][10] by using antibody-secreting hybridoma cell cultures. However, the large-scale production of mAbs by using mammalian cells in optimized batch, or fedbatch culture systems (that is with continuous feeding of GLC, and GLN substrates, or nutrients with variable concentrations/feeding rates [10][11][12] is limited by engineering problems related to associated to the modelbased optimization, that is: a) The presence of multiple (often contrary) objectives [13]; b) An important degree of uncertainty coming from multiple sources [14][15][16][17][18][19], such as: model inaccuracies, variability of biomass activity, presence of disturbances in the operating variables; nonlinear dynamics of the bioprocess of low reproducibility.
And, most important, as remarked by Dorka [12]; Li et al. [20], the large-scale production of monoclonal antibodies (mAb) by mammalian cells in batch and fedbatch culture systems is limited by the unwanted decline in cell viability and reduced productivity that may result from changes in culture conditions. Therefore, it becomes imperative to gain an in-depth knowledge of the factors affecting cell growth and cell viability that in turn determine the antibody production [20][21].
To support engineering calculations, several attempts have been reported to obtain adequate dynamic models to predict the behaviour of both batch, or fed-batch systems as a function of the extra-cellular nutrient/metabolite concentrations. Such model formulations will aid in identifying and, eventually, in controlling the dominant factors involved in the optimization of the mAb production. More sophisticated bioprocess kinetic models try to also explicitly include intra-cellular factors related to the cellular metabolic fluxes [11,12,22], associated to the so-called central carbon metabolism (CCM) [23]. Other optimization approaches are focus on improving the fed-batch operation by also accounting for the uncertainty in the model structure and its parameters, that is, the so-called robust (stochastic) optimization; see the review [24]. Most of the developed kinetic models account for the inhibition/limitation effects on the reaction rates due to various intra-/extra-cellular metabolites, substrates, or by-products [25].
By adopting an adequate kinetic model (experimentally checked) from literature [24], this paper is aiming at insilico investigating the influence of various operating variables (that is, the initial glucose GLC, glutamine GLN, and viable biomass X v ) on the performances of a TPMAB with immobilized hybridoma culture on porous alginate beads [26,27]. As an element of novelty, an increased attention is given to the major role played by the initial load, and the net evolution of the viable biomass (growth, and decay) during the batch in the optimal operating setpoint location. Such an in-silico analysis is done by adding a term for enhanced biomass growth and decay to the Liu and Gunawan model [24], thus deriving optimal operating points.
This in-silico analysis opens the possibility to optimize the bioreactor performances by placing the setpoint (SP) in the most favourable location, by adjusting the substrate and biomass initial load according to the characteristics of the modified / improved biomass with an enhanced evolution (with the below k net ≠ 0). This in-silico analysis is also a valuable tool to determine an optimal operation of the bioreactor based on the timevarying characteristics of the biomass culture during the batch [18,28], or during a batch-to-batch operation mode (not approached here [29] ), or by adjusting the biomass load in the bioreactor according to the time-varying characteristics of the biomass culture, in a semi-continuous operating mode (not approached here [30] ).
A separate preliminary determination of the biomass kinetic characteristics (see the below net k net ≠ 0) in an incipient stage of the bioprocess can lead to determine the optimal operating policy of a batch (this paper), or of a fed-batch bioreactor [28,31]. Optimization analysis is realized gradually, by considering one control variable (i.e. initial concentration of GLC, or of GLN), or by simultaneously considering two control variables (i.e. initial concentrations of GLC, and GLN), or three control variables (i.e. initial concentrations of GLC, GLN, and of biomass X ν,o = X t,o ). Calculations are made with considering the biomass basic characteristics of [24] (with the below k net = 0), or considering a modified biomass with an enhanced evolution (k net ≠ 0).
The parametric uncertainty associated to the biomass dynamics models (not specified by [24]), is neglected here to not complicate the computational analysis. Due to such a reason, the predicted optimal operating policy by the math model cannot be realized with a perfect certainty. However, when significant inconsistecies vs. experimental evidence are reported, such checks will lead to the process lumped model updating/completions, and re-estimation of some of their parameters for correcting its adequacy in order to perform futures bioreactor analyses, and to eventually correct its optimal operating policy. To partly avoid such problems various alternatives are suggested in the literature, such as an observer-based robust control strategy for controlling overflow metabolism cultures operated in fed-batch mode in order to maximize the biomass productivity [18,[28][29]31], even if their implementation is difficult and costly.

Bioprocess / bioreactor dynamic model Bioprocess kinetic model
A typical stirred tank bioreactor is equipped with temperature, pressure, agitation, pH, and DO control systems [1][2]. The relatively simple dynamic ideal model of the TPMAB describing the key-species dynamics is presented in Table 2. The ideal model main hypotheses are the followings: I) isothermal, iso-pH, iso-DO operation; ii) nutrients [25] are added initially and during TPMAB operation to ensure the optimal maintenance of the biomass, in a recommended C/N/P ratio of ca. 100/5/1 wt. [5]; iii) aeration in excess over the batch to ensure an optimal biomass maintenance; iv) perfectly mixed liquid phase, of quasi-constant volume; v) negligible mass resistance to the transport of nutrients/substrates/products in the porous alginate beads; vi) negligible mass resistance to the transport of oxygen into the liquid and beads; vii) substrates, and the solid carrier (of a millimeter size) including the immobilized biomass are initially added in the batch bioreactor; viii) the solid particles (of uniform characteristics) are considered uniformly distributed in the homogeneous liquid phase, due to the well-mixing conditions.
The batch bioreactor is considered here in the alternative with using immobilized biomass from several reasons: I) due to the reported biomass higher stability [26,32], and ii) to eventually consider its fed-batch operation {not approached in this paper [30]}. The adopted kinetic model for mAb production by using hybridoma cell culture is presented in Table 2 together with the associated rate constants. This dynamic model [24], was experimentally validated by [33], and others. Due to such a reason, one expect that the present in-silico analysis of the batch bioreactor to provide results of a satisfactory predictive value.
From the mathematical point of view, the bioreactor dynamic model translates in a set of differential mass balances for every considered species, of the form: (1) with the reaction rate r j expressions given in Table 2. To get the key-species (j) dynamics over the batch time (t f ), the model (1) is solved with a proposed initial condition C j,o = C j (t=0), under the best medium conditions of Table 1. To get the solution with enough precision, a low-order stiff integrator (ODE23S) of the Matlab™ package was used. Biomass dynamics model As the biomass dynamics plays the principal-role in the all reported kinetic models for the production of mAb by using hybridoma cell cultures, the rate expressions proposed in the literature for the biomass growth/death are of a Monod-type, by including terms accounting for inhibitions and limitations given by metabolites, substrates, or by-products [11,24,[34][35][36], i.e: (2) Terms f(ξ(t)), g (ξ(t)) are complex functions accounting for the biomass growth, or death inhibition respectively, due to various extra-cellular species, i.e. GLC, GLN, LAC, AMM, etc. [12,24,[34][35][36], or intra-cellular metabolites [22] related to the so-called central carbon metabolism (CCM). More structured models by accounting for CCM are also developed by [37].
Other authors [25,35] proposed empirical kinetic forms for the cell growth and death. Regardless the approach, the initial cell density and the medium characteristics clearly play the central role on the hybridoma growth, metabolism, and mAb production yield [22,38].
To better highlight this prevailing role plays by the biomass growth/death rates on how to choose the bioreactor optimal operating conditions, one supplementary k net term has been introduced in this paper in the biomass balance, the eqn.(2) becoming: (3) where the lumped k g X ν , and k d X ν terms, include the enhanced biomass specific growth rate (k g ), and the enhanced biomass specific death rate (k d ). These two  [24] terms account for conditions others than those considered in the usual terms f(ξ(t)), g (ξ(t)). Several values of the net biomass evolution rate k net = k g -k d will be tested in this study, according to the literature information. The biomass enhanced evolution might be explained by lot of operating variables not included in the common dynamic models of eq.(2) form. Such variables could include [9,36]: the biomass X n quality, immobilization type, etc [9,26,27,39].
To in-silico study how the location of the optimal setpoint in the parametric space is directly influenced by the biomass evolution, a term k net accounting for the enhanced (modified / improved) biomass characteristics (growth, and decay) was added to the Liu and Gunawan [24] model, and the derived optimal operating points were compared.
Common values for the biomass maximum specific growth rate (k g ) and maximum specific death rate (k d ) are given in Table 3. According to these reported values, the following values for the k net will be tested in the present study (in 1/h units to be in agreement with the similar constants of Table 2): A separate determination of the biomass kinetic characteristics (that is the net k net ≠ 0) in an incipient stage of the bioprocess can lead to determining an optimal operating policy of a batch or of a fed-batch bioreactor. The bioreactor nominal setpoint SPN, presented in Table 1, is those of [24]. For this reference setpoint, the biomass evolution is governed by only f(ξ(t)), g (ξ(t)) terms of eq.(2). When significant inconsistecies between predicted optimal operating policy vs. experimental evidence are reported, such checks will lead to the process lumped model updating/completions, and re-estimation of some of their parameters for correcting its adequacy in order to perform futures bioreactor analyses, and to eventually correct its optimal operating policy.

Formulation of the optimization problem
For the present case study, to determine the optimal running conditions of the TPMAB of Table 1, that maximize the mAb production by the used hybridoma cell culture, several control variables could be considered. These parameters are physical, chemical and biological in nature [20], as followings: I) inoculum size, that is initial cell density, X ν0 [22]; ii) substrates initial concentrations, [GLC](0), [GLN](0) [24]; iii) an optimal feeding policy of the fed-batch bioreactor with substrates and nutrients to control the biomass evolution (not for the batch case here) [12]; iv) an optimal feeding policy with immobilized viable biomass X ν,0 (t) of the fed-batch bioreactor, with an incremental addition policy to compensate the exponential decay of the viable biomass in the bioreactor due to its degradation and/or washout from the alginate-carrier (not for the batch case here) [30].
Because not all the mentioned control variables are included in the adopted bioprocess dynamic model [24], the present paper will be focus on determining the following three control variables under certain operating alternatives: I) substrates initial concentrations, [GLC](0), [GLN](0); ii) inoculum size (that is initial cell density, X ν,0 (t); due to inoculum preparation reasons, the total initial biomass will be considered equal to the viable one, X ν,0 (=X t,0 ).
In mathematical terms, the optimization problem can be formulated as following: where the mAb(t) time-evolution is determined by solving the bioreactor dynamic model eqn.(1) with a proposed initial condition C j,o =C j,j (t=0), and for an imposed batch time (t f ), and optimal medium conditions (temperature, pH, DO, nutrients) of Table 1. Notations: C = species conc. vector; C o = species initial conc.; k= model rate constants.

In-silico analysis results
Optimization analysis is realized gradually, by considering one control variable (that is the initial concentration of GLC, or GLN), or simultaneously two (initial concentrations of GLC, and GLN), or simultaneously three control variables (initial concentrations of GLC, and GLN, and of X ν,0 =X t,0 ) with the biomass basic characteristics of [24]. Finally, the optimal policy given by searching over the three-control variables was reconsidered by using a modified / improved biomass with an enhanced evolution (with ), and the results compared.
The reference setpoint for our in-silico analysis is the bioreactor nominal setpoint SPN [24] presented in Table 1. For this SPN, the biomass evolution is governed by only the common terms of eq.(2), with no enhanced process (that is k net =0). The key-species dynamics is plotted in Figure 1, with a realized Max[mAb](t) = 1254.6 (mg/L) at the batch end t f =103 h (Table 4) in concordance with those reported by [24]. To summarize our calculations, the following optimization alternatives have been performed for k net =0 case: Optimization type 1-D. Control variables: [GLC](0), or [GLN](0), or X ν,0 =X t,0 . Results: SP1, and SP2. Displayed results in Figure 2 and Table 4 Figure 2D and Table 4 Displayed results in Figure 3 and Table 4. Observation: Optimal policy SP4 indicates the best bioreactor performance. For 2-D and 3-D optimization problems (5), the effective multimodal optimization solver MMA of Maria [40,41] has been used implemented in the Matlab™ computation platform, to prevent the search collapse in a possible local optimum.
For the biomass enhanced evolution (with k net ≠ 0) case, several alternatives have been investigated, corresponding to certain medium conditions, not explicitly considered in the adopted dynamic model [24], or due to a modified biomass. The enhanced evolution refers either to an enhanced biomass growth with a maximum specific growth rate k g higher than µ max in Table 2, or to an enhanced biomass decay due to its degradation and/or washout from the carrier with a maximum specific decay rate k d higher than µ d,max in Table 2.
Such medium conditions not considered in the bioreactor model could be each of the following parameters: the temperature, the gas flow rate, agitation speed, DO, dissolved CO2, pH, osmolality, waste by-products, etc. [20].
To account for such alternatives, k ne values of eqn. (4) have been tested. Simulation results have been comparatively presented in Figure 4, in terms of the realized Max(mAb) dependence on the initial conditions [all three control variables (5)], for k net =0 (SPN), k net =+0.1 (enhanced biomass growth), and k net =-0.1 (quick biomass decay).
By analysing all the optimization result alternatives , several conclusions can be derived:    [24], for an enhanced biomass evolution, that is for: (black line 1) k net =0 (SPN of [24], (blue line 2) k net =+0.1 (enhanced biomass growth), (red line 3) k net =-0.1 (quick biomass decay). Notation (0) indicates the initial condition. Plots of are made: (A) function of initial , for = 4.9 mM, = 2×10 8 Cell/L; (B) function of initial , for = 29.1 mM, = 2×10 8 Cell/L; (C) function of initial , for = 4.9 mM, = 29.1 Mm; (D) function of both initial , and initial , for = 2×10 8 Cell/L Fig. 3. Key species dynamics for the optimal setpoint SP4 (with k net =0) of the batch bioreactor, with the initial conditions specified in Table 4 a)the 1-D optimal setpoints SP1-SP2 reported better performances compared to the SPN of [24], even if only one control variable was considered; b)the 3-D optimal setpoint SP4 reported better performances compared to SPN, or SP1-SP3, because all control variables were simultaneously considered (default k net =0); c)the results displayed in Figure 4, for k net ≠ 0 proved the key role of the biomass evolution regarding the bioreactor performance. Thus, an enhanced biomass growth with k net =+0.1 reported performances better with several orders of magnitude compared to the SPN with k net =0. On the contrary, an enhanced biomass decay with k net =-0.1 reported performances with at least one-order of magnitude worse compared to the SPN with k net =0. d)By comparing the 2-D optimal operating region location Max(mAb) given in Figure 2D for SPN, and in Figure  4D for k net =-0.1 (quick biomass decay), it clearly appears that location of the optimal setpoint is moving according to the net evolution of the biomass (k net ≠ 0), even if the initial load of biomass (X ν,0 = X t,0 ) is the same (those of Table 1).
To partly avoid optimal policy implementation problems, various alternatives are suggested in the literature. For instance, an observer-based robust control strategy for controlling overflow metabolism cultures operated in fedbatch mode in order to maximize the biomass productivity [31,46], or on-line optimal operation based on the time-varying characteristics of the biomass culture during the batch [18,28,30,46], or during a batch-to-batch operation mode [29], or by adjusting the biomass load in the bioreactor according to the time-varying characteristics of the biomass culture, in a semi-continuous operating mode [30].

Conclusions
This paper proves in a relatively simple yet eloquent way how a lumped, but enough detailed and adequate dynamic model of an important bioprocesses can support in silico engineering evaluations, if the cell nano-scale metabolism (including metabolic by-products) is somehow reflected by the bioreactor macro-scale dynamic model.
A separate determination of the k net ≠ 0.1in an incipient stage of the bioprocess can lead to determining an optimal operating policy of a fed-batch bioreactor. Once validated, such valuable model-based predictive analyses will lead to develop a model-based control and optimization of the batch or fed-batch (not approached here) bioreactor operation [11].
In such a way, the significant large experimental and computational effort to elaborate such hybrid cell/ bioreactor models of good quality is fully justified through the benefits of subsequent in silico analyses (as is also the case here).