A Framework for the Development of Integrated and Computationally Feasible Models of Large-Scale Mammalian Cell Bioreactors

Abstract

Industrialization of bioreactors has been achieved by applying several core concepts of science and engineering. Modeling has deepened the understanding of biological and physical phenomena. In this paper, the state of existing cell culture models is summarized. A framework for development of dynamic and computationally feasible models that capture the interactions of hydrodynamics and cellular activities is proposed. Operating conditions are described by impeller rotation speed, gas sparging flowrate, and liquid fill level. A set of admissible operating states is defined over discretized process parameters. The burden on a dynamic solver is reduced by assuming hydrodynamics at its fully developed state and implementation of compartmental modeling. A change in the conditions of operation is followed by hydrodynamics switching instantaneously to the steady state that would be reached under new conditions. Finally, coupling the model with optimization solvers leads to improvements in operation.

1. Introduction

1.1. The Importance of Reliable Unit Operation Models

Global sales of biopharmaceutical products reached $228 billion in 2016 [1]. Compared to microbial and yeast-based production systems, mammalian cells possess the cellular machinery to manufacture and secrete large proteins with the necessary post-translational modifications. Mammalian cell cultures are responsible for half of the revenue generated by the biotechnology industry, which is expected to grow by 15% annually [2].

The average commercial-scale titer for mammalian-expressed products has increased by 10-fold since the early 1990s and reached 2.5 g/L in recent years [3]. This has mainly been achieved via clone selection, cell line screening, host cell engineering, improving vectors, and gene amplification. At the same time, the demand and regulations for production of therapeutic proteins from mammalian cells have also been increasing. As a result, the delivery of affordable products at consistent quality is still a challenge for the industry. The industry has responded to this challenge by expanding the manufacturing capacity and improving the operational agility and efficiency of the supply chain. Bioreactors as large as 25,000 L are used, which has increased the capacity of manufacturing sites up to 200,000 L [4]. However, scale-up methodologies have remained based on the overall characterization of the components of the system such as aeration and agitation mechanisms. Considering the total cost of launching an asset and net cash inflow an asset is forecast to deliver, the return on innovation, as the engine of the industry, has been declining. Deloitte Center for Health Solutions monitored the Research and Development (R&D) performances of 12 major biopharmaceutical companies. The rate of return declined from 10.1% in 2010 to 3.7% in 2016 [5]. The response of the industry has been mainly strategic [6]. The industry has tried to increase the value of the drug pipeline through mergers and acquisitions and the identification of priority customers and emerging markets. It also has pursued a reduction in the cost of launching new assets by developing explicit therapy area focus, balancing in-house and outsourced activities, defining specific missions for the enterprise, utilization of new technologies and advanced analytics in R&D, and data exploitation.

Fulfilling strategic plans depends on the execution of tactical decisions [7]. For example, strategic R&D management involves project selection, budgeting, and commercialization, while tactical R&D addresses the scheduling and resource management necessary for the accomplishment of the project. Integrated decision-making as a means for efficient use of data and knowledge has been highlighted in the “National Strategic Plan for Advanced Manufacturing” published by the Executive Office of the President in 2012. Enterprise-wide decision-making demands the integration of operational activities of planning, scheduling, and control [8]. Control of unit operations needs mechanistic models. These models capture complex reactions and transport phenomena and may demand many CPU hours due to heterogeneity and dynamics of certain systems. Maintaining computational feasibility has been approached by substituting the original detailed model with a reduced model. Reduced models have been developed either by reduction of the order of the model based on evaluation of the significance of its components or development of surrogate models using data obtained through careful experimental design.

1.2. Bioreactor

The goal in the operation of a bioreactor is to enhance the growth, viability, and productivity of organisms by adjusting their environment. Figure 1 presents a cause and effect diagram of bioreactor operation. The main operational parameters are impeller rotation speed, gas sparging flow rate, pH, and temperature. Inefficient mixing creates spatial gradients in mechanical shear, volume fraction of gas phase, dissolved oxygen (DO), carbon dioxide and metabolites concentrations, pH, and temperature. As organisms move inside the reactor they experience fluctuations in environmental conditions that affect metabolism, yield, and quality of product [9]. Due to low solubility in water, high densities of organisms quickly consume all the oxygen in a saturated culture and produce enough carbon dioxide to have inhibitory effects. Therefore, the addition and removal of gases are inevitable parts of fermenter operation. Bioreactor systems can operate under different modes: batch, fed-batch, continuous, and perfusion [10]. Considerations regarding the set-up, operation, and control of bioreactors can be found in the literature [11,12,13,14,15].

