Central Upwind Scheme for Solving Multivariate Cell Population Balance Models

Microbial cultures are comprised of heterogeneous cells that differ according to their size and intracellular concentrations of DNA, proteins and other constituents. Because of the included level of details, multi-variable cell population balance models (PBMs) offer the most general way to describe the complicated phenomena associated with cell growth, substrate consumption and product formation. For that reason, solving and understanding of such models are essential to predict and control cell growth in the processes of biotechnological interest. Such models typically consist of a partial integro-differential equation for describing cell growth and an ordinary integro-differential equation for representing substrate consumption. However, the involved mathematical complexities make their numerical solutions challenging for the given numerical scheme. In this article, the central upwind scheme is applied to solve the single-variate and bivariate cell population balance models considering equal and unequal partitioning of cellular materials. The validity of the developed algorithms is verified through several case studies. It was found that the suggested scheme is more reliable and effective.


Introduction
In mammalian cell culture, individual cells exhibit heterogeneity due to differences in their cellular metabolism and cell-cycle dynamics [1].In the step-by-step cell cycle action each cell of the population entity grows to a certain size (approximately double to its original size) and then divides into two identical daughter cells.Basically cell division is an exponential process, the two daughter cells further divide into four daughter cells, then four into eight and so on [2].The environmental factors such as oxygen, pH, temperature and substrate concentration greatly affect the cell growth rate.At any point in time t, in heterogeneous population different cells exist at different stages of the cell cycle.The different cells could be differentiated according to their size, DNA and RNA contents, protein contents and other inter-cellular properties.Thus, the developments of accurate mathematical models are needed that account for heterogeneities present at the single-cell level.
Various mathematical models exist in the literature for describing dynamics of biological systems.These models consider the cell population as a "continuum" or lumped biophase, thus assuming that it behaves as a homogeneous entity.Cell population balance models (PBMs) are the only models that take into account the fact that cell properties are distributed among the cells of a population, such as protein content, DNA content, etc.Moreover, they have capability to effectively describe the internal chemical structure of the single cell by incorporating intracellular chemical reactions involving multiple chemical species.Therefore, these models provide the most accurate way of describing the complicated phenomena associated with cell growth, nutrient uptake and/or product formation in microbial populations.They typically consist of multidimensional integrto-partial differential equations for describing the dynamics of the state distribution function, nonlinearly coupled with integro-ordinary differential equations accounting for substrate consumption and (or) product formation.Because of their obvious mathematical complexity, the numerical solutions of such models are challenging for the given numerical scheme.In mid 1960s, these models were used for the first time to simulate cell dynamics [3] [4].Cell PBMs are either single-variable or multi-variable, depending on numbers of variables involved in the system.Multivariable models explain the complicated phenomena of cell development, product formation and substrate consumption as compared to the single-variable models.Despite of these details, bioreactor PBMs are associated with many difficulties.Firstly, for majority cell system the growth rate function, the cell division rate and the partitioning probability density function are not known, that can be overcome by flow cytometry analysis.The on-line flow cytometry, when combined with suitable cell population models, develops the computerbased system that provides control of cellular distribution.Secondly, a great computational cost is required to numerically solve multi-variable cell PBMs.
Initially, Hulburt and Katz [5] introduced population balance models in chemical engineering, and these models were later on developed by Randolph and Lason [6].By the time, their applied nature recognized them to be the vital part of research.The large computational time and unavailability of exact solution were the primary obstacles in the development of these balances.The analytical solutions are possible in simplified situations only.Thus, researchers started to find its numerical solution instead of the exact solution.A large number of numerical techniques are present for approximating population balance equations (PBEs) in chemical and biological engineering, such as the weighted residual method [7], the Monte Carto simulation [8], the method of moments [9] [10], the finite difference method [11]- [13], the spectral methods [14], the finite element method [15] [16], and the high resolution finite volume scheme [17] [18], etc.
In this article, the high-resolution non-oscillatory central-upwind scheme is applied to solve the single-variate and bivariate cell population balance models considering the equal and unequal partitioning of mother cells [19] [20].The basic idea of such schemes is that they use information of local propagation speeds and the approximated solution is obtained as cell averages.Further, such schemes have an upwind nature, as they take care of directions of wave propagation by measuring the one-sided local speeds.Several case studies are carried out and the results of central-upwind scheme are compared with those obtained from the first order upwind scheme.

