ABSTRACT

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.

Keywords

Cell Population Balance, Cell Growth, Substrate Consumption, Central Upwind Scheme, Equal and Unequal Partitioning of Cells

1. 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.

2. 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 is a function that represents the state of entire population with the amount of biomass x at time t between x and. The zero moment and the first moment respectively describe the total number of cells per unit bio-volume and concentration. They are defined as [13]

(1)

(2)

The number density function is defined as

(3)

The computational domain is taken to be, where and represent the minimum and maximum values for the amount of biomass components of daughter cells, respectively. Let is smallest value of the mass components of parent cells, where Moreover, the minimum and maximum values of all mass components are considered to be and. Here, indicates the growth rate, represents the nutrient consumption process, represents division rate of a cell and D denotes the dilution rate. The partitioning function illustrates the cellular division process of a parent into daughter cell. It gives the probability that a parent cell having state x gives birth to a daughter cell with state y. The following relation satisfies the normalization condition.

(4)

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 must be zero for all states of daughter cells which are larger than the states of parent cells,

(5)

Furthermore, the probability of a dividing cell with physiological state y to produce a daughter cell of state must be equal to the probability of producing a daughter cell of state, i.e.

(6)

Due to the above-mentioned hypothesis and process description, the dynamics of the state distribution function is described by the general cell population balance equation [13]

(7)

where

(8)

with initial conditions

(9)

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

(10)

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

(11)

Under this assumption, Equation (8) reduces to [13]

(12)

Thus, the integral term disappears. Since the maximum and minimum values of the dividing cell are

and, respectively, the corresponding values for each newborn cell must be and

, respectively. Thus, the birth term in Equation (12) exists only for the domain. The substrate consumption in time is described by an ordinary integro-differential equation of the form [13]

(13)

with initial conditions

. (14)

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, denotes the consumption rate and is the saturation concentration. In this model, coupling is the only reason for nonlinearity. For a constant abiotic environment, the problem becomes linear.

3. 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 are partitions of the given interval. Moreover, denotes the constant width of each grid interval, represent the cell centers, and indicate the cell boundaries. We refer,

. (15)

Moreover,

(16)

Let for. For every, the averaged values of the conservative variables and are defined as

. (17)

After integrating Equation (7) over the interval, we get

. (18)

The numerical fluxes are given as [19]

(19)

Here, N⁺ and N⁻ shows the piecewise linear reconstruction for N, namely

(20)

Here, are the first order approximations of 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, for example

, (21)

where and

(22)

Here, denotes central differencing and MM is the min-mod nonlinear limiter. Further, the local one sided speed at is given as [19]

. (23)

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

, (24)

. (25)

where, represent right hand sight of Equation (18) and represents solution at next time step and is time steps.

4. 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 is expressed as [13]

, (26)

where

(27)

The initial condition is given as

. (28)

Here, , , and the boundary conditions are defined as

(29)

The above cell population balance model is coupled with the mass balance of substrates [13]

(30)

with initial condition

. (31)

5. 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 which is concealed by cells

for and. We take

. (32)

. (33)

At any time t, the cell averaged values for conserved variables are given as

, (34)

and

. (35)

The piecewise linear interpolation is described as [19] [20]

, (36)

where, and and are the approximations of the x and y-derivatives of N at the cell centres To calculate these partial derivatives, we use the minmod limiter (MM) with,

(37)

(38)

On integrating Equation (26) over the control volume, the two dimensional scheme can be expressed as [19] [20]

(39)

Here

(40)

(41)

where, we have

, (42)

(43)

The local speeds are computed as

(44)

(45)

6. 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.

6.1. 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.

•  The dilution rate D = 0.

•  Due to constant abiotic environment, Equation (13) for substrate consumption is not considered.

•  Growth rate functions are consider as:

1) Constant growth rate:

(46)

2) Linear growth rate:

(47)

3) Quadratic growth rate:

(48)

where, m₀ 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]

. (49)

The division function is given as [13]

. (50)

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.

6.2. 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]

1) Test problem 1: Constant growth

2) Test problem 1: Linear growth

3) Test problem 1: Quadratic growth

(51)

where represents the maximum concentration of the substrates and K is the Michaelis constant describing the substrate concentration due to which the reaction rate is equal to half of.

•  The cell consumption rate is given as [13] -[15]

. (52)

Here y is constant mass. Take doubling time is and the simulation time is 20h. Figures 13-16 display the numerical results for equal and unequal partitioning. The numerical results of the mass concentrations are similar to those available in the literature [13] . In equal partitioning, the number density function represents the constantly oscillatory behavior which disappears in unequal partitioning. From these results it can be noticed that central scheme has comparable accuracy and the first order scheme is diffusive.

6.3. Test Problems 3: Bivariate Case

Here, the following assumptions in two dimensions are considered.

•  The growth rates of two single cells are shown by the following equations.

. (53)

•  Here, equal partitioning function is only considered.

•  The division rate Γ is given by

. (54)

•  The state distribution function for first variable is define as

. (55)

•  The number density function with respect to the first variable x is

. (56)

Figure 17 shows the comparison of mass concentration and middle point of number density function obtained by central upwind and first order schemes. From this figure it is evident that central schemes confine the periodic behavior of the mid-point of the number density function. The reason, however, is that partitioning density function and growth rates are the same as used in Test problem 1 and mass Equation (13) is not considered because growth rates are independent of substrate concentration s. Due to numerical diffusion in the first order scheme, it is unable to maintain the period oscillations. Figure 18 shows the 2D plots of number density function and state distribution function obtained from the central scheme. By integrating the number density function and state distribution function with respect to second variable y, we can obtain these results (c.f. Equations (55) and (56)). Thus, the proposed scheme is also accurate in the bivariate case.

7. Conclusion

The central-upwind scheme was applied to solve the seingle-variate and bivariate cell population balance models considering equal and unequal partitioning of cellular materials. The proposed scheme can be applied for any growth rate, cell division rate and partitioning probability density function due to its general nature. The scheme gives minimal numerical diffusion and avoids oscillations that can be normally observed in the second and

higher order schemes. The numerical test problems verify the usefulness of the suggested-scheme which is computationally efficient, accurate, and easily applicable.

References