The environment refers to the physical and chemical stimuli acting on organisms. Physical stimuli of bubble size distribution and dissipation rate of energy directly affect the viability of cells. Hydrodynamic forces acting on an organism can be quantified using shear stress and energy dissipation rate. Some of the effects of mechanical stress are lysis, change in distribution and number of surface receptors, alteration in production of specific proteins, rate of DNA synthesis, rate of certain metabolic processes, and induction of apoptosis [16]. Chemical stimuli include concentrations of metabolites (e.g., glucose, glutamine, lactate, and ammonia), dissolved oxygen concentration, and pH. Concentrations of other supplements that have functions such as species transport enhancement, growth stimulation, shear protection, and surface charge modification can also be considered. The dynamics of organisms and environment interaction should be modeled based on an understanding of the response times of biological and physical mechanisms. Prokop compared the characteristic times of physical and biological mechanisms [17]. Characteristic time is a measure of the time needed by the mechanism to smooth out a change to a certain extent. It was observed that physical mechanisms of mixing, oxygen transfer to the liquid phase, and diffusion are faster compared to the biological mechanisms of oxygen and substrate consumption.

1.3. Existing Culture Models and Their Need for Improvement

Robust mechanistic bioreactor models facilitate improvement in equipment utilization, medium design, and feeding strategy, and explain causes of scale-up problems such as productivity and byproduct formation [18]. Assuming spatial homogeneity and neglecting the effects of physical environmental parameters have been the basis of most culture models [19,20,21,22,23,24,25,26,27,28,29]. This reduces the problem of simulating relevant portions of cellular activities that control the production of a product of interest. Metabolic models are divided into two main groups based on whether the cell mass composition is considered variable (structured) or fixed (unstructured). Structured models are based on a set of biochemical reactions, which create stoichiometric relations to represent the metabolism of the organism. The reactions contain both extracellular components and intracellular components. The flux of the intracellular metabolites can be determined by assuming a pseudo-steady state inside the cell. Unstructured models utilize a reduced number of reactions to capture metabolism macroscopically. Unstructured models do not have the ability to capture the effects of the growth condition on cellular composition, and as a result they cannot accurately predict balanced growth for a culture in transient condition. In our earlier work examples of different metabolic models and their capabilities have been reviewed [30]. On the unreliability of existing cell culture models, it should be noted that their predictions are independent of whether a shake flask or a large-scale bioreactor is used for the cultivation. These models, in which hydrodynamics is often excluded, usually consider all causes of cell loss in one parameter, e.g., growth rate. The performance of the model deteriorates when the process conditions change as the parameters of unstructured models are highly dependent on strain, cultivation medium, and fermentation conditions [31]. Consequently, existing models do not possess satisfactory predictive power, especially when used outside the calibration range, and literature data can only be used for qualitative studies. In addition to narrow confidence intervals, achieving a satisfactory fit of experimental data often requires considering extra terms or assumptions, which are theoretically difficult to explain. Overall, despite the practical and commercial applications of animal cells, there are only a few reports on their kinetics of growth and production. More specifically, there is no literature report on the kinetic parameters for Chinese hamster ovary (CHO) cell batch culture related to mAb production [32]. Therefore, in this work, instead of final protein concentration, the time-integrated value of biomass is the subject of maximization. The specific attributes of large-scale mammalian cell cultures, e.g., high level of spatial heterogeneity and sensitivity of organisms to physical environmental stimuli, demand a modeling framework that captures both the biology and the hydrodynamics of the system and also their interactions. This work aims to improve on the current state of lumped-parameter bioreactor modeling by capturing the interactions of system components. The selection of components is based on experimental observations. The formulation is devised in order to take into account the inherent dynamics of the system, maintain computational tractability, and couple the model with optimization solvers. In the next section, the modeling and integration of hydrodynamics and biological processes are discussed and the application of compartmental modeling for maintaining computational feasibility is explained. In Section 3, the integrated model is used to find a near-optimal operation scenario. Section 4 is a discussion of the challenges and potential in the area of modeling and the optimization of a bioreactor operation.

2. Development of a Dynamic, Integrated, and Computationally Feasible Bioreactor Model

This work seeks to improve the reliability of bioreactor models by capturing the effects of hydrodynamics on the performance of a bioreactor. Computational Fluid Dynamics (CFD) simulation is employed to calculate the attributes of flow required for this purpose. Biological processes are inherently dynamic; therefore, solving the problem requires dynamic CFD simulations. Dynamic CFD simulations demand the discretization of time using step sizes in the range of 0.01 to 0.1 s [33,34,35,36,37], even for a small reactor. A CPU time of 12 s per CPU for each time step has been reported [34], and 0.5 to 2 s per CPU has been observed for every iteration [36]. Considering the reported CPU times and the fact that bioreactors are usually operated in fed-batch mode for up to two weeks, dynamic CFD simulation of the entire operation is computationally unfeasible. To tackle this problem, in our earlier work a two-step framework was developed [38]. This method is based on the assumption that the effects of metabolic activities on hydrodynamics are negligible [18,39]. The simulation is run in two steady-state and dynamic steps. First the steady state of the two-phase flow inside the bioreactor for specific values of impeller rotation speed and gas sparging flowrate is calculated. At this stage cells are considered merely components of the liquid phase without any biological function. Then the problem is solved dynamically to obtain the evolution of biophase over time. In the dynamic phase of the simulation only species conservation is considered. During the dynamic run, the process parameters of impeller rotation speed and gas sparging flowrate are fixed and the flow remains at its fully developed steady-state condition. Therefore the values of velocity, gas volume fraction, kinetic energy, and energy dissipation rate remain approximately constant. Obtaining the solution in two stages replaces the problem with two smaller systems of equations. This allows for the specification of larger time step sizes and improves the convergence of the simulation.