Single-Variate Cell Population Balance Model
A single-variate cell population balance model that expresses the ratio of nutrients is considered in this section.Let ( ) , N x t is a function that represents the state of entire population with the amount of biomass x at time t between x and d x x + .The zero moment ( ) 0 N t and the first moment ( ) N t respectively describe the total number of cells per unit bio-volume and concentration.They are defined as [13] ( ) ( ) The number density function is defined as The computational domain is taken to be This means that the amount of biochemical components is conserved at cell partitioning.Particularly, daughter cells cannot have larger amounts of biochemical components than the parent cells.Thus, the probability density function ( ) , , P x y s must be zero for all states of daughter cells which are larger than the states of parent cells, ( ) Furthermore, the probability of a dividing cell with physiological state y to produce a daughter cell of state x must be equal to the probability of producing a daughter cell of state y x − , i.e.

P x y s P y x y s
Due to the above-mentioned hypothesis and process description, the dynamics of the state distribution function ( ) , N x t is described by the general cell population balance equation [13] ( ) where

Q x t x s N x t DN x t y s N y t P x y s x
with initial conditions Let us define boundary b of the state as a point where at least one quantity of biochemical components obtains either its minimum value or its maximum value.Then, the boundary conditions are given as Equation ( 8) has three terms, the first term shows accumulation, the second term accounts for the growth into larger cells, and the third term is the source term.In the source term (c.f.Equation ( 8)), the first term represents loss of cells due to the partitioning of cells which leads to the birth of daughter cells.The second term is the dilution rate.The integral birth term is multiplied by the factor 2 for representing the division of a parent cell into two daughter cells.The equal partitioning, where each mass of a mother cell is equally divided into two daughter cells, can be mathematically described by replacing the partition probability density function by a Dirac delta function ( ) , , . 2 Under this assumption, Equation ( 8) reduces to [13] ( ) Thus, the integral term disappears.Since the maximum and minimum values of the dividing cell are , respectively.Thus, the birth term in Equation ( 12) exists only for the domain ,min ,max 2, 2  .The substrate consumption in time is described by an ordinary integro-differential equation of the form [13] ( ) The first term in Equation ( 13) denotes the inlet minus outlet rates from the reactor, and the integral term describes the rate of loss of substrate due to cell growth.Here, ( ) , q x s denotes the consumption rate and f s is the saturation concentration.In this model, coupling is the only reason for nonlinearity.For a constant abiotic environment, the problem becomes linear.

Numerical Scheme for Single-Variate Cell PBM
Here, the semi discrete central-upwind scheme is derived to numerically approximate the PBE model in Equation (7) [19] [20].Before applying the scheme, we need to discretize the computational domain.Let n be the number of discretization points and Moreover, x ∆ denotes the constant width of each grid interval, i x represent the cell centers, and cate the cell boundaries.We refer, Moreover, ( ) 2 and . 1 After integrating Equation ( 7) over the interval The numerical fluxes are given as [19] Here, N + and N − shows the piecewise linear reconstruction N  for N, namely Here, x i N are the first order approximations of ( ) N x t and are calculated using a nonlinear limiter that confirms the non-oscillatory nature of reconstruction [19] [20].The computation of these slopes is given by family of discrete derivatives parameterized with where Here, ∆ denotes central differencing and MM is the min-mod nonlinear limiter.Further, the local one sided speed at is given as [19] ( ) ( ) The second order TVD Runge-Kutta scheme is used to solve Equation ( 18) to achieve second order accuracy in time.That upgrades N in the following two stages ( ) ( ) ( ) where, ( )