Deconstructing the problem into steady-state flow and dynamic metabolism problems and solving them sequentially seems to successfully address the batch operation with constant process parameters [38]. Although this simulation procedure exploits off-the-shelf software packages, it is still limited to simple case studies. The improvement in computational feasibility in order to simulate fed-batch operation and study effects of composition and schedule of feeding on the performance of bioreactor is achieved through the development of a compartmental model. Compartmental modeling facilitates time–space decomposition, which significantly reduces computation. Therefore, it has been widely used for modeling the hydrodynamics of stirred tanks [40,41,42,43]. In this methodology the reactor is divided into well-mixed zones that do not contain segregated regions. Then fluxes between the zones are calculated based on the expected flow patterns resulting from CFD simulations or experimental data [18,40,41,44]. A number of states of operation are defined by discretizing the process parameters of impeller rotation speed, gas sparging flow rate, and operating volume. Data on fully developed steady state flow under all pre-defined operation states are obtained and stored in flow matrices. Flow matrices contain information on inlet and outlet fluxes, gas volume fraction, dissipation rate of mechanical energy, and gas superficial velocity of compartments and characterize the flow inside the reactor under specific operating conditions. The change in operating conditions is simulated by replacing the flow matrices of the current state of operation with those associated with the new state. This is based on the assumption that the time for the flow to reach a new steady state as a result of a change in the operational conditions is negligible compared to the total processing time. It should be noted that compartmental modeling provides an approximation of the solution and as such predicts more homogenous distribution for species since each compartment is homogeneous, and does not account for diffusive mass transfer. On the other hand, it makes it possible to take into account hydrodynamics in dynamic analysis of reactor performance and couple the model with optimization solvers. In this work, integration of hydrodynamics with metabolism refers to capturing the effects of dissolved oxygen (DO) concentration, bubbles, and turbulent eddies on the metabolic activities and viability of cells. Biological processes are captured through unstructured modeling. This utilizes a reduced number of reactions to macroscopically capture cellular kinetics.

2.1. Development of CFD Simulations

Computational fluid dynamics (CFD) simulations are developed in ANSYS^® Fluent^® 15.0.7 for the prediction of spatial variations of environmental parameters. Conservation laws of mass, momentum, and energy are usually used to describe a single phase flow, gas or liquid. If the thermodynamic, transport, and chemical properties of a component need to be specified, the field equations may be accompanied by the constitutive equations of state, stress, chemical reactions, etc. The presence of interfacial surface in a multi-phase flow complicates the mathematical formulation of the problem. To derive the field and constitutive equations of a multi-phase flow, such as inside a bioreactor, local characteristics have to be considered. This is not straightforward due to unknown motions of multiple deformable interfaces, variable fluctuations due to turbulence and moving interfaces, and discontinuity of properties at the interface. Obtaining local mean values of flow properties has been shown to be an efficient way to eliminate instantaneous fluctuations. Three averaging methodologies have been developed: Eulerian, Lagrangian, and Boltzmann statistical averaging. In the Eulerian approach, time and space coordinates are independent and other variables are expressed with respect to them. In the Lagrangian averaging methodology, particle coordinates replace spatial coordinates. If the purpose of modeling is studying the group behavior of particles, the Eulerian approach is preferred. However, if the behavior of individual particles is of interest, the Lagrangian description has a clear advantage [45]. Tracking individual bubbles increases the computation. Additionally, it would only improve model predictive power if the extent of the interactions between individual bubbles and the liquid phase could be quantified. These interactions involve growth, breakage, and agglomeration of bubbles and energy dissipation due to bubble rupture. Therefore, in this study gas and liquid phases are treated as continua and Eulerian averaging is used. The Eulerian multiphase model creates sets of momentum and continuity equations for each phase and couples them through exchanging pressure and interphase coefficients [46]. Turbulence of flow is calculated using the k-ε viscosity model, which has been widely used for stirred tanks [47]. It is a robust model that gives reasonably accurate results for a wide range of turbulent flows [48]. A k-ε model consists of two transport equations, one each for the turbulent kinetic energy (k) and the energy dissipation rate (ε). The motion of the impeller is captured using a multiple reference frame (MRF). To implement the MRF model, the geometry is broken up into stationary and moving zones. The MRF model approximates the flow in the moving zone around the impeller by freezing the motion of the moving part in a specific position and observing the instantaneous flow field. To use flow variables of one zone for calculation of fluxes at the boundary of the adjacent zone, a local reference frame transformation is performed at the interface between cell zones. In the absence of large-scale transient effects due to weak impeller–wall interactions, the MRF approach provides a reasonable approximation of the flow [48].

2.2. Development of the Integrated Model

The state of operation is defined by operating volume, impeller rotation speed, and gas sparging flowrate. For every admissible state a separate CFD simulation is developed to obtain the steady-state flow under associated operating conditions. After the solution is converged, the flow information of computational cells is extracted. Compartments are formed through agglomeration of computational cells. Inlet and outlet fluxes, gas volume fraction, dissipation rate of mechanical energy, and gas superficial velocity of compartments are calculated using data obtained from CFD simulations. The surface between two neighbor compartments is composed of faces of computational cells. The mass flowrates reported for these faces are summed to calculate the flowrates between neighbor compartments. The net flow between two neighbors is not necessarily zero, but the net flow for each compartment should be zero in order to maintain mass conservation. Fluent does not report flowrates for faces at the boundaries. Boundary faces are located at the surface of meshing zones. Although very small, these missing values cause mass imbalance and introduce errors. These errors are smoothed out by making the smallest possible changes to the calculated flowrates. It is assumed that the calculated value for mass flowrate from compartment j to i, $F_{ij}$, has error ${\epsilon }_{ij}$. The sum of squared errors is minimized by solving the problem explained by Equations (1) to (3). The second constraint makes sure that zero elements will remain zero, i.e., flowrates between non-neighbor compartments remain zero. The effect of agitation on distribution of cells between compartments is only understood if the sedimentation of cells is captured. The rate of sedimentation is estimated as r²/4 (mm/h) and r is cell radius in μm [49]. Cell radius is estimated assuming spherical shape, density equal to that of water, and average mass of 1.1165 × 10⁻⁶ mg [20].

$$\underset{{\epsilon }_{\mathrm{ij}}}{\mathrm{minimize}}{\displaystyle \sum }_i{\displaystyle \sum }_j{\epsilon }_{\mathrm{ij}}{}^2$$

(1)

subject to:

$${\displaystyle \sum }_i{F}_{ij}-{\epsilon }_{ij}={\displaystyle \sum }_i{F}_{ji}-{\epsilon }_{ji}$$

(2)

$$-F_{ij}\le {\epsilon }_{ij}\le F_{ij}$$

(3)

It has been reported that cells can be damaged when the power input is greater than 22,500 W%·m⁻³ [50,51]. Power input is calculated by multiplying the density of the liquid phase (kg/m³) by the turbulent energy dissipation rate (m²/s³). The turbulent energy dissipation rate is the rate of absorption of kinetic energy that breaks up large eddies. This is then converted to heat by viscous forces [52]. A rate of cell damage of 3.4% min⁻¹ has been reported for cells in high shear regions [51]. The volume fraction of high shear region in compartments is calculated using the data obtained from CFD simulations. Since cells are assumed to be homogenously distributed inside compartments, the volume fraction of the high shear region is equal to the fraction of cells exposed to shear beyond their tolerable threshold. Therefore, the rate of loss of viable cells under operating condition $op$ and in compartment $c$ is calculated using Equation (4):

$$k_{d,shear}{}_{op,c}=0.034.volumefractionofhighshearregio{n}_{op,c}.$$

(4)

Interaction with bubbles has also been reported as one of the sources of cell loss. Cells attach to bubbles, rise with them to the surface, become trapped in the foam layer, and perish. Also, the maximum energy dissipated due to bubble rupture is two or three orders of magnitude higher than the tolerable threshold for cells [51]. No value has been reported in the literature for the rate of cell loss due to interaction with bubbles. One motive for integrated modeling is to capture uncertainty where it occurs. The rate of cell loss due to interaction with bubbles is estimated by assuming an interaction vicinity around bubbles. It is assumed that a fraction of the cells in a compartment that are in the vicinity of bubbles are lost over the average lifespan of a bubble in the compartment. The volume of bubble vicinity is calculated using the reported average bubble diameter of 0.00289 m [53] and assuming a particular value for critical distance from the surface of bubbles. The critical distance is assumed to be equal to the cell radius. The number of bubbles is calculated using gas holdup in the compartment obtained from CFD simulations and the reported value for average bubble diameter. Average bubble lifespan is calculated using the gas holdup and air sparging flowrate. Equation (5) shows the estimated rate of loss of cells due to interactions with bubbles under operating condition $op$ and in compartment $c$:

$$k_{d,bubble}{}_{op,c}=\frac{\ln (1-\frac{volumeoftheinteractionvicinity\times numberofbubble{s}_{op,c}}{volum{e}_c})}{bubblelifespa{n}_{op,c}}.$$