Bivariate Cell Population Balance Model
This section considers the bivariate cell population balance model for equal partitioning.Here, two property coordinates are used.The evolution of ( ) , , g x s N x y t g x s N x y t N x y t Q x y t t x y where The initial condition is given as Here, , N x y ∈  , and the boundary conditions are defined as The above cell population balance model is coupled with the mass balance of substrates [13] ( )

Numerical Scheme for Bivariate Cell PBM
Let n x and n y represents large integers in the x and y-directions, respectively.We assume a Cartesian grid with a rectangular domain ,min max ,min max , , and 1 y j n ≤ ≤ .We take ( ) x y x y At any time t, the cell averaged values ( ) , i j N t for conserved variables are given as ( ) The piecewise linear interpolation is described as [19] [20] ( , , On integrating Equation (26) over the control volume i j C , the two dimensional scheme can be expressed as [19] [20] Here , 2 2 , 2 2 where, we have , 2 2 , . 2 2 The local speeds are computed as , .

Numerical Case Studies
In this section, a few single-variate and bivariate case studies are considered.The suggested numerical schemes are qualitatively and quantitatively analyzed.

Test Problem 1: Single-Variate Case
The following presumptions are considered for this problem • The Gaussian division probability density function is taken as the initial distribution.
• Due to constant abiotic environment, Equation ( 13) for substrate consumption is not considered.
• Growth rate functions are consider as: 1) Constant growth rate: 2) Linear growth rate: 3) Quadratic growth rate: ( ) where, m 0 is average mass at time t and U d is the average doubling time.In the unequal partitioning, the number density function does not show a periodic behavior.In that case, the partitioning mechanism is calculated by beta distribution [13] [21] ( ) ( ) The division function is given as [13] ( Here, f is division probability function which is a truncated Gaussian distribution and it is assumed that a balanced growth state with time independent mass is reached [21].
In the case of constant, linear and quadratic growth rates doubling time for equal and unequal partitioning are taken to be U d = 5h, U d = 2h and U d = 5h and simulation times are 30h, 20h and 30h, respectively.Figures 1-4 show the constant growths of state distribution function and number density function in two and three dimensions obtained by the comparison of upwind central scheme with first order backward difference method.Figures 5-8 describe the linear growth obtained by the same comparison of upwind central scheme with first order backward difference method.Similarly, Figures 9-12 show quadratic growth.The values of parametric are given in Table 1.From the figures, it is clear that in both cases of equal partitioning and unequal partitioning, solution reaches the fastest for quadratic growth and slowest for the linear growth and it is in between in constant growth.It is observed from figure that when time is double the mass concentration become double and also it is obvious from above comparison that first order scheme is more diffusive while central scheme shows accuracy.
The results verify that central scheme can be used to capture such profiles more efficiently and accurately.

Test Problem 2: Single-Variate Case
In this problem, the following assumptions are taken into account.Other parametric values are listed in Table 2.
In this case growth rates are also depending on the substrate concentration.Therefore, Equations ( 7) and ( 13) are solved together.
• The growth rate is taken as [13]- [15]                      higher order schemes.The numerical test problems verify the usefulness of the suggested-scheme which is computationally efficient, accurate, and easily applicable.
of the x and y-derivatives of N at the cell centres ( )

Figure 2 .
Figure 2. Test problem 1 (single-variate case): State distribution functions for constant growth rate.

Figure 3 .
Figure 3. Test problem 1 (single-variate case): Number density functions for constant growth rate.

Figure 4 .
Figure 4. Test problem 1 (single-variate case): Number density functions for constant growth rate in 3-D.

Figure 6 .
Figure 6.Test problem 1 (single-variate case): State distribution functions for linear growth rate.

Figure 7 .
Figure 7. Test problem 1 (single-variate case): Number density functions for linear growth rate.

Figure 8 .
Figure 8. Test problem 1 (single-variate case): Number density functions for linear growth rate in 3-D.
the averaged values of the conservative variables