(5)

The integrated model predicts viable cell density (VCD) in compartment $c$ using Equation (6). $\mu$ and ${\mu }_d$ are the metabolic rates of growth and death. They are calculated as functions of metabolites’ concentrations using the metabolic model. $F^{sed}$ is the rate of sedimentation. $c^{\prime }$ and $c^{″}$ are compartments below and above $c$. $N$ is the total number of compartments. It can be seen that the hydrodynamic part of Equation (6) is a linear system of independent ODEs of rank $N$. The discretization of space provided by compartmental modeling allows for the application of proper orthogonal decomposition (POD) for the development of a reduced-order model (ROM) [54]. In order to achieve this, the full rank system is first solved to create time-series snapshots of distribution of cells over compartments under a specific operating condition. Then a set of orthonormal bases are generated through eigen-decomposition [55]. Bases with no significant impact on the solution profile are truncated to obtain the reduced rank model. A system with $N=16$ was used to evaluate the performance of the ROM developed with this methodology. The ROM showed satisfactory performance, while the reduction in the rank of the system was small. The impact of the initial condition used for generating snapshots became apparent when fewer basis functions were used for approximation. For large reactors with a greater number of compartments, however, it is recommended to investigate the application of this methodology.

$$\begin{array}{c}\frac{d{X}_c}{dt}=\left(\mu -{\mu }_d-k_{d,shear}{}_{op,c}+k_{d,bubble}{}_{op,c}-F^{sed}{}_{c^{\prime },c}-\frac{{{\displaystyle \sum }}_{i=1}^N{F}_{i,c}}{{\rho }_l.volum{e}_c}\right)X_c\\ +F^{sed}{}_{c,c^{″}}{X}_{c^{″}}+{\displaystyle \sum }_{i=1}^N\frac{F_{c,i}}{{\rho }_l.volum{e}_c}{X}_i\end{array}$$

(6)

Contrary to cells, metabolites are assumed to be homogenously distributed at all times. This is due to the fact that fast diffusion dominates mass transfer in the small reactor considered for case study. For larger bioreactors, local diffusive mass transfer has lower importance relative to convection, so the flux matrix may be used for calculation of distribution of metabolites [18]. The dissolved oxygen (DO) concentration is also assumed to be homogenously distributed. In the calculation of DO concentration, mass transfer and cellular uptake are considered. The ratio of rates of oxygen uptake to carbon dioxide production by cells has been reported to vary within a narrow range around 1 [56]. So the oxygen uptake rate (OUR) is assumed to be 0.35 pmol·cell⁻¹·h⁻¹, which has been reported for carbon dioxide production [57]. The overall volumetric mass transfer coefficient, k_La, is calculated using Equation (7), in which $U_G{}_{op,q}$ is superficial gas velocity (m/s) for computational cell $q$ under operating condition $op$ [33]. Volumetric mass transfer coefficient is the product of liquid phase mass transfer coefficient; k_L (m/s) and specific interfacial area; a (m²·m⁻³). Mass transfer stops after reaching the saturation concentration at 37 °C. The reported saturation mass fraction of oxygen is 3.43 × 10⁻⁵ for [58]. The unstructured model is assumed to predict metabolite uptake and production rates when the culture is oxygen-saturated. Experimental data show the dependence of these rates on the concentration of dissolved oxygen [59]. The reported data are used to calculate correction factors for metabolites’ uptake and production rates at different concentrations of DO (Figure 2). Uptake and production rates of metabolites predicted by the unstructured model are multiplied by correction factors to take into account the effects of mass transfer mechanism on metabolism of cells.

$${(k_{L}a)}_{op}=\frac{{{\displaystyle \sum }}_{q}0.412U_G{{}_{op,q}}^{0.809}.gasvolumefractio{n}_{op,\mathrm{q}}.volum{e}_q}{{{\displaystyle \sum }}_{q}volum{e}_q}$$

(7)

The behavior of the bio-phase is computed through the incorporation of an unstructured metabolic model and the consideration of the effects of environmental parameters on viable cell density. For the purpose of this study, a metabolic model developed in the literature is adapted to represent cellular growth and death rates as functions of the concentrations of metabolites. As discussed, due to the lumped nature of unstructured models, the estimated values of their parameters have narrow confidence intervals. Moreover, these values lose their meanings when the model is used to predict the dynamics of a different bioreactor system. An important group of parameters are threshold metabolites’ concentrations, which separate growth and death rates into different regimes. Threshold concentrations are determined by observing how viable cell density reacts to concentrations of metabolites. Capturing the dynamic behavior of the system sometimes requires considering multiple phases for cellular growth and death, during which cells react differently to environmental stimuli. Xing et al. [60] assumed the death phase begins after viable cell density declines by 10% from its peak value. It was assumed that cellular growth did not happen during this phase. The integration with hydrodynamics incorporates additional sources of cell loss into the model, which impacts the viable cell density profile. For these reasons, the metabolic model is merely adapted to explain the integration of physical and biological processes. This paper proposes a framework for capturing the interaction of system components and uncertainty where it occurs. The results presented in this paper are only meant to demonstrate the capabilities of the modeling framework.

2.3. Coupling the Model with Nonlinear Solvers

Maximization of bioreactor yield is achieved through manipulation of process parameters based on the optimal operating policy. Reduction of the order of the model through compartmental modeling provides the formulation and modeling environment necessary for coupling the model with nonlinear or mixed-integer nonlinear programming solvers. Optimization algorithms that use reduced models are categorized based on their level of dependence on the original detailed model for the calculation of gradients [61]. In the proposed formulation, reduced models are developed before calling the optimization algorithm since minimal communication between the algorithm and the original model is necessary to reduce the computational complexity. Advanced nonlinear optimization algorithms are capable of handling numerous decision variables and constraints. However, the inherent dynamics of this problem makes it challenging. The classical approach to dynamic optimization problems takes advantage of Pontryagin’s maximum principle and maximizes the control Hamiltonian over the set of all admissible controls [62]. The application of this approach becomes difficult for larger systems with state constraints [63], so direct approaches based on parameterization of variables have been preferred [64,65]. The method of collocation has been proposed for parameterization of variables [66,67]. The stiffness of the system of ordinary differential equations (ODE) is evaluated for different initial values. The stiffness ratio is calculated using the eigenvalues of the Jacobian matrix at different points in time [68]. It is observed that the order of magnitude of the stiffness ratio varies between 6 and 28 throughout the integration. Therefore, an approach based on the discretization of time using fixed step sizes does not provide a good approximation of the solution unless the step size is very small, i.e., less than 10⁻⁹ h. Instead, integration is carried out using appropriate solvers for stiff, nonlinear ODEs and the problem is formulated for the application of the interior point method [69]. Equations (8)–(11) represent the solution to the mathematical optimization problem. $T_i$ and $C_i$ are the time and composition of the ith feeding. In addition to initial nutrient concentrations, schedule, and composition of feeding, it also finds optimal criteria for setting aeration and agitation rates. Aeration is stopped or started based on the DO level. For the adjustment of impeller rotation speed, a measure of homogeneity is defined based on the relative standard deviation (RSD) of distribution of cells over compartments.

$$\underset{C_0,T_i,C_i,{\mathrm{RSD}}_{\mathrm{cri}},{\mathrm{DO}}_{\mathrm{cri}}}{\mathrm{maximize}}{{\displaystyle \int }}_0^{t_f}\mathrm{biomassdt}$$

(8)

Subject to:

Process Model: Equation (6)

Control Bounds:

$$0<\dots <T_{i-1}<T_i<T_{i+1}<\dots <t_f$$

(9)

$$C_i{}_{\min }\le C_i\le C_i{}_{\max }$$

(10)

$$0\%\le {\mathrm{RSD}}_{\mathrm{cri}},{\mathrm{DO}}_{\mathrm{cri}}\le 100\%.$$

(11)

3. Case Study

A 3 L bioreactor with Rushton impeller and sparger is considered. The operation lasts for two weeks. The dimensions of the reactor are shown in Figure 3. It is operated at liquid fill levels of 130, 155, 180, and 205 mm. Feeding is simulated by change of liquid fill level. So, overall, three steps of feeding are allowed. For the discretization of space into computational cells, the space is divided into several meshing zones. This helps Fluent to generate hexahedral cells where possible. It also reduces the total number of computational cells and improves the convergence of the solution. Meshing for the highest operating volume results in 363,175 computational cells, 947,988 faces, and 226,971 nodes. Agitation rates of 150, 225, and 300 RPM and two options for aeration are considered: sparging air at 0.01 vessel volume per minute (vvm), and no aeration.

Table 1. Physical properties calculated from CFD simulations.

No Aeration			Aerated System
Fill Level (mm)	Impeller Rotation Speed (RPM)	Power Input (W%·m−3)	Fill Level (mm)	Impeller Rotation Speed (RPM)	Power Input (W%·m−3)	Volumetric Mass Transfer (h−1)
130	150	3.8	130	150	3.5	12.8
	225	11.0		225	11.1	16.6
	300	24.5		300	26.5	14.9
155	150	3.2	155	150	3.1	10.9
	225	9.5		225	9.7	11.3
	300	21.6		300	23.6	16.1
180	150	2.7	180	150	2.4	14.4
	225	8.4		225	7.9	14.8
	300	19.2		300	19.6	20.4
205	150	2.3	205	150	2.7	11.9
	225	7.0		225	6.5	9.6
	300	16.1		300	15.5	13.7

shows the physical properties calculated from CFD simulations. The reactor is divided into compartments as shown in Figure 4. The data exported from CFD simulations are used to calculate flow matrices for all states of operation. The unstructured model developed by Xing et al. [60] is adapted to predict the behavior of the bio-phase. The metabolic model captures the effects of concentrations of glucose, glutamine, lactate, and ammonia on cellular rates of growth and death in a CHO culture. The rate of utilization of glutamine for essential metabolic functions, i.e., maintenance, is calculated from Equation (12). The values of the model parameters are shown in

Table 2. Unstructured model parameters.

Parameter	Value	Unit
$m_{Glc}$	$6.92\times {10}^{-11}$	$\mathrm{mmol}·{\mathrm{cell}}^{-1}·{\mathrm{h}}^{-1}$
$a_1$	$3.2\times {10}^{-12}$	$\mathrm{mmol}·{\mathrm{cell}}^{-1}·{\mathrm{h}}^{-1}$
$a_2$	$2.1$	$\mathrm{mM}$
${\mu }_{max}$	$0.029$	${\mathrm{h}}^{-1}$
${\mu }_d{}_{max}$	$0.016$	${\mathrm{h}}^{-1}$
$K_{Glc}$	$0.084$	$\mathrm{mM}$
$K_{Gln}$	$0.047$	$\mathrm{mM}$
$K{I}_{Lac}$	$43$	$\mathrm{mM}$
$K{I}_{Amm}$	$6.51$	$\mathrm{mM}$
$K{D}_{Lac}$	$45.8$	$\mathrm{mM}$
$K{D}_{Amm}$	$6.51$	$\mathrm{mM}$
$d_{Gln}$	$7.2\times {10}^{-3}$	${\mathrm{h}}^{-1}$
$Y_{X/Glc}$	$1.69\times {10}^8$	$\mathrm{cell}·{\mathrm{mmol}}^{-1}$
$Y_{X/Gln}$	$9.74\times {10}^8$	$\mathrm{cell}·{\mathrm{mmol}}^{-1}$
$Y_{Lac/Glc}$	$1.23$	$\mathrm{mmol}·{\mathrm{mmol}}^{-1}$
$Y_{Amm/Gln}$	$0.67$	$\mathrm{mmol}·{\mathrm{mmol}}^{-1}$
$D{O}_{eq}$	$1.0699$	$\mathrm{mM}$
$OUR$	$3.5\times {10}^{-10}$	$\mathrm{mmol}·{\mathrm{cell}}^{-1}·{\mathrm{h}}^{-1}$

.

$$m_{Gln}=\frac{a_1[Gln]}{a_2+[Gln]}$$

(12)

Cellular rates of growth and death are calculated by knowing the Monod constants ($K_{metabolite}$) and using Equations (13) and (14), respectively. Growth and death rates are then used in Equation (6) to calculate viable the cell densities in compartments:

$$\mu ={\mu }_{max}\frac{\left[Glc\right]}{K_{Glc}+\left[Glc\right]}\frac{\left[Gln\right]}{K_{Gln}+\left[Gln\right]}\frac{K{I}_{Lac}}{K{I}_{Lac}+\left[Lac\right]}\frac{K{I}_{Amm}}{K{I}_{Amm}+\left[Amm\right]}$$

(13)

$${\mu }_d={\mu }_d{}_{max}\frac{\left[Lac\right]}{K{D}_{Lac}+\left[Lac\right]}\frac{\left[Amm\right]}{K{D}_{Amm}+\left[Amm\right]}.$$

(14)

The estimation of overall volumetric mass transfer coefficient under present operation conditions is used in Equation (15) to calculate the concentration of DO, where X is the viable cell density.

$$\frac{d\left[DO\right]}{\mathrm{dt}}={(k_{L}a)}_{op}\left(D{O}_{eq}-\left[DO\right]\right)-X.OUR$$

(15)

Concentrations of metabolites are calculated from Equations (16)–(19) by knowing the values of yield parameters ($Y$’s). The impact of DO on uptake and production rates ($D{O}_{metabolite}$) is estimated through spline interpolation of the experimental data shown in Figure 2. The model also captures the chemical degradation of glutamine.

$$\frac{d\left[Glc\right]}{dt}=-D{O}_{Glc}(\frac{\mu -{\mu }_d}{Y_{X/Glc}}+m_{Glc})X$$

(16)

$$\frac{d\left[Gln\right]}{dt}=-D{O}_{Gln}\left(\frac{\mu -{\mu }_d}{Y_{X/Gln}}+m_{Gln}\right)X-d_{Gln}\left[Gln\right]$$

(17)

$$\frac{d\left[Lac\right]}{dt}=D{O}_{Lac}.Y_{Lac/Glc}(\frac{\mu -{\mu }_d}{Y_{X/Glc}}+m_{Glc})X$$

(18)

$$\frac{d\left[Amm\right]}{dt}=D{O}_{Amm}.Y_{Amm/Gln}\frac{\mu -{\mu }_d}{Y_{X/Gln}}X$$

(19)

The definition of states of operation creates a finite set of admissible actions for every point in time. Depending on the present operational conditions, it is possible to increase or decrease impeller rotation speed, stop or start air sparging, and feed or not feed. It is also possible to continue under the current conditions. Changing the state of operation is limited to once every 2 h. The concentrations of glucose and glutamine, initially and in feed, are constrained to be under 100 and 10 mM, respectively. The solid lines in Figure 5 show the system operated under the policy obtained from solving the optimization problem described in Equations (8)–(11) (system 1). Although the existence of multiple local optima cannot be ruled out, the solution shows improvement in the yield of operation through manipulation of feeding schedule and composition, agitation, and aeration rates. Nutrient concentrations stay at their upper bounds due to the limited number of feeding steps. The criterion for aeration is 47% of saturated DO concentration, i.e., aeration starts if DO concentration falls below this value and stops otherwise. The obtained criterion for agitation is 0.02%. The impeller rotation speed is increased if the RSD of distribution of cells over compartments is greater than the agitation criterion. In the opposite case, where RSD has a smaller value, the lower agitation rate is selected for the next 2 h. For comparison, the dashed lines in Figure 5 show the system operating with the same initial and feed compositions but a uniform feeding schedule (system 2). System 2 is consistently aerated and agitated at 300 RPM, i.e., the criteria for aeration and agitation are 100% DO saturation and RSD of 0%, respectively. Schedule of feeding should be determined with consideration of capacity of the reactor and duration of operation. Early addition of feed, despite increasing cellular population, causes accumulation of lactate and ammonia in the system, which further inhibits growth. Late feeding, on the other hand, results in poor utilization of nutrients and reactor capacity. Solving the optimization problem results in feeding times of 215, 265, and 300 h after the start of cultivation. Even with a limited number of feeding steps, through the manipulation of the feeding schedule, the cells in system 1 are provided with enough nutrients to maintain growth throughout the operation. Inclusion of mass transfer mechanism in the model leads to an improvement in aeration. The DO concentration drops quickly when aeration stops because of fast consumption by cells. Stopping and starting aeration according to the near-optimal policy prevents loss of viable cells due to unnecessary aeration while guaranteeing that DO is not depleted. The impeller rotation speed switches between 150 and 225 RPM for most of the operation to maintain RSD of distribution of cells at 0.02%. The three sharp peaks in RSD values show disturbances in spatial homogeneity of cells caused by feeding. The declining trend in net growth that happens toward the end of the process is due to the fact that the cellular growth rate, which is reduced by the inhibitory effects of metabolites, cannot compete with the cell loss due to the effects of hydrodynamics. The results demonstrate that integrated modeling is able to capture the behavior of the system using a mechanistic understanding of the reactor without the need for unnecessary assumptions.

4. Conclusions and Future Directions

The implementation of mechanistic models in the biotechnology industry has been hindered by a lack of universality of cell culture models. Integrated modeling, followed by experimental design and parameter estimation, can lead to quantification of the extent of effects of mechanical shear and bubble interactions on the viability of cells. It also gives better estimations of the cellular rates of growth and death as functions of metabolite concentrations. A model that represents the system well can be linked to proper optimization algorithms to recommend low-cost improvements for the operation of a bioreactor. Furthermore, computationally feasible unit operation models facilitate the integration of control, scheduling, and planning as a leading step toward integrated decision-making.

The incorporation of a buffer system into the model results in better representation of the system. It makes it possible to calculate pH and capture its effects on biological processes. This is achieved at the cost of increasing the nonlinearity and rank of the ODE system. The discrete space of admissible operating conditions can be expanded through finer discretization of process parameters. It can also be replaced by a continuous space through multi-scale surrogate modeling. Data obtained from CFD simulations for computational cells can be used for the development of surrogate models for hydrodynamics [70]. The main challenge in this area is to devise an efficient algorithm to explore the large sampling space. To select a subset of computational cells, Zhao et al. proposed a sampling method based on Latin hypercube designs (LHDs) [71]. After decomposing the complete data into disjoint equally spaced blocks, a subsample is obtained by collecting blocks according to a randomly generated LHD. This method is called LHD-based block bootstrap. It takes into account the spatial dependency and therefore improves the accuracy of estimations. Integration with hydrodynamics introduces new parameters to the cell culture model, i.e., sedimentation rate, tolerable shear threshold, rate of cell damage due to shear, bubble radius, interaction distance with bubbles, and mass transfer coefficients. The uncertainties in the values of these parameters can be included in the dynamic analysis of the operation. The sensitivity of the solution profile to model parameters can be determined and the uncertainty in sensitive model parameters can be considered while solving for the optimal operating policy. Although the proposed approach cannot guarantee that the global solution is obtained, since a local optimization algorithm is utilized, the convergence to the global optimal solution can be improved using initialization strategies in the interior point method.