Numerical solution of multi-variable cell population balance models: I. Finite difference methods

Computers and Chemical Engineering 25 (2001) 1411– 1440 www.elsevier.com/locate/compchemeng

Numerical solution of multi-variable cell population balance models: I. Finite difference methods Nikolaos V. Mantzaris a,b, Prodromos Daoutidis a,b, Friedrich Srienc a,b,* a b

Department of Chemical Engineering and Materials Science, Uni6ersity of Minnesota, Minneapolis, MN 55455, USA Biological Process Technology Institute, Uni6ersity of Minnesota, 1479 Gortner A6e. c 240, St. Paul, MN 55108, USA Received 10 August 2000; received in revised form 23 April 2001; accepted 1 May 2001

Abstract Multi-variable cell population balance models represent the most accurate and general way of describing the complicated phenomena associated with cell growth, substrate consumption and product formation due to the level of detail included in them. Therefore, the ability to solve and understand such models is of fundamental importance in predicting and/or controlling cell growth in processes of biotechnological interest. However, due to the fact that such models typically consist of first-order, partial integro-differential equations coupled in a nonlinear fashion with ordinary integro-differential equations, their solution poses a serious challenge. In this work, we have developed several finite difference algorithms for the solution of the problem in its most general formulation (i.e. for any set of single-cell physiological state functions). The validity of the developed algorithms was verified by comparing their results with those of three specific test problems for which several solution characteristics are known. Moreover, the numerical schemes were compared to each other with respect to their key numerical features, such as stability, accuracy and computational speed. Solutions of the cell population balance model with up to three state variables were obtained using a Pentium II 450 MHz PC in tractable CPU times. © 2001 Elsevier Science Ltd. All rights reserved. Keywords: Cell population balance; Cell growth; Substrate consumption; Numerical solution; Finite differences; Numerical stability

1. Introduction It is an obvious fact that a cell population consists of individual cells. Each cell of the population undergoes the so-called cell cycle, during which it grows, and after a certain point, it divides and partitions its cellular material into two daughter cells, each of which enters its own cell cycle. As a result, at any point in time, the cells of the population exist at different stages of the cell cycle. Thus, each cell contains different quantities of proteins, DNA, RNA and other cell properties. Consequently, since cell properties are really distributed among the cells of the population, a cell population is a heterogeneous system. The vast majority of mathematical models that have been suggested in the literature to describe the dynamics of biological systems treat the cell population as a * Corresponding author. Tel.: + 1-612-6249776; fax: +1-6126251700. E-mail address: [email protected] (F. Srienc).

‘continuum’ or a lumped biophase, thus assuming that it behaves as a homogeneous entity. The only mathematical models that ‘respect’ the heterogeneous and distributed nature of cell growth by recognizing the simple fact that a cell population consists of individual cells are the cell population balance models, also known as corpuscular models (Roels, 1983). Moreover, the mathematical formulation of such models naturally allows the incorporation of information about stage-tostage transition and partitioning of cellular material upon cell division. This type of information is simply not present in continuum models. Hence, due to the level of detail built into their mathematical formulation, cell population balance models represent the most accurate way of describing the complicated phenomena associated with cell growth, nutrient uptake and product formation. Furthermore, contrary to continuum models, which can predict only average population properties, cell population balance models are able to predict entire cell property distributions. Thus, when

0098-1354/01/$ - see front matter © 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 7 0 9 - 8

1412

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

used in conjunction with modern control approaches for partial differential equations (e.g. Christofides & Daoutidis, 1997), they can address the control of a wider range of objectives of biotechnological interest. Depending on the number of cell properties or variables that are used to describe the state of the single cell at any point in time, cell population balance models are classified in single- or multi-variable models. Multivariable models can effectively describe the internal chemical structure of the single cell by accounting for intracellular chemical reactions, and are obviously more detailed than single-variable models, which use the quantity of only one cell property in order to distinguish a cell of the population from its fellows. Moreover, depending on the number of cell-cycle stages that are included in the cell population balance formulation, these models are assorted in single- and multi-staged models. Finally, if the property or properties that are used to describe the state of the single cell obey the mass conservation law, then the model is referred to as mass structured, whereas if age is used to differentiate each cell from other cells of the population, then the model is referred to as age structured. Despite the generality, accuracy, predictive power of cell population balance models and the fact that they have been formulated since the mid-1960s (Eakman, Fredrickson, & Tsuchiya, 1966; Tsuchiya, Fredrickson, & Aris, 1966; Fredrickson, Ramkrishna, & Tsuchiya, 1967), the use of these models for modeling, control, design and optimization of bioprocesses has been sparse. This can be attributed to two reasons. First, in order to be applied for any purpose, cell population balance models require information at the single-cell level. In particular, the so-called intrinsic physiological state functions (single-cell growth rates, single-cell stage-to-stage transition rates and partitioning function) need to be known. The experimental determination of these functions is very hard, mainly due to the fact that it requires measurements at the single-cell level. Second, contrary to continuum models that are mathematically formulated as a set of nonlinear ordinary differential equations, cell population balance models typically consist of first-order partial integrodifferential equations describing the dynamics of cell property distributions, coupled in a nonlinear fashion with ordinary integro-differential equations accounting for substrate consumption and/or product formation. As a result, the development of numerical algorithms that enable the accurate approximation of the solution of cell population balance models is a challenging task. The evolution of flow cytometry (Srienc, 1993) has contributed significantly in providing a basis for obtaining information at the single cell level. The analysis of such information with inverse population balance modeling techniques (Ramkrishna, 1994) has enabled, in some cases, the determination of the intrinsic physio-

logical state functions (Block, Eitzman, Wangensteen, & Srienc, 1990; Kromenaker & Srienc, 1991; Srienc & Dien, 1992; Sweeney, Srienc, & Fredrickson, 1994; Hatzis, Fredrickson, & Srienc, 1997; Natarajan & Srienc, 1999). Despite this progress, solutions of cell population balance models are still hard to obtain. Trucco (1965a,b) has analytically solved the age-structured cell population balance model using the method of characteristics. We further note the work of Hjortso and Bailey (1983) who studied transient responses of budding yeast populations by using age-structured population models. The same type of models was employed by Hjortso and Nielsen (1994, 1995) in order to investigate microbial oscillations. Moreover, several reports have been presented on the numerical solution of different age-structured cell population balances, which cannot be solved analytically (Kostova, 1988, 1990; Kostova & Marcheva, 1989; Kim, 1996; Kim & Park, 1995a,b; Chiu, 1990). However, age-structured models are of limited practical value due to the fact that age is very difficult to measure experimentally in microbial populations. On the other hand, some properties of cells, such as volume, total protein content, and DNA content are not difficult to measure, even at the single-cell level. Therefore, the use of such properties in the formulation of mass-structured models is more meaningful. Only recently, with an approach presented by Liou, Srienc, and Fredrickson (1997), analytical solutions of mass structured, and age-mass structured cell population balances have been obtained. The disadvantage of this approach is that it can be applied only for some simple single-cell growth-rate expressions. Evidently, the generality of cell population balance models is best maintained by the numerical approach. Reports on the numerical solution of such mass-structured models are sparse. We first note the work of Subramanian and Ramkrishna (1971), who employed a combination of the weighted residual method and the successive approximation method. However, this approach has significant memory and computational time demands, and it is limited to the case of linear growth rate where the cell population balance and the substrate concentration equation can be decoupled. In general, this is not possible. Moreover, Sulsky (1994) has addressed a specific nonlinear mass-structured population balance model, with the use of classical finite difference schemes as well as the method of characteristics. For the particular formulation, she has suggested a very efficient adaptive technique to be used in conjunction with the method of characteristics and established via numerical simulations the superiority of the proposed algorithm over classical finite difference methods. However, the model under consideration does not include changes in the environmental conditions, which, when incorpo-

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

rated into the mathematical formulation, can dramatically alter the dynamics of the cell population as well as the behavior of numerical schemes. We also note the work of Godin, Cooper, and Rey (1999a,b), who employed a first-order accurate finite element technique with linear basis functions for the solution of the problem under conditions of changing substrate concentration. The common characteristic of the above approaches is that they are limited to specific mathematical formulations where the intrinsic physiological state functions are of a particular functional form. In order to overcome these limitations, we have presented a one-step, time-explicit, finite difference method, which can treat the single-variable cell population balance model under conditions of changing substrate concentration, for any set of intrinsic physiological state functions (Mantzaris, Liou, Daoutidis, & Srienc, 1999). However, despite its generality, the proposed algorithm has certain limitations that will be discussed in the following. The aforementioned reports refer to the single-variable formulation, which cannot incorporate any internal chemical structure of the single cell and, hence, has obvious modeling and predictive limitations. It is clear that such limitations can be surpassed with the use of the multi-variable formulation, which can model entire metabolic pathways of interest or account for intracellular compartments. Due to the obvious mathematical complexity of such formulations, which consist of multi-dimensional first-order partial integro-differential equations, and the lack of models to describe the intrinsic physiological state functions, there have been no attempts reported in the literature dealing with their numerical solution. There are two fundamental goals of this work: first, to develop different numerical schemes, which belong to different general classes of numerical methods, such as finite difference, spectral and finite element methods, for the numerical solution of multi-variable cell population balance models; second, to reveal the relative advantages and disadvantages of each of the numerical algorithms as well as those of the corresponding classes of numerical methods with the ultimate objective to obtain the most accurate approximations in the smallest possible computational (CPU) times. We consider two levels of comparison. On the first level, we investigate numerical schemes that belong to the same class of numerical methods, and then, we compare the relative effectiveness of each class of numerical methods for the solution of the problem under consideration. The algorithms are compared to each other in terms of numerical stability, accuracy and computational time. Despite the fact that an analytical solution of multivariable cell population balance models is not available, there are specific cases where some information about the solution and/or its nature is known. We consider

1413

three such test problems and use them in order to test the validity and performance of the proposed algorithms. The focus of the present paper is the development and comparative evaluation of finite difference algorithms, while in the subsequent papers, we examine the cases of spectral and finite element methods. Finite difference methods have been used extensively for the numerical solution of first order (hyperbolic) partial differential equations (see Strikwerda, 1989; LeVeque, 1992). This can be mainly attributed to the fact that compared to spectral methods and especially finite element methods are significantly easier to implement. However, they are characterized by several disadvantages: (a) they are problematic in handling complicated geometries, (b) they are sensitive to nonlinearities and (c) in order to yield accurate approximations of first derivatives, they frequently require the ‘invention’ of numerical boundary conditions (i.e. conditions not imposed by the physics of the problem). The first two relative disadvantages of finite difference methods are less important for the problem under investigation due to the following facts: (a) the topology of the physiological state space can be thought as that of a hypercube, which does not require special treatment from the numerical approximation point of view, and (b) the only source of nonlinearity that appears in the system of equations under consideration originates from the coupling of the substrate mass balance with the cell population balance. However, the requirement of numerical boundary conditions by some finite difference algorithms can significantly affect the convergence properties of the algorithm. The issue associated with numerical boundary conditions can be stated as follows. For a given problem and a given numerical approximation that requires extra conditions for the approximation of derivatives at the boundaries of the domain, find those numerical boundary conditions that are capable of retaining the consistency, stability and accuracy properties of the numerical scheme. Despite the progress made by Gustaffson, Kreiss, and Sundstro¨ m (1972) and Osher (1969a,b), who developed a theory (GKSO theory) that applies to a special class of hyperbolic equations and systems, the problem is still far from solved for more general formulations. As a result, appropriate numerical boundary conditions for a specific numerical scheme are still found in practice by trial and error. Despite these possible obstacles, finite difference methods represent an attractive way of approaching the problem under consideration. There are several properties that any finite difference approximation of a certain problem should have. First, the numerical scheme should be consistent with the continuous problem, in the sense that the difference operator converges to the differential operator as the grid spacings go to zero.

1414

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

Moreover, the scheme should also be numerically stable, i.e. the numerical solution over a finite time interval should remain bounded as the grid is refined. Numerical stability and consistency guarantee convergence to the solution of a well-posed continuous problem (see Zachmanoglou & Thoe, 1986 for a definition of a well-posed problem) according to the Lax– Richtmyer equivalence theorem (Richtmyer & Morton, 1967). However, for all practical purposes, it is as important to know the rate of convergence or order of accuracy of the scheme, which is determined by the finite difference formulae that are used to approximate the time and spatial derivatives, as well as the lower-order terms. For detailed mathematical definitions of consistency, numerical stability and accuracy of finite difference methods for partial differential equations, see Strikwerda (1989) In the following, first, we present the general mathematical formulation of the problem. Second, we develop three test problems that will be used in order to examine the validity and performance of the proposed algorithms. Third, we present the numerical schemes, and finally, we comparatively evaluate their stability, accuracy and computational speed. 2. Mathematical formulation of the problem Consider a spatially homogeneous population of cells growing in either a batch or a continuous stirred tank reactor. Let us first define the nutrient concentration vector S = [S1,S2,…,Ss ] as the one containing the concentrations of the s substances in the abiotic environment (substrates), which cells consume in order to grow. Moreover, let us assume that each cell in the population is comprised of a certain number, r, of biochemical components. The r-dimensional vector x= [x1,x2,…,xr ], which contains the amounts of these components in each cell, is called the physiological state vector of this cell. We define the domain of admissible r states as G=[xn,min, xmax]¦ R + , where xn,min and xmax denote, respectively, the vectors containing the minimum and maximum values for the amounts of the r biomass components of the newborn cells. Finally, let xmin denote the vector containing the minimum values of the amounts of the r biomass components of the dividing cells. Clearly, xn,min Bxmin. For the sake of simplicity and without loss of generality, in the developed numerical algorithms that will be presented in the following, we assumed that the minimum and maximum values of the quantities of all biomass components are: xn,min =xmin =0 and xmax =1. If different values for xn,min, xmin and xmax are specified, the proposed numerical algorithms can be applied with a trivial linear transformation that will bring the physiological state space from the domain [0,1] to the domain [xn,min, xmax].

The state of the entire population is described by a time-dependent function N(x,t) which is such that N(x,t)dx represents the number of cells per unit biovolume that at time t have physiological state representation between x and x+ dx. The total number of cells per unit biovolume (cell density) and concentration of the ith biomass component at time t are, respectively, expressed as the zeroth and first moments of the state distribution function at time t:

& &

xmax

Nt (t)=

N(x,t)dx

(1)

xn,min xmax

Nb,i (t)=

xi N(x,t)dx, i =1,…,r.

(2)

xn,min

The integration in both Eqs. (1) and (2) is performed over the entire r-dimensional space of admissible states, G. The sum from 1 until r of all expressions defined in Eq. (2) yields the total biomass concentration at time t. A cell population balance model includes information about nutient uptake, growth, division and birth at the single-cell level. These processes are mathematically described by a set of functions known as intrinsic physiological state functions which, in general, depend on the physiological state of the cell x and the state of the abiotic environment S. The nutrient consumption process is characterized by the s-dimensional consumption rate vector, q(x,S) = [q1(x,S), q2(x,S),…, qs(x,S)] where qi (x,S), i= 1,…,s denotes the single-cell rate of consumption of the ith substrate. The growth process is represented by the r-dimensional single-cell growth rate vector r(x,S) = [r1(x,S), r2(x,S),…, rr (x,S)] where ri (x,S), i= 1,…,r denotes the rate of increase in the amounts of the ith cellular constituent. Finally, the cell division and birth processes are described by the division rate, Y(x,S) and the partition probability density function p(x,y,S), respectively. The latter must satisfy the normalization condition:

&

xmax

p(x,y,S)dx = 1

(3)

xmin

which guarantees that it is a probability density function. It should also be such that the amount of each one of the r biochemical components is conserved at cell division. In particular, since no daughter cell can have greater amounts of any component than the dividing cell from which it originates, the partitioning function, p, should be zero for all daughter cell states that are greater than the states of the corresponding mother cell, i.e.: p(x,y,S) =0, Öxi \ yi, i= 1,…,r.

(4)

Moreover, the probability of a dividing cell with physiological state vector 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.:

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

p(x,y,S) =p(y−x,y,S).

(5)

If the population grows in a CSTR we denote by D the dilution rate and by Sf the vector with the concentrations of the s substances of the abiotic environment at the feed stream of the reactor. We finally assume that no cell death occurs and that cells grow in one stage. Under the assumptions and the process description presented above, the dynamics of the state distribution function N(x,t) were shown to be described by the general cell population balance equation (Fredrickson et al., 1967; Ramkrishna, 1979): (N(x,t) + 9x [r(x,S)N(x,t)]+G(x,S)N(x,t) +DN(x,t) (t =2

&

xmax

G(y,S)p(x,y,S)N(y,t) dy.

(6)

x

The initial condition for the state distribution function is denoted as follows: N(x,0) = N0(x).

(7)

The first term in Eq. (6) denotes accumulation. The second term accounts for the loss of cells with the physiological state vector representation x due to the fact that they grow into bigger cells. The third term represents loss of cells with physiological state vector x due to division leading to the birth of smaller cells. The fourth term is the dilution term describing the rate by which cells exit the reactor. The first term on the right-hand side is the rate of birth of cells with the physiological state vector x originating from the division of all bigger cells. The integration in this term is performed in all r dimensions of the physiological state space and has a lower limit of x due to the fact that cells of physiological state x cannot be born from cells with amounts of biochemical components less than x. Moreover, the factor of two multiplying the integral birth term accounts for the fact that each division event leads to the production of two daughter cells. For a detailed discussion on the statistical foundation of Eq. (6) and the assumptions made in order to derive Eq. (6), see Fredrickson et al (1967). Besides an initial condition, appropriate boundary conditions for the first-order partial differential Eq. (6) are required. As such, we used the so-called regularity conditions first proposed by Behnken, Horowitz, and Katz (1963) and later discussed in the context of cell theory by Eakman et al. (1966) and Fredrickson et al. (1967). From the biological point of view, these conditions can be thought as defining the boundary of the physiological state space. It is due to this fact that the expression ‘containment conditions’ employed by Fredrickson and Mantzaris (submitted for publication) is more appropriate than the previously used ‘regularity conditions’. If we define the boundary, B, of the physiological state space as the set of points where at least

1415

one of the quantities of the r biochemical biomass components obtains either its maximum or minimum value (xmax or xn,min, respectively), then we can mathematically express the containment conditions as: r(x,S)N(x,t)=0, Öx B.

(8)

For a detailed derivation of Eq. (8) from the cell population balance Eq. (6) (see Appendix A). There is a special form of the cell population balance Eq. (6) that deserves special attention due to the fact that it alters the mathematical form of the general formulation. Let us assume that each one of the r biomass components of a dividing cell is partitioned equally into the two daughter cells. This is mathematically described if we choose the partition probability density function to be expressed as the product of the following Dirac delta functions: r



p(x,y,S) = 5 l xi − i=1



yi . 2

(9)

Substituting Eq. (9) into the general Eq. (6) yields the special form of the cell population balance equation under the assumption of equal partitioning: (N(x,t) + 9x [r(x,S)N(x,t)]+G(x,S)N(x,t)+DN(x,t) (t =2r + 1G(2x,S)N(2x,t).

(10)

Notice that now the partial integro-differential equation has become a partial functional differential equation. In addition to that, the assumption of equal partitioning has a consequence on the minimum and maximum values of the physiological state vector for each newborn cell. Since the minimum and maximum values of the physiological state vector for any dividing cell are xmin and xmax, respectively, the corresponding values for each newborn cell must be xn,min = xmin/2 and xmax/2, respectively. Thus, the functional birth term exists in the cell population balance Eq. (10) only in the lower half of the physiological state space, i.e. for x



n

xmin xmax , . 2 2

Also notice that the birth term is now multiplied by a factor of 2r + 1 instead of 2. This is a direct consequence of the form of the partitioning function (Eq. (9)). It guarantees that the following two facts that are independent of the form of the partitioning function hold: (a) the total rate of cell density loss due to division is equal to one half of the total rate of cell density production due to birth and (b) the total rate of biomass component loss due to division is equal to the total rate of biomass component production due to birth, for all biomass components. This can be easily seen by checking the validity of the following two equations, that should hold for any partitioning mechanism:

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1416

&

xmax

xmin

&

1 G(x,S)N(x,t)dx = 2

&

xmax 2

xmin 2

2r + 1G(2x,S)N(2x,t)dx (11)

xmax

xmin

xi G(x,S)N(x,t)dx

&

= 2r + 1

xmax 2

xmin

xi G(2x,S)N(2x,t)dx Öi= 1,…,r.

(12)

2

The cell population balance equation (Eq. (6) or Eq. (10) if the equal partitioning assumption is made) is coupled with the equations describing the dynamics of the abiotic environment. These are the mass balances of the s substrates that are included in the model: dS =D(Sf − S) − dt

&

xmax

q(x,S)N(x,t)dx

(13)

xn,min

The corresponding initial conditions for the abiotic environment, S, are expressed as: S(0) =So.

(14)

The first term in Eq. (13) describes the inlet minus outlet rate from the reactor, whereas the integral term represents the rate of loss of substrate leading to cell growth. In summary, the problem under consideration consists of the cell population balance partial integro-differential Eq. (6), or the partial functional differential equation (Eq. (10)), coupled with the set of the s ordinary integro-differential equations (Eq. (13)), subject to the initial conditions Eqs. (7) and (14) and to the boundary conditions (Eq. (8)). The unknowns of the problem are obviously the state distribution function N(x,t) and the vector of the substrate concentrations, S. Notice that the coupling between the cell population balance equation and the ordinary integro-differential equations occurs only through the dependence of the intrinsic physiological state functions on the concentrations of the substrates. This coupling is the only source of nonlinearity in the model. If the assumption of constant abiotic environment is made, then the problem consists only of Eq. (6) or Eq. (10) subject to the initial condition (Eq. (7)) and boundary conditions (Eq. (8)), and it is linear. Up to this point, only three fundamental assumptions were made. In particular: (a) the cell population is spatially homogeneous. The heterogeneity is only due to the fact that different cells contain different quantities of the biomass components, (b) no cell death occurs and (c) cells grow in a single stage, or there is no cell-cycle structure in the mathematical model. From the numerical point of view, none of these assumptions is really restrictive. In particular, spatial heterogeneity can be mathematically incorporated in the model with a term very similar to the growth term (Ramkrishna, 1985). Moreover, cell death can be modeled with a linear term analogous to the division term (Diekmann, Heijmans, &

Thieme, 1984; Gyllenberg & Webb, 1990). Finally, the incorporation of cell cycle stages can be easily performed in the way that has very recently been suggested by Fredrickson (1999) or by earlier work by Hatzis, Srienc, and Fredrickson (1995). Therefore, the assumptions made do not significantly alter the mathematical nature of the problem under consideration. Hence, the numerical algorithms that will be presented here can, in principle, be applied to the even more general formulations that include spatial heterogeneity, cell death and multiple cell cycle stages with only minor modifications. However, the model presented here will be clearly inadequate to handle ‘compartmentalization’ of subcellular activities with transport resistance between compartments (Fredrickson et al., 1967). There are no reports in the literature that can successfully account for such compartmental organization at the single-cell level and in the corpuscular framework. Moreover, if one accounts for randomness in the growth process by including a diffusion term in the cell population balance equation, the mathematical formulation of the problem becomes significantly different since the cell population balance equation is transformed from hyperbolic in nature to parabolic. Hence, the numerical techniques should be fundamentally different as well. Finally, if nonlinear terms are included in the formulation, as suggested in several mass structured cell population balance models (e.g. Gyllenberg & Webb, 1990), the proposed algorithms should be altered in order to account for this type of nonlinearities.

3. Test problems The problem under consideration is not, in general, analytically solvable. However, there are specific cases of the generalized formulation presented in the previous section, where some aspects or characteristics of the solution are either known or can be mathematically derived. We used these special cases to comparatively evaluate the performance of each developed numerical algorithm based on its ability to capture these characteristics. We considered three such test problems. The parameter values that were used in the simulations for all three test problems and for number of state variables up to three are summarized in Table 1.

3.1. Test problem 1 For the development of this test problem, we made four assumptions: 1. Batch growth (i.e. D= 0). 2. The abiotic environment is constant. Thus, the substrate concentration mass balances (Eq. (13)) are not included in the model.

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

3. The single-cell growth rates are directly proportional to the quantity of each component of the physiological state vector, and the proportionality constant is the same for each component. The value of this constant was chosen to be equal to ln 2/Td, where Td is the doubling time, and is treated as a model parameter. Based on these assumptions, the growth-rate vector can be expressed by the following linear equation: r(x,S) = r(x) =

ln 2 x. Td

(15)

4. Dividing cells partition their material equally to the two daughter cells. Thus, the model consists of the cell population balance equation (Eq. (10)). Under these assumptions, there are two features of the solution that a successful numerical scheme should reproduce: first, the biomass concentration of each component of the physiological state vector should double every doubling time; second, the number density function, defined as: n(x,t)=

N(x,t) Nt (t)

Table 1 Parameter values used in the simulations Test problem 1 Doubling time: Td = 2 h Means and standard de6iations of gaussian di6ision probability density function: Single variable: v1 = 0.575, |1 = 0.125 Two variables: v1 =v2 =0.575, |1 = |2 = 0.125 Three variables: v1 =v2 =v3 = 0.575, |1 = |2 = |3 = 0.125 Partitioning function (single-6ariable): symmetric beta distribution (q= 40) Means and standard de6iations of Gaussian initial condition: Single variable: v1 = 0.2875, |1 = 0.0675 Two variables: v1 =v2 =0.14375, |1 = |2 = 0.03375 Three variables: v1 =v2 =v3 = 0.0958, |1 = |2 = |3 = 0.0225 Test problem 2 vmax = 1 h−1 Ks =0.2 g l−1 Y = 0.5 g g−1 Division probability density function: same as Test problem 1 Partitioning function (single-6ariable): same as Test problem 1 Initial condition: same as Test problem 1 D =0.5 h−1 Sf =10 g l−1 So = 2 g l−1 Test problem 3 k1 =0.8 l g−1 h−1 k2 =0.5 Y =0.5 g g−1 Division probability density function: same as Test problem 1 Initial condition: same as Test problem 1 D =0.4 h−1 Sf =2.0 g l−1 So = 1 g l−1

1417

should follow a periodic behavior with an oscillation period equal to the doubling time. The first feature is independent of the assumption on the partitioning mechanism and can be easily shown by taking the first moment of the cell population balance Eq. (6) or Eq. (10), applying the boundary condition (Eq. (8)) and the fact that biomass is conserved at cell division, using the growth-rate expression (Eq. (15)) and performing a simple integration in time. While the equal partitioning assumption does not affect the biomass concentration dynamics, it has a striking effect on the dynamics of the number density function of the states. In particular, for the case of the single-variable model, Diekmann et al. (1984), have mathematically shown that under the assumptions stated above, the number density function does not reach a stable steady state as it does for other single-cell growth-rate expressions. Furthermore, based on the work of Greiner and Nagel (1985) on semigroup theory, Gyllenberg and Webb (1990) have shown that the number density function exhibits a perfect periodic behavior with a period of oscillation exactly equal to Td for the growth-rate expression (Eq. (22)) and for the single-variable case. Finally, Liou et al. (1997) and Mantzaris et al. (1999) have verified this result via numerical simulations. A nice biological explanation of this phenomenon can be found in the work of Diekmann et al. (1984). They stated that under the assumptions stated above and for the single-variable case, ‘the growth model behaves like a multiplicating machine that copies the size distribution’, or, in other words, ‘the equal size relation is hereditary and extends over the generations’. Despite the fact that the proofs (analytical and numerical) presented to date for this periodic behavior are all related to the single-variable case, we hypothesized that the same behavior should appear for a higher dimensionality of the physiological state space as well. The basis of our hypothesis is that the biological explanation of the periodicity is really independent of the number of components of the physiological state vector. It must be noted that the periodicity check is a very severe test for the numerical algorithms since we demand to capture a characteristic related to the entire distribution and not just a moment of the distribution, such as the biomass concentration. We applied this test to problems with up to three state variables. For the single-variable case, we also considered the more general case of unequal partitioning, where the number density function does not exhibit a periodic behavior as shown by Mantzaris et al. (1999). In this case, we assumed that the partitioning mechanism is independent of the substrate concentration and can be captured with a symmetric beta distribution with a parameter of q:

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1418

p(x,y)=

   q−1

1 1 x B(q,q) y y

1−

q−1

x y

.

(16)

The features that we are interested in capturing for this test problem with the developed algorithms are independent of the division rate. Despite the fact that we used several different functional forms that have been presented in the literature, all the results that will be shown here correspond to the following division rate model:

:

;

summing over all biomass components, we obtain the following system of equations: dNb vmaxS = ·Nb − DNb dt Ks + S

(20)

dS 1 vmaxS = D(Sf − S)− N dt Y Ks + S b

(21)

where hi (xi ) is the division probability function with respect to the ith component of the physiological state vector. The model described by Eq. (17) has been extensively used in the single-variable case (r= 1) by several researchers (Subramanian & Ramkrishna, 1971; Liou et al., 1997; Mantzaris et al., 1999, etc.). In this work, for the division probability density functions, we considered truncated Gaussian distributions with a mean of vi and a standard deviation of |i. Notice that according to this model, the dependence of the division rate on the substrate concentration is only via the division rate dependence on the single-cell growth rates. Although, in this test problem, the growth rate is independent of the substrate concentration, we used this general model for the other test problems where there is a dependence on the substrate concentration.

where Nb(t) is the total biomass concentration at time t. Eqs. (20) and (21) constitute the famous Monod model (Monod, 1942). The previous analysis shows that under assumptions (3.2.a)–(3.2.d), the cell population balance model and the Monod model should make identical predictions for the biomass and substrate concentrations. This agreement is totally independent of the division rate function and partitioning function. Notice also that the agreement on the total biomass and substrate concentrations is also independent of the number of components included in the physiological state vector. Despite the fact that the analytical solution of the system (Eqs. (20) and (21)) is not obtainable, this system of the two nonlinear ordinary differential equations can be integrated easily and rapidly using standard time-integration methods. We applied this test for our numerical algorithms for up to three state variables. For the single-variable simulations, we considered both cases of equal and unequal partitioning. For higher numbers of state variables, we applied this test only for the equal partitioning case (Eq. (10)). We used the division rate model described by Eq. (17) for all simulations and all dimensions.

3.2. Test problem 2

3.3. Test problem 3

For this test problem, we made the following assumptions: 1. Growth on a continuously stirred tank reactor (CSTR). 2. Growth on a single substrate (s =1) with concentration denoted by S. 3. The single-cell growth rate vector is proportional to the physiological state vector and follows Michaelis–Menten kinetics, i.e.:

Despite the fact that, in the previous test problem, we included a change in the substrate concentration, we still kept a linear law for the single-cell growth rate dependence on the physiological state vector. In order to test the algorithms under more severe nonlinear growth rate expressions, we considered a third test problem. This test problem was designed for the twovariable case only (r=2), where the mass of each cell comprises two components with masses x1 and x2. The following assumptions were made: 1. Growth on a continuously stirred tank reactor (CSTR). 2. Growth on a single substrate (s= 1) with concentration denoted by S. 3. The two single-cell growth rates are expressed by the equations:

r

i=1



& 

hi (xi )

G(x,S) = % 1−

r

xi

hi (y)dy

i=1

xi,min r

= % ki (xi ) % ri (x,S) i=1

r(x,S) =

r

% ri (x,S)

(17)

i=1

vmaxS x Ks +S

(18)

4. The single-cell consumption rate of the single substrate is expressed by the equation: q(x,S) =

1 vmaxS r %x Y Ks +Si = 1 i

(19)

where Y is a constant biomass yield. Substituting expressions Eqs. (18) and (19) into Eqs. (A4) and (13), respectively, applying the definition of the biomass component concentrations (Eq. (2)) and

r1(x1,x2,S)= k1(x1 + x2)S−k2 r2(x1,x2,S)= k2

x1x2 . x1 + x2

x1x2 x1 + x2

(22) (23)

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

4. The single-cell consumption rate of the single substrate is expressed by the equation: q(x1,x2,S)=

1 k S(x1 +x2) Y 1

(24)

Substituting the sum of Eqs. (22) and (23) into Eq. (A4), Eq. (24) into Eq. (13) and using the definition of the biomass component concentrations (Eq. (2)), we obtain the following equations describing the dynamics for the biomass component, total biomass and substrate concentrations:

& &

dNb,1 = k1SNb − k2 dt

x1,max

x1,min

x2,max

x2,min

x1x2 N(x1,x2,t)dx1dx2 x1 +x2 (25)

−DNb,1

& &

dNb,2 = k2 dt

x1,max

x1,min

x2,max

x2,min

x1x2 N(x1,x2,t)dx1dx2 −DNb,2 x1 +x2 (26)

dNb =k1SNb −DNb dt

(27)

dS 1 = D(Sf − S)− k1SNb. dt Y

(28)

Notice that the dynamics of the total biomass and substrate concentrations are independent of the two biomass component concentrations. Also, the application of the proposed numerical algorithms to the full model (functional partial differential equation and integro-ordinary differential equation) should agree with the predictions of the simple nonlinear ODE system of Eqs. (27) and (28) on the biomass and substrate concentrations. Eqs. (27) and (28) are the same equations as those of the well-known continuum model proposed by Williams (Williams, 1967; Bailey & Ollis, 1986). However, the equations describing the biomass component concentrations according to the continuum formulation are different when compared to those resulting from the corpuscular framework. In particular, according to the original Williams model, the lumped biomass component concentration Nb,1,Nb,2 dynamics are expressed by the equations: dN( b,1 N( b,1N( b,2 = k1SNb − k2 −DN( b,1 dt N( b,1 +N( b,2

(29)

dN( b,2 N( b,1N( b,2 = k2 −DN( b,2. dt N( b,1 + N( b,2

(30)

The difference between the system represented by Eqs. (29) and (30) and the corresponding system represented by Eqs. (25) and (26) is a direct consequence of the assumption of homogeneous growth that is inherent in all continuum models. Notice that the different terms between the two sets of equations are nonlinear with respect to the states. However, the first terms of the right-hand side of Eqs. (29) and (25) are the same and

1419

are linear with respect to the states. Thus, differences in the predictions of corpuscular and continuum models on average population properties are expected to arise in the presence of nonlinearities (with respect to the states) in the kinetic expressions. For the simulations, we used Eq. (17) for the division rate and assumed equal partitioning (Eq. (10)).

4. Development of finite difference algorithms We define a multi-dimensional equally spaced grid consisting of M+ 1 grid points in each dimension of the physiological state space. The distance between two neighboring grid points in each dimension is h= 1/M. In general, h can be different for each dimension of the physiological state space. The first grid point (i=1) always corresponds to the zero value, whereas the last grid point (i = M+ 1) corresponds to the value of 1. Thus, for each coordinate of the physiological state space: xi = (i −1)h. We also use dt to denote the time step. The solution of the cell-population balance equation at the nth time step tn = n dt is approximated by the numerical solution at the prespecified grid points. Thus, for the single-variable case and using the notation ‘  ’ to denote the numerical approximation, we have: N(xi,tn ): N0 (xi,tn ) N0 ni . Also, for each of the concentrations of the substrates, we have: S(tn ): S0 (tn ) S0 n. For simplicity in the notation, in the following, we will present the developed numerical algorithms for the case where the dimension of the physiological state space is 1 (r=1) and the abiotic environment consists of a single substrate (s= 1).

4.1. A hybrid scheme A good starting point for the algorithmic development of finite difference schemes is a numerical scheme that was developed and presented in our previous work (Mantzaris et al., 1999). This scheme was applied to the single-variable problem and consisted of a combination between the Lax–Friedrichs and leapfrog schemes. The scheme is one-step and time-explicit. We will use the results of the application of this scheme for comparison purposes with the rest of the methods that were developed in this work.

4.2. The Lax–Wendroff scheme In order to obtain more accurate approximations, we also considered applying the idea that Lax and Wendroff had for the numerical solution of systems of equations satisfying conservation laws (Lax & Wendroff, 1960). They first expressed the solution of the equations at the (n+1)th time step (t= tn + 1 =(n+ 1) dt) in a second-order Taylor series expansion with

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1420

respect to time around the solution at the nth time step (t = tn = n dt). By applying this to the solution of our system of the partial and ordinary differential equations, and neglecting the terms with order higher than three O(dt 3), we obtain the approximate solution at the (n+ 1)th time step and at the ith grid point: N0 ni + 1 = N0 ni + dtN0 nt,i + S0 n + 1 = S0 n + dtS0 nt +

dt 2 n N0 tt,i 2

dt 2 n S0 2 tt

(31) (32)

where N0 nt,i, N0 ntt,i denote the approximations of the first and second time derivatives of the number of cells, respectively, at the nth time step and at the ith grid point, and S0 nt , S0 ntt are the approximations of the substrate concentration first and second time derivatives, respectively, at the nth time step. Using the general form of the cell population balance Eq. (6) for a single variable and a single substrate, the substrate mass balance (Eq. (13)) and the model (Eq. (17)) for the division rate, we obtain: Nt (xi,tn )= − [r(xi,S(tn ))N(xi,tn )]x

&

− k(xi )r(xi,S(tn ))N(xi,tn ) −DN(xi,tn ) 1

k(y)r(y,S(tn ))p(xi,y)N(y,tn )dy

+2

xi

St (tn )= D(Sf − S(tn )) −

&

(33)

1

q(x,S(tn ))N(x,tn )dx.

(34)

0

Differentiating Eqs. (33) and (34) w.r.t. time and applying the chain rule, we obtain the following equations for the second time derivatives: Ntt (xi,tn )= − St (tn )[rs (xi,S(tn ))N(xi,tn )]x −[r(xi,S(tn ))Nt (xi,tn )]x − St (tn )k(xi )rs (xi,S(tn ))N(xi,tn ) −k(xi )r(xi,S(tn ))Nt (xi,tn ) −DNt (xi,tn ) + 2St (tn )

&

&

1

k(y)rs (y,S(tn ))p(xi,y)N(y,tn )dy

xi 1

+2

xi

k(y)r(y,S(tn ))p(xi,y)Nt (y,tn )dy

&

1

Stt (tn )= − DSt (tn )− St (tn ) −

&

(35)

qs (x,S(tn ))N(x,tn )dx

0

4.3. Leapfrog schemes

1

q(x,S(tn ))Nt (x,tn )dx

For the approximation of the integral terms, we applied the Trapezoid rule, while for the approximation of all first order spatial derivatives, we used the classical second-order centered difference formulae. We directly approximated the derivative of the product of the two functions instead of performing the spatial differentiation first, as it is common in approximations of spatial derivatives for conservation law equations (LeVeque, 1992). The approximation of the first spatial derivatives cannot be performed for the first (i= 1) and last (i= M+1) grid points. This is due to the fact that at i= 1, the (i− 1)th grid point does not belong to the domain, whereas at i =M+ 1, the (i+ 1)th grid point does not belong to the domain. In the case where the value of the growth rate function is different than zero at the boundary points, we can directly apply the containment conditions (Eq. (8)) and set the value of the state distribution function at these boundary points equal to zero. However, if the growth-rate function is zero at any of these points, then the containment conditions cannot be used to find the corresponding values for the state distribution function. There are several ways to deal with this difficulty. One way is to change the formulae for the first derivative approximations at these two grid points by employing a forward difference formula for the first grid point and a backward difference for the last grid point. Since the centered difference formulae used are second order and in order to maintain the overall accuracy of the scheme, both formulae for the first and last grid points should be second order as well. An alternative way is to use the same formulae for the boundary points, but apply the containment conditions (Eq. (8)) for the non-existent points, i.e. for i= 1, where the non-existent point is i− 1= 0, set: r(x0,S0 n)N0 n0 =0, while for i=M+1, where the non-existing point is i+ 1= M+ 2, set: r(xM + 2,S0 n)N0 nM + 2 = 0. The idea is the same for the approximation of the two other first spatial derivatives that appear in the derivation of the Lax– Wendroff scheme. We applied both ideas in the developed algorithms and saw that the application of the containment conditions at the non-existent points as numerical boundary conditions yields significantly more stable approximations for every test problem and for any number of variables.

(36)

0

where the subscript s in rs (x,S(tn )) and qs (x,S(tn )) denotes differentiation of the growth and nutrient uptake rates, respectively, w.r.t. the substrate concentration. Thus, the Lax – Wendroff idea applied to this problem consists of approximating the first (Eqs. (33) and (34)) and second time derivatives (Eqs. (35) and (36)) of the population balance equation and substrate mass balance.

The Lax –Wendroff scheme is a representative example of a so-called dissipative scheme (see Strikwerda, 1989 for a definition). This type of schemes can be very useful in decreasing the frequency of high-frequency oscillations that might appear in the solution of certain problems. If such high-frequency oscillations exist, they can significantly reduce the accuracy of an approximation. However, the mathematical nature of hyperbolic partial

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

differential equations is non-dissipative. Thus, when approximating the solution of such equations, we ideally want a numerical scheme that can maintain the non-dissipative nature of the continuous problem. Due to the fact that the appearance or not of high-frequency oscillations in the numerical solution of the problem under consideration depends on the specifics of the formulation and cannot be known a priori, it is worth examining the behavior of non-dissipative as well as dissipative schemes. Perhaps the best-known time-explicit non-dissipative schemes that are used for the solution of hyperbolic partial differential equations are the so-called leapfrog schemes. These are two-step algorithms where the approximation at each time step requires the use of the approximate solution at the two previous time steps. The time derivative at the nth time step is approximated by a second-order accurate, with respect to time central difference formula, which can be applied only for n\ 1. Thus, for initiation purposes, we need a one-step numerical scheme. As such, we used the Lax– Wendroff scheme presented earlier. According to leapfrog schemes, the first spatial derivative in the partial differential equation (Eq. (33) in this case) is approximated by a central difference formula. The order of this approximation (2nd, or 4th, or 6th, etc.) is one characteristic that differentiates leapfrog schemes from each other. It is obvious that an increase in the order of accuracy for the approximation of the spatial derivatives increases the overall accuracy of the scheme. Thus, we can achieve the same overall accuracy using coarser grids in space for the same time step. As a result, given a certain value for the time step, the same accuracy can be achieved in less computational time. However, it was found that the increase in the order of accuracy for the approximation of the first spatial derivatives significantly decreases the stability limits of the algorithm. Thus, for a given number of grid points, the maximum time step that can be used in conjunction with higher-order approximations of the spatial first derivatives is smaller than that of lower order approximations, thus resulting in an increase in the computational time requirements. Therefore, with lower-order accurate approximations of the first spatial derivatives and in order to achieve certain desirable levels of accuracy, we might need more grid points, which results in increased computational time, but we are allowed to use greater time steps, leading to decreased computational time. As is apparent from the above observations, there must be an optimum order of approximation for the spatial derivatives above which the decrease in the stability limit becomes more important than the increase in the accuracy of the method. We have found that the best result is obtained for a 6th-order accurate approximation of the first spatial derivatives. For higher-order approximations, the increase in computational time caused by the decrease in the stability limit is not compensated by the higher-accuracy w.r.t. space.

1421

The issue of deciding on the type of first derivative approximation close to the boundaries of the physiological state space domain is present here as in the case of the Lax–Wendroff algorithm. Since the optimal approximation of the first spatial derivative is 6th order, numerical boundary conditions are required for the spatial derivative approximation at the 2nd, 3rd, (M− 1)th and Mth grid points. Contrary to the Lax–Wendroff scheme, it was found that the more stable approximations are obtained if we reduce the order of spatial approximation close to the boundaries of the physiological state space, instead of using the same formulae and set the non-existent terms equal to zero (i.e. apply the containment conditions at the non-existing points as numerical boundary conditions). This choice obviously reduces the accuracy of the spatial approximation close to the boundaries but increases the stability limit of the algorithm. It must be noted that the choice of numerical boundary conditions influences the performance of leapfrog schemes more than any other type of method developed in this work. We have also considered applying the containment conditions as numerical boundary conditions in conjunction with the formula for the approximation of the first spatial derivative, as well as quasi-characteristic extrapolation of numerical boundary conditions suggested by Strikwerda (1989). Another characteristic that differentiates leapfrog schemes from each other is the way that the approximation of the lower-order terms is performed. Several different ways have been suggested in the literature. In order to illustrate them, let us consider the approximation of the dilution term DN(xi,tn )at grid point xi and at the nth time step that appears in Eq. (33). One way is to just use the numerical solution at the nth time step, i.e. DN(xi,tn ): DN0 ni . The second way is to use the average of the lower-order term approximations at the nth and (n− 1)th time steps, i.e. DN(xi,tn ):

DN0 ni + DN0 ni − 1 . 2

Finally, it is also possible to average between the approximations at the nth and (n+1)th time steps, i.e. DN(xi,tn ):

DN0 ni + DN0 ni + 1 . 2

The drawback of the latter approach is that the scheme becomes time-implicit, which leads to increased computational demand, whereas the other two approaches maintain the explicit nature of the leapfrog scheme. We have applied both of the other two approaches for the approximation of the lower-order terms, and it was seen that the averaging between the current (nth) and previous ((n − 1)th) time steps significantly increases the numerical stability of the scheme. Moreover, notice that, as opposed to the Lax–Wendroff scheme, which

1422

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

requires two evaluations of the integral birth term at each time step, in the case of leapfrog schemes, only one evaluation is required if, at each time step, we store the integral birth term.

4.4. High-order time-explicit schemes As will be illustrated in Section 5, the leapfrog scheme with the 6th-order accurate approximation of the first spatial derivative proved to be more successful than both the hybrid and Lax– Wendroff algorithms for every test problem. The major drawback of leapfrog schemes, however, was the significant decrease in the stability limits caused by the use of the more accurate, 6th-order approximation of the first spatial derivative. We speculated that this might be due to the fact that leapfrog schemes are, by construction, only second-order accurate in time. Thus, we also considered the possibility of deriving finite difference schemes that are high-order accurate w.r.t. both space and time. We studied two different types of time integration techniques in order to achieve high-accuracy w.r.t. time: one-step and multi-step time-explicit methods. One-step time-explicit methods only use the solution estimated at the previous time step in order to approximate the solution at the current time step, while in time-explicit multi-step methods, the approximations at more than one previous time step are required. It is obvious that multi-step methods have greater storage requirements than one-step methods.

4.4.1. One-step Runge– Kutta schemes Canuto, Hussaini, Quarteroni, and Zang (1988) have illustrated that as the order of accuracy of one-step Runge– Kutta (RK) methods increases, the stability limit of the approximation increases. Despite the fact that their work considered spectral approximations of the much simpler advection equation with periodic boundary conditions, we have found that their conclusion is also true for finite difference approximations of the generalized cell population balance problem considered here. Increased stability limits allow the use of greater time steps for a given spatial discretiztion (fixed number of grid points in space) which in turn reduces the CPU time. However, the drawback of one-step methods is that as the order of accuracy of the time derivative approximation increases, the number of right-hand side evaluations at each time step increases as well, thus leading to an increase in CPU time. It is clear that there must be an order of accuracy for the RK approximation, where the decrease in CPU time that results from the increased stability limit compensates, in an optimal way, the increase in CPU time caused by the increased number of right-hand side evaluations. We have constructed 2nd, 3rd, 4th and 5th order RK schemes by applying the well established

methodology of Taylor series expansions used for deriving such schemes. For the problem under consideration, it was found that 4th order RK methods lead to the optimal result. On the other hand, as is the case for leapfrog schemes, for any given RK approximation scheme in time, the stability limits of the overall algorithm decrease as the order of approximation of the first spatial derivatives increases. Thus, high order approximations of the first spatial derivatives improve the overall accuracy but shrink the stability region of the scheme. Therefore, again there must exist an optimal compromise between the increased accuracy that can lead in decreased CPU times with the use of coarser grids, and the decreased stability region that can lead in increased CPU time requirements. By employing centered difference approximations of different order for the first spatial derivative, we have found that for the test problems considered here where the solutions are smooth enough, the optimal result is obtained for 10th order accurate, centered difference approximations of the first spatial derivatives. It is obvious that this approximation requires the application of appropriate numerical boundary conditions for grid points close to the two boundaries of the domain (i= 2, 3, 4, 5, M-3, M-2, M-1, M). For this type of algorithm the application of the containment conditions at all non-existent grid points was found to yield the best result in terms of numerical stability. For the first and last grid points (i= 1, M+ 1) we apply the containment conditions in Eq. (8), thus setting the corresponding values of the state distribution function equal to zero, if the values of the growth rate are non zero. In the case where the growth rate is indeed zero at the boundary points, the state distribution function was estimated using the scheme with the containment conditions applied as numerical boundary conditions wherever they were required. At each of the right-hand side evaluations that are performed at each step of the time integration, the lower order terms were evaluated by simply substituting the numerical approximation at the corresponding grid point and time step. Since the Runge–Kutta algorithm is 4th order accurate, there exist four right-hand side evaluations at each time step and hence, four evaluations of the integral birth term as well. This is expected to significantly increase the CPU time requirements of the scheme in the case of unequal partitioning.

4.4.2. Multi-step schemes Multi-step time-explicit schemes can achieve high accuracy in time with fewer function evaluations per time step than one-step schemes. Thus, it is possible that the use of such schemes for the time integration can produce a significant reduction in the CPU time required for obtaining very accurate solutions to the

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

problem, when compared to one-step Runge– Kutta algorithms. Contrary to the case of one-step methods, we have found that the stability limit of multi-step algorithms decreases as the order of accuracy of the time derivative approximation increases. This result is in agreement with a similar result that Canuto et al. (1988) obtained from studying time integration techniques for spectral approximations to the advection equation. Therefore, the increase in the order of accuracy of the time derivative approximation is unavoidably accompanied by a reduction in the stability limit, which could result in increased CPU time requirements. The advantage however, is that this increase in the order of accuracy does not require more function evaluations per time step as it is the case with one-step methods. The derivation of all multi-step algorithms is based on numerical quadrature rules. Depending on the number of previous time steps that are used in the evaluation of the solution at each time integration step we distinguish between single-step (one previous step) and multi-step (more than one previous steps) multi-step algorithms. Furthermore, depending on the order of the polynomial that is used to approximate the function f in the integration interval, we obtain different orders of accuracy for the time integration technique. Finally, there exist time explicit and time-implicit multi-step schemes. For a detailed discussion on multi-step algorithms, see Young and Gregory (1988). We have focused on explicit schemes, which may in general have lower stability limits when applied to partial differential equations, but are significantly less computationally expensive than implicit schemes. We have found that for the problem under consideration, single-step algorithms of a certain order yield more stable approximations than multi-step algorithms of the same order of accuracy. Perhaps, the best-known family of time-explicit, single-step, multi-step algorithms is the Adams– Bashforth (AB) family (Gerald & Wheatley, 1989). The corresponding time-implicit, single-step, multi-step algorithms are known as Adams– Moulton (AM) methods. All AB schemes require only one function evaluation per time step, but generally have very low stability limits. We have considered up to 5th order of accuracy AB schemes. An attractive alternative to AB schemes is offered by the Adams– Predictor – Corrector (APC) methods, which were found to yield increased stability limits at the expense of one extra function evaluation per time step. In these schemes, at each time step, an AB step of order n −1 is used as the predictor, in conjunction with an AM step of order n as the corrector. The overall order of accuracy of such schemes is n. It is of course possible to use an AB predictor step of the same order of accuracy as the AM corrector step. However, it was found that using an

1423

n− 1 order accurate AB predictor step yields more stable approximations. By comparing AB and APC algorithms of the same order, we have seen that the increased stability limit of APC schemes offers a more significant reduction in CPU time than the increase in CPU time caused by the extra function evaluation per time step. Moreover, the best compromise between accuracy w.r.t. time and stability was offered by the 4th-order APC method, although 3rd- and 5th-order schemes were tested as well. Their derivation was performed using the method of undetermined coefficients. The 4th-order APC method that was found to be optimal consists of an AB3 predictor step and an AM4 corrector step (see Young & Gregory, 1988). This APC scheme can be applied only for n\ 2. Thus, for initiating the APC algorithm, another one-step algorithm is required. In order to retain the 4th-order accuracy w.r.t. time that the APC scheme has, we used the 4th order RK algorithm for the first two time steps. As is the case with leapfrog and Runge–Kutta algorithms, for an APC algorithm of a certain order, the stability limit is reduced when more accurate approximations of the first spatial derivatives are used. The optimal result w.r.t. CPU time requirements for very accurate approximations was obtained for 8th-order accurate centered difference approximations of the first spatial derivative. Close to the boundaries of the computational domain where the 8th order approximation of the first spatial derivative cannot be applied (i=2, 3, 4, M− 2, M− 1, M), lower-order accurate formulae were used, as opposed to RK schemes where the best choice was found to be the application of the containment conditions as numerical boundary conditions.

4.5. Time-implicit schemes The most important disadvantage of time-explicit schemes is their low stability limit, which may impose severe restrictions on the maximum time step that can be used for a given number of grid points in space. Moreover, it was seen that for all the developed timeexplicit algorithms, there is an unfavorable relationship between the order of accuracy of the spatial derivative approximation and the stability limit: the higher the order of accuracy of the spatial derivative approximation, the lower the stability limit. Time-implicit schemes are generally more stable. However, their implementation requires the solution of a system of equations at each time step as opposed to time-explicit schemes where the unknown at each grid point in space at each time step is decoupled from the neighboring unknowns at the same time step. The system of equations that needs to be solved is linear (nonlinear) if the continuous problem is linear (nonlinear). Hence, in most cases, time-implicit algorithms are computationally more demanding than time-explicit algorithms.

1424

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

In order to compare the efficiency of time-implicit methods versus that of the developed time-explicit methods for the problem under consideration, we implemented the Crank– Nicholson, 2nd-order accurate time-implicit algorithm for the time integration in conjunction with the Newton– Raphson technique for the solution of the nonlinear (in general) system of equations at each time step. The Crank– Nicholson scheme can also be thought of as an Adams– Moulton 2nd-order time-implicit multi-step method. The great advantage of this time integration technique is that it is unconditionally stable (Strikwerda, 1989). Therefore, one can use any consistent formula for the approximation of the spatial derivative and the lower-order terms in conjunction with this time-integration method. In order for the algorithm to be comparable with the most accurate time-explicit scheme developed, and since the order of the spatial derivative approximation is completely dissociated from stability, we have considered the 10th-order approximation of the spatial derivatives. Similar to the Runge– Kutta 4– 10 scheme, the containment conditions were applied at the non-existent grid points, as numerical boundary conditions.

5. Results In this section, we will comparatively evaluate the performance of the numerical algorithms that were presented above. For simplicity, we will adopt the following abbreviations: CN2 –10: 2nd-order Crank– Nicholson for the time integration and 10th-order accurate approximation of the first spatial derivatives. RK4 –10: 4th-order Runge– Kutta algorithm for the time integration and 10th-order accurate approximation of the first spatial derivatives. APC4 –8: 4th-order Adams– Predictor – Corrector for the time integration and 8th-order accurate approximation of the first spatial derivatives. LF2 –6: 2nd-order Leapfrog algorithm for the time integration and 6th-order accurate approximation of the first spatial derivatives. LW2 –2: The Lax – Wendroff algorithm, which is second-order accurate in both space and time.

5.1. The single-6ariable case We applied the developed algorithms for test problems 1 and 2 for both cases of equal and unequal partitioning and for the parameter values presented in Table 1. The final time for all simulations was Tf = 20 h. For the first test problem, this final time corresponds to 10 doubling times, whereas for the second problem, simulations with the Monod ODE model have shown

that after 20 h, the system is at steady state. In order to examine whether the developed algorithms are indeed convergent, we performed simulations using very fine grid spacings (dt =5×10 − 4, dx = 2× 10 − 3, u=0.25). Fig. 1 shows the simulation results for all algorithms for these grid spacings and for the first test problem for both cases of equal and unequal partitioning. First, notice that all numerical algorithms produce results that are virtually indistinguishable. For simulations performed with each of the algorithms using finer grid spacings, no appreciable change in the numerical approximation was observed. These observations indicate that all numerical algorithms indeed converge to the solution of the continuous problem. Second, notice that the algorithms correctly predict that the biomass concentration doubles every doubling time (2 h). This is independent on the assumption of equal or unequal partitioning. Thus, all algorithms satisfy the conservation of biomass at cell division as they are supposed to. Third, in the case of equal partitioning, all algorithms capture the perfect periodic behavior of the middle point of the number density function, which is vanished in the case of unequal partitioning where the middle point of the number density function reaches a steady state. This fundamental difference between the cases of equal and unequal partitioning is better shown in Fig. 2, which shows the predicted dynamics of the entire number density function as well as those of the state distribution function. Finally, all algorithms agree on their predictions for the state distribution function at t= Tf (Fig. 1c and f) Fig. 3 shows the same characteristics of the solution for the second test problem and for the same grid spacings. First, as in the case of the first test problem, notice that all numerical algorithms converge to the same solution. Second, the dynamics of the biomass and substrate concentrations predicted by the numerical algorithms are the same as those predicted by the Monod model described by the ordinary differential equations Eqs. (20) and (21). Moreover, in the case of equal partitioning, the middle point of the number density function exhibits a sustained oscillatory behavior, which is not present in the case of unequal partitioning. Thus, despite the fact that the biomass concentration (first moment of the state distribution function) reaches a steady state, the number density function does not do so in the case of equal partitioning. Finally, notice that in the case of equal partitioning, the state distribution function is bimodal as opposed to the case of unequal partitioning. Figs. 1 and 3 show that all the developed algorithms converge to the same solution. Since vastly different schemes produce identical results, this converged solution must be a very accurate approximation of the solution to the continuous problem. However, the results shown in Figs. 1 and 3 do not reveal the relative

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1425

Fig. 1. Single-variable case. Illustration of convergence. Test problem 1. Simulation results for dt=5 ×10 − 4 and dx= 2 ×10 − 3. Plots (a), (b) and (c): equal partitioning. Plots (d), (e) and (f): unequal partitioning. Solid line, open circles: Hybrid scheme. Dashed line, filled circles: the Lax–Wendroff scheme. Solid line, open squares: the Leapfrog 2 – 6 scheme. Dashed lines, filled squares: the Adams – Predictor – Corrector 4 –8 scheme. Solid line, open triangles: the Runge –Kutta 4 – 10 scheme. Dashed line, filled triangles: the Crank – Nicholson 2 – 10 scheme.

1426

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

Fig. 2. Single-variable case. Test problem 1. The converged numerical approximation (dt= 5 × 10 − 4 and dx= 2 ×10 − 3) describing the dynamics of the number density function n(x,t) and state distribution function N(x,t). Top plots (a and b): equal partitioning. Bottom plots (c and d): unequal partitioning.

advantages and disadvantages of the developed algorithms. In order to compare the accuracy of the schemes for the first test problem, we performed simulations with coarser grids and with a ratio u= dt/dx close to the stability limit of the least stable numerical scheme. Thus, for the first test problem, we chose dt = 3.571×10 − 2 and dx =2.94 ×10 − 2 (u =1.214) for the case of equal partitioning, whereas for the case of unequal partitioning, the corresponding choice was dt = 3.846×10 − 2 and dx = 5 ×10 − 2 (u = 1.3). For the second test problem, using the same reasoning, we chose dt = 9.434×10 − 3, dx = 2.27 ×10 − 2 (u =0.415) for equal partitioning and dt =1.515 × 10 − 2, dx= 4.54× 10 − 2 (u= 0.33) for unequal partitioning. The

results of the comparison are shown in Figs. 4 and 5, for the first and second test problems, respectively. Also, we show the converged solution presented in Figs. 1 and 3, respectively, in order to indicate the ‘distance’ from the converged result. First, notice that for both test problems and for both cases of equal and unequal partitioning, all numerical algorithms predict the dynamics of the biomass and substrate concentrations very accurately, despite the fact that the chosen grids were much coarser than those used in the previous simulations. Small deviations are observed only for the hybrid scheme predictions of the biomass concentration for the first test problem and for the case of unequal partitioning as well as for the

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1427

Fig. 3. Single-variable case. Illustration of convergence. Test problem 2. Simulation results for dt=5 ×10 − 4 and dx= 2 ×10 − 3. Plots (a), (b) and (c): equal partitioning. Plots (d), (e) and (f): unequal partitioning. Solid line: the results from the time integration of the 2 × 2 Monod model (Eqs. (20) and (21)). Solid line, open circles: hybrid scheme. Dashed line, filled circles: the Lax – Wendroff scheme. Solid line, open squares: the Leapfrog 2–6 scheme. Dashed lines, filled squares: the Adams –Predictor – Corrector 4 – 8 scheme. Solid line, open triangles: the Runge – Kutta 4 – 10 scheme. Dashed line, filled triangles: the Crank –Nicholson 2 –10 scheme.

1428

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

Fig. 4. Single-variable case. Comparison between the developed numerical schemes. Test problem 1. Plots (a), (b) and (c): equal partitioning: dt= 3.571 × 10 − 2 and dx= 2.94× 10 − 2. Plots (d), (e) and (f): unequal partitioning: dt= 3.846 × 10 − 2 and dx= 5 ×10 − 2. Solid line: converged results. Solid line, open circles: Hybrid scheme. Dashed line, filled circles: the Lax – Wendroff scheme. Solid line, open squares: the Leapfrog 2 –6 scheme. Dashed lines, filled squares: the Adams –Predictor –Corrector 4 – 8 scheme. Solid line, open triangles: the Runge – Kutta 4 – 10 scheme. Dashed line, filled triangles: the Crank –Nicholson 2 –10 scheme.

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1429

Fig. 5. Single-variable case. Comparison between the developed numerical schemes. Test problem 2. Plots (a), (b) and (c): equal partitioning: dt= 9.434 ×10 − 3 and dx= 2.27×10 − 2. Plots (d), (e) and (f): unequal partitioning: dt= 1.515 × 10 − 2 and dx= 4.54 ×10 − 2. Solid line: converged results. Solid line, open circles: Hybrid scheme. Dashed line, filled circles: the Lax –Wendroff scheme. Solid line, open squares: the Leapfrog 2 – 6 scheme. Dashed lines, filled squares: the Adams – Predictor– Corrector 4 – 8 scheme. Solid line, open triangles: the Runge –Kutta 4–10 scheme. Dashed line, filled triangles: the Crank – Nicholson 2 – 10 scheme.

1430

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

predictions of biomass and substrate concentrations for the second test problem and for the case of unequal partitioning. However, this agreement is not observed in the predictions of the middle point of the number density function, where the algorithms exhibit significant differences for both test problems and for the case of equal partitioning. The RK4– 10, APC4 – 8 and CN2 –10 schemes are clearly the most accurate, since they produce a periodic approximation for the first test problem, virtually indistinguishable from the converged result. The LF2 –6 approximation overpredicts the amplitude of the periodic solution as time increases. This is a direct consequence of the algorithm’s strictly non-dissipative nature. On the contrary, notice that the LW2– 2 algorithm, which is dissipative (of order 4 for the one-way, constant coefficient, wave equation as Strikwerda (1989) shows), tends to smooth the periodic nature of the solution to the continuous problem. As a result, as time increases, the LW2– 2 fails to capture the periodic behavior of the solution. This ‘smoothing’ effect is even more pronounced in the case of the Hybrid algorithm, which is dissipative as well. The performance of the numerical algorithms is analogous for the second test problem and for the equal partitioning case, where the converged solution exhibits a nonperiodic, oscillatory behavior. For both test problems and for the case of unequal partitioning, the less accurate algorithms exhibit transient deviations from the converged result; however, they manage to predict the steady-state result accurately. Thus, as expected based on the construction of the numerical algorithms, the schemes can be ordered in decreasing accuracy as follows: RK4 – 10\ APC4 –8 \ CN2 – 10 \ LF2 – 6\ LW2 – 2\Hybrid. Notice the oscillations produced by the LF2 –6 algorithm for the state distribution function at t =Tf close to the right boundary of the domain in both test problems. These are due to the fact that, for illustration purposes, the simulations were performed slightly above the stability limit of the algorithm. As is often the case, the most accurate algorithm is not necessarily the best. In order to obtain a more in depth and complete picture of the relative advantages and disadvantages of the developed algorithms, we compared them with respect to three additional features: first, their stability limit, which is defined as the maximum value of the ratio, u, of the time step dt over the space step dx, for which the numerical algorithm converges as the grid is refined; second, the CPU time requirements of each of the algorithms for the same grid spacings (dt =5 × 10 − 4, dx = 2 × 10 − 3, u=0.25). This CPU time provides an absolute measure of the computational work that the algorithm requires. However, the ultimate goal is to obtain an algorithm that can provide very accurate approximations in the smallest possible CPU time. Thus, the CPU time of greater importance, is that required to obtain accurate

approximations. This was the third and most decisive feature used to compare the algorithms to each other. In order to quantify what an ‘accurate approximation’ is, a certain measure of accuracy is required. As such, we used a measure of the distance of a given approximation from the converged solution depicted in Figs. 1 and 3. Let n˜c(0.5,tj ) denote the converged number density function at time tj and at the middle point of the space grid. Let also n˜dt,dx (0.5,tj ) denote the numerical approximation at time tj and at the middle point of the space grid obtained with grid spacings dt and dx. Finally, if J is the total number of points in time where we compare the converged and the numerical results, we define the distance from the converged result as: J

ddt,dx = % n˜c (0.5,tj )− n˜dt,dx (0.5,tj ) .

(37)

j=0

It was found that when the value of the distance, d, drops below the value of 0.1, the numerical solution is virtually indistinguishable from the converged result at any point in time and space. Thus, we are interested in finding the CPU time required by each algorithm in order for the value of d to drop below the value of 0.1. The results of these comparisons w.r.t. stability, CPU time required for the same grid spacings and CPU time required for accurate approximations, for both test problems and for both cases of equal and unequal partitioning, are shown in Table 2. Based on the results shown, one can make the following remarks: 1. The LF2 –6 algorithm is the least stable, whereas the RK4 –10 the most stable time-explicit scheme for both test problems and for both cases of equal and unequal partitioning. 2. Numerical stability is not influenced by the mathematical nature of the birth term (integral or functional), since the stability limit is the same for equal and unequal partitioning for every time-explicit algorithm and for both test problems. 3. The stability limit for the first test problem is approximately 2.8 times higher than that for the second test problem for all time-explicit numerical algorithms. This value is approximately equal to the inverse of the ratio of the maximum values for the single-cell growth rate for the two test problems (ln 2/2 and 1 for the first and second test problems, respectively). 4. Despite the fact that the hybrid algorithm requires a smaller CPU time than the LW2–2 scheme for the same grid spacings in the case of equal partitioning, the situation is reversed in the case of unequal partitioning for both test problems. This is due to the fact that the hybrid scheme requires, by construction (see Mantzaris et al., 1999), three evaluations of the birth term, as opposed to the LW2–2, which requires only two.

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

5. The RK4 –10 algorithm is the most computationally expensive time-explicit algorithm, whereas the LF2– 6 algorithm is the least expensive for the same grid spacings. This is expected, since the RK4–10 scheme requires four right-hand side evaluations per time step, while the LF2– 6 scheme requires just one. 6. Despite the fact that the RK4– 10 is the most computationally expensive time-explicit scheme in absolute terms, it is the most accurate and most stable from all the rest of the developed time-explicit algorithms. As a result, operating at the highest possible value of u, the algorithm can be applied with very coarse grid spacings and obtain very accurate approximations. Thus, for both test problems and for both cases of equal and unequal partitioning, the RK4 –10 scheme requires the smallest CPU time for accurate approximations. In particular, for the first test problem, this is achieved in 0.16 s for both equal and unequal partitioning. For the second test problem, convergence is achieved in 0.27 s for equal partitioning and 0.16 s for unequal partitioning. The RK4 – 10 algorithm is consistently more efficient

1431

than the APC4 – 8 algorithm, although the differences are very small. 7. The time-implicit CN2–10 algorithm performs relatively well only for the linear problem and the case of unequal partitioning. In all other cases, the CPU time required for accurate approximations is higher than that of the least efficient time-explicit scheme (hybrid). When comparing the performance of the CN2 –10 scheme with that of the hybrid scheme, it becomes clear that the more stable and accurate method is not necessarily the most efficient. 8. The CPU time required for the same grid spacings (4th column) for all time-explicit algorithms is generally higher in the case of unequal partitioning than in the case of equal partitioning. However, the CPU time requirements for accurate approximations (5th column) in the two cases are of the same order of magnitude. This is due to the properties of the solution to the continuous problem for equal and unequal partitioning. In the latter case, the converged solution shown in Figs. 1 and 3 seems to be such that it allows coarser grid spacings. On the

Table 2 Comparison between finite difference schemes (single-variable case) Scheme

Accuracy rank

Stability

CPU time (s) for dt= 5×10−4, dx=2×10−3

CPU time (s) for accurate approximations

Value

Rank

Value

Rank

Value

Rank

1.70 3.66 1.20 2.25 5.25 8

5 3 6 4 2 1

13.51 15.43 9.89 21.37 42.35 :30 h

2 3 1 4 5 6

1.15 1.00 0.44 0.22 0.16 1.42

5 4 3 2 1 6

Test problem 1 — unequal partitioning Hybrid 6 1.70 LW2–2 5 3.66 LF 2–6 4 1.20 APC 4–8 2 2.25 RK 4–10 1 5.25 CN2–10 3 8

5 3 6 4 2 1

19.64 17.12 12.10 21.25 41.19 :30 h

3 2 1 4 5 6

2.64 0.77 0.40 0.20 0.16 0.60

6 5 3 2 1 4

Test problem 2 — equal partitioning Hybrid 6 LW2–2 5 LF 2–6 4 APC 4–8 2 RK 4–10 1 CN2–10 3

0.61 1.46 0.42 0.82 1.92 8

5 3 6 4 2 1

14.39 18.51 10.98 24.40 48.01 :90 h

2 3 1 4 5 6

4.01 3.68 0.85 0.36 0.27 31.81

5 4 3 2 1 6

Test problem 2 — unequal partitioning Hybrid 6 0.61 LW2–2 5 1.46 LF 2–6 4 0.42 APC 4–8 2 0.82 RK 4–10 1 1.92 CN2–10 3 8

5 3 6 4 2 1

19.71 17.17 12.13 21.35 41.47 :90 h

3 2 1 4 5 6

1.82 0.44 0.36 0.17 0.16 2.69

5 4 3 2 1 6

Test problem 1 — equal partitioning Hybrid 6 LW2–2 5 LF 2–6 4 APC 4–8 2 RK 4–10 1 CN2–10 3

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1432

contrary, in the case of equal partitioning, in order to capture the periodic or oscillatory behavior of the solution, finer grid spacings are required, thus resulting in the observed result. Based on the observations made above and the results presented in Figs. 4 and 5, it is clear that the RK4 –10 and APC4 – 8 schemes are the two most efficient schemes. Despite the fact that, in absolute terms, the APC4 – 8 scheme is approximately twice as fast as the RK4 –10 scheme (two right-hand-side evaluations per time step versus four), the latter is more than twice as stable and is more accurate w.r.t. space. As a result, the RK4 –10 can be thought as the most efficient, although at such low CPU times as those presented in the 5th column of Table 2, significant differences cannot be observed. Moreover, despite the low CPU time requirements per time step for the LF2– 6 scheme, its performance is much worse than that of RK4– 10 and APC4 –8 schemes due to its significantly lower stability limit. The LW2 –2 and hybrid algorithms are obviously inferior due to their decreased accuracy. However, the LW2 –2 algorithm is clearly better than the hybrid scheme due to its improved stability properties. Finally, although the time-implicit CN2–10 scheme is unconditionally stable, it is less efficient than the time-explicit algorithms due to its increased CPU time requirements, which are even more pronounced in the case where cell growth is coupled with substrate consumption.

5.2. The two-6ariable case Using the conclusions drawn from the single-variable simulations, we considered the expansion of the two most efficient algorithms (RK4– 10 and APC4 – 8) in the case where the physiological state space is two-dimensional (r =2). Fig. 6 shows the results of the application of the RK4 –10 and APC4 – 8 schemes for the first test problem, whereas Fig. 7 corresponds to the second test problem. All simulations were performed using the same grid spacings (dt =10 − 3, dx1 =dx2 =5 × 10 − 3). No significant change in the numerical solution for either scheme was observed for finer grid spacings. Notice that for both test problems, both schemes converge to the same solution for the biomass and substrate concentrations, the middle point of the number density function and the state distribution function w.r.t. the first state variable, defined as follows: N(x1,t)=

&

1

N0 (x1,x2,t)dx2.

(38)

0

Thus, N(x1,t) is the average (in the x2 dimension) cell density of cells with content x1 at time t. The average in the x2 dimension number density function that is shown in Figs. 8 and 10 as a function of time is defined as follows:

n(x1,t)=

& && 1

0

1

0 1

N0 (x1,x2,t)dx2 .

(39)

N0 (x1,x2,t)dx2dx1

0

Due to the fact that the single-cell growth rates and single-cell division probability density functions were taken to be exactly the same in both space dimensions, the results are the same for the average in the x1 dimension cell density and number density function of cells with content x2 at time t (N(x2,t), n(x2,t), respectively) and for this reason are not shown. The numerical solutions correctly predict the perfect periodic behavior of the number density function in the case of batch growth and constant substrate concentration (test problem 1). Furthermore, for the second test problem, both algorithms predict a sustained oscillatory behavior for the number density function, similar to the single-variable case, while they both agree with the predictions of the Monod model for the biomass and substrate concentrations. The properties of the two algorithms when applied to the two test problems are summarized in Table 3. First, notice that for both test problems, the stability limits of both schemes drop approximately by a factor of two, when compared to the single-variable case. Moreover, similar to the single-variable case, the ratio of the stability limits for the two test problems is the inverse of the ratio of the maximum of the total single-cell growth rate (sum of the two single-cell growth rates) for the two test problems. Based on these observations, we found that the stability limit of the RK4–10 and APC4–8 algorithms, for any of the two test problems and for either the single- or the two-variable case, can be more generally expressed by the following two inequalities:

 

r

n n

max % ri (x,S) uRK4 − 10 5 1.8 x,S i = 1 r

max % ri (x,S) uAPC4 − 8 5 0.78. x,S i = 1

(40) (41)

Substitution of the parameter values for the single-cell growth rates given in Table 1, and the stability results shown in Tables 2 and 3 in Eqs. (40) and (41), shows that Eqs. (40) and (41) are approximately satisfied for any case. Thus, the RK4–10 scheme is more than twice as stable than the APC4–8 scheme. Despite the fact that the RK4– 10 scheme is approximately twice more computationally expensive for the same grid spacings (dt =10 − 3, dx1 = dx2 = 5× 10 − 3) (see 2nd column of Table 2), the higher stability limit and the increased accuracy, make the RK4–10 scheme a more efficient algorithm than the APC4–8, as the 3rd column of Table 3 shows. In particular, for the first (linear) test problem, the RK4–10 algorithm produces an approximation virtually indistinguishable from the converged result in just 6.1 s using dt=5.55× 10 − 2 and dx1 =

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1433

Fig. 6. Two-variable case. Illustration of convergence. Test problem 1. Simulation results for dt =1 × 10 − 3 and dx = 5 ×10 − 3. Solid line, open circles: the Runge – Kutta 4 –10 scheme. Dashed line, filled circles: the Adams – Predictor– Corrector 4 – 8 scheme. (a) Total biomass concentration; (b) middle point of the number density function; (c) the number density function and w.r.t. the first state; (d) the state distribution w.r.t. the first state as functions of time.

dx2 =2.27×10 − 2, while the APC4 – 8 scheme achieves convergence in 11.5 s with dt =2. × 10 − 2 and dx1 = dx2 =2× 10 − 2. For the second (nonlinear) test problem, very accurate solutions are obtained in similar CPU times. The RK4 – 10 scheme requires 17.1 s (dt = 2×10 − 2 and dx1 = dx2 =2.174 × 10 − 2), whereas the APC4 –8 scheme achieves convergence in 32.4 s (dt = 7.692×10 − 2 and dx1 =dx2 =2 ×10 − 2). We have also applied the most efficient RK4 –10 scheme for the nonlinear single-cell kinetics corresponding to the Williams model (test problem 3). The results of the simulation for dt =5 ×10 − 3 and dx1 =dx2 = 5× 10 − 3 (CPU time: 31.2 min) are shown in Fig. 8. The results of the integration of the continuum ODE

model consisting of equations Eqs. (27)–(30) are also shown for comparison purposes. Notice that, as expected by theory, there is a perfect agreement in the prediction of the two models for the biomass and substrate concentrations. Notice also that the dynamics of the two number density functions defined as in Eq. (39) are oscillatory. This is a consequence of the equal partitioning assumption. Contrary to the first two test problems where the dynamics of the number density functions are identical in the two space dimensions, in the case of the Williams model, the two number density functions have different dynamics due to the fact that the two single-cell growth rates are also different (see Eqs. (22) and (23)).

1434

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

5.3. The three-6ariable case The two-variable simulations have clearly shown the superiority of the RK4–10 algorithm over the rest of the developed algorithms. Very accurate solutions can be obtained for the two test problems by operating on the stability limit and using less than 50 grid points in each space dimension. Moreover, from the information obtained from the single- and two-variable simulations, it was also possible to extract an approximation of the stability limit of the RK4– 10 algorithm. We have expanded the RK4– 10 algorithm in order to solve the three-variable problem. Using the estimate in Eq. (40), we predicted the stability

limits for the two test problems. For both test problems, we used 51 grid points in each of the three space dimensions. Based on the stability limit predictions, the time steps were chosen to be dt= 3.125× 10 − 2 for the first test problem and dt =1× 10 − 2 for the second test problem. The CPU times required for the two simulations were 32.8 min for the first test problem and 111.7 min for the second test problem. The results are shown in Figs. 9 and 10 for the first and second test problems, respectively. Notice that as in the lower dimensional problems, the numerical scheme produces the expected results for both test problems, thus verifying again the validity of the RK4 –10 numerical algorithm.

Fig. 7. Two-variable case. Illustration of convergence. Test problem 2. Simulation results for dt = 1 × 10 − 3 and dx =5 ×10 − 3. Solid line: Monod model predictions for the biomass and substrate concentrations. Solid line, open circles: the Runge – Kutta 4 – 10 scheme. Dashed line, filled circles: the Adams – Predictor– Corrector 4 –8 scheme. (a) Total biomass and substrate concentrations; (b) middle point of the number density function; (c) the number density function and w.r.t. the first state; (d) the state distribution w.r.t. the first state as functions of time.

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1435

Fig. 8. Two-variable case. Test problem 3. Simulation results for dt = 5 × 10 − 3 and dx= 5 ×10 − 3. Solid line: Williams model predictions for the biomass and substrate concentrations. Dashed line: The Runge – Kutta 4 – 10 scheme. Dashed line, filled circles: the Adams – Predictor– Corrector 4 –8 scheme. (a) Total biomass concentration; (b) total substrate concentration; (c) and (d) the number density function and w.r.t. the first n(x1,t) and second n(x2,t) states averaged in the x2 and x1 space dimensions, respectively.

Table 3 Comparison between finite difference schemes (two-variable case) Scheme

Accuracy rank

Stability

CPU time (min) for dt= 1×10−3, dx=5×10−3

CPU time (s) for accurate approximations

Value

Rank

Value

Rank

Value

Rank

Test problem 1 — equal partitioning APC 4–8 2 RK 4–10 1

1.15 2.6

2 1

81.5 134.3

2 3

11.5 6.1

2 1

Test problem 2 APC 4–8 2 RK 4–10 1

0.41 0.96

2 1

88.3 147

2 3

32.4 17.1

2 1

1436

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

6. Summary and conclusions In this work, we have developed several finite difference algorithms for the numerical solution of the generalized multi-variable cell population balance equation coupled with the equations describing consumption of nutrients. All algorithms are general in the sense that they can be applied for any set of intrinsic physiological state functions (i.e. single-cell growth rates, single-cell division rates, single-cell rate of substrate consumption and partition probability density function). The developed algorithms were compared to each other w.r.t. stability, accuracy and computational speed. The ultimate criterion used to evaluate the efficiency of each

algorithm was the CPU time required in order to obtain approximations sufficiently close to the converged result. It was shown that the time-explicit algorithms are in general far more efficient than the time-implicit algorithms. Despite the fact that the specific time-implicit algorithm implemented was unconditionally stable, it requires the solution of a linear or nonlinear system of equations at each time step. As a result, the CPU time requirements are significantly increased when compared to those of the time-explicit algorithms especially in the case where the problem is nonlinear. This clear superiority of the time-explicit algorithms can be attributed to the mild nature of the nonlinearities that appear in the

Fig. 9. Three-variable case. Test problem 1. Simulation results for dt = 3.125 × 10 − 2 and dx=2 ×10 − 2 using the Runge– Kutta 4 – 10 scheme. (a) Total biomass concentration; (b) middle point of the number density function; (c) the number density function and w.r.t. the first state; (d) the state distribution w.r.t. the first state as functions of time.

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

1437

Fig. 10. Three-variable case. Test problem 1. Simulation results for dt=1 ×10 − 2 and dx =2 ×10 − 2 using the Runge – Kutta 4 – 10 scheme. (a) Total biomass and substrate concentration; (b) middle point of the number density function; (c) the number density function and w.r.t. the first state; (d) the state distribution w.r.t. the first state as functions of time.

problem under consideration, which do not result in severe stability restrictions. Consequently, some of the developed time-explicit algorithms are not just very accurate but have also high stability limits. Thus, the unconditional stability of the time-implicit algorithm does not offer a significant advantage over the timeexplicit schemes. By applying formulae of different order of accuracy for the approximation of the first spatial derivatives in conjunction with leapfrog, multi-step and one-step Runge–Kutta techniques for the time integration, it was seen that the higher the order of accuracy of the first spatial derivative approximation, the lower the

stability limit of the overall algorithm. Moreover, the higher the order of accuracy of the time derivative approximation in multi-step algorithms, the lower the stability limit is. The opposite holds for one-step Runge– Kutta schemes. Thus, for multi-step schemes, high accuracy in time integration results in a lower stability limit, while it leaves the CPU time requirements the same. However, for one-step Runge–Kutta schemes, a high accuracy in time integration results in higher stability limits. However, it leads to higher CPU time requirements as well, due to the higher number of right-hand side evaluations per time step.

1438

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

Stability is also influenced by the choice of numerical boundary conditions. Several different options were considered. The most sensitive family of schemes is the leapfrog family. The best result for these schemes is obtained by completely avoiding the use of any numerical boundary condition. This was achieved by using lower-order approximations of the first spatial derivatives close to the boundaries, which results in a reduction of the overall accuracy of the scheme. The same idea was also found to be the most effective in the case of the higher-order multi-step schemes. On the contrary, the application of the containment conditions as numerical boundary conditions to all non-existent grid points appearing in the spatial derivative approximations close to the boundaries yields the higher stability limit for all one-step Runge– Kutta algorithms. The most efficient algorithm was found to be the 4th-order Runge–Kutta algorithm with a 10th order approximation of the first spatial derivatives. The high accuracy and stability properties of this algorithm enabled the simulation of cell population balance formulations containing up to three state variables in a Pentium II 450 MHz PC in tractable CPU times. This algorithm can achieve convergence in 0.16 s for the linear problem and in 0.27 s for the nonlinear problem in the single-variable case. Furthermore, the CPU time requirements for accurate approximations are comparable for the equal and unequal partitioning cases, despite the fact that in the unequal partitioning case, the evaluation of an integral term per right-hand side evaluation per time step is required. This similarity in CPU time requirements is due to the fact that in the unequal partitioning case, the solution characteristics are such that they allow coarser grid spacings than in the equal partitioning case. In the two-variable case, convergence is achieved in 6.1 s for the linear problem and in 17.1 s for the nonlinear problem, while in the three-variable case, the algorithm requires 32.8 min and 111.7 min for the linear and nonlinear problems, respectively. The single-variable simulations revealed that the RK4– 10 algorithm provides a major improvement over the previously developed hybrid scheme, which was the starting point of this work. Moreover, the RK4– 10 scheme was shown to be the most accurate and the most stable from all the time-explicit algorithms developed. Based on the simulation results, a general formula for the stability limit that holds for the test problems considered and for any dimension of the physiological state space was derived. It must be noted, however, that the use of a 10th-order accurate approximation of the first spatial derivatives and a 4th order accurate approximation of the time derivatives requires the corresponding smoothness from the solution to the continuous problem. In general, whether the solution to the continuous problem is smooth enough or not, cannot be known a priori. Thus,

depending on the specific application of the cell population balance model, it may be that the upper bound in the order of accuracy of the derivative approximations is imposed not by a compromise between accuracy and numerical stability but by the smoothness of the solution to the continuous problem. Acknowledgements We thank the Graduate school of the University of Minnesota for awarding N.V.M a dissertation fellowship. Financial support by the National Science Foundation through the grant NSF/CTS-9624725, NSF/BES-9708146 and NSF/EES-9319380 is also gratefully acknowledged. Appendix A. Derivation of the boundary conditions (Eq. (8)) We will derive these conditions for the single-variable case and for a single-substrate. The zeroth moment of the state distribution function Nt (t) (cell density), increases due to birth and decreases due to division and exit from the bioreactor, while it remains unaffected from growth. Based on these observations, the dynamics of Nt (t) must be described by the equation: dNt (t) = dt

& & & xmax

xn,min

n

xmax

G(y,S)p(x,y,S)N(y,t)dy dx

2

x

xmax

G(x,S)N(x,t)dx − DNt (t).

−

(A1)

xmin

The first term of the right-hand side expresses the total rate of cell number increase due to birth, while the second expresses the total rate of loss due to division. Taking the zeroth moment of Eq. (6) over the entire one-dimensional physiological state space, we obtain: dNt (t) + dt

& &

( [r(x,S)N(x,t)]dx + xn, min(x xma x

xmax

G(x,S)N(x,t)dx + DNt (t)

& &

xmin

xmax

=

xmax

2

xn,min

n

G(y,S)p(x,y,S)N(y,t)dy dx.

x

(A2)

By comparing Eqs. (A1) and (A2), it is obvious that in order for the cell population balance equation to satisfy Eq. (A1) imposed by the physics of the problem, the following condition needs to be satisfied:

&

( [r(x,S)N(x,t)]dx = 0 xn, min(x or r(xn, min,S)N(xn, min,t)= r(xmax,S)N(xmax,t). xmax

(A3)

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

Moreover, since biomass must be conserved at cell division, the dynamics of the biomass concentration Nb(t) should only by affected by the growth rate and the dilution rate. Thus, the following equation should hold: dNb(t) = dt

&

+

r(x,S)N(x,t)dx −DNb(t).

(A4)

xn,min

& & & &

( [r(x,S)N(x,t)]dx (x

xma x

x

xn, min xmax

xG(x,S)N(x,t)dx +DNb(t)

x min

xmax

=

n

xmax

x 2

xn,min

G(y,S)p(x,y,S)N(y,t)dy dx

x

(A5)

By applying the properties in Eqs. (3) and (5) of the partition probability density function, it can be shown that the division term in Eq. (A5), which expresses the rate of biomass loss due to division, should cancel with the birth term in Eq. (A5) expressing the rate of biomass gain due to birth, i.e.

&

xmax

xG(x,S)N(x,t)dx

& &

xmin

xmax

=

xmax

n

G(y,S)p(x,y,S)N(y,t)dy dx.

x 2

xn,min

x

Furthermore, integrating by parts the second term in Eq. (A5) and applying the above equality expressing conservation of biomass at cell division, Eq. (A5) yields: dNb(t) + dt =

&

&

( [xr(x,S)N(x,t)]dx xn, min(x xma x

xmax

r(x,S)N(x,t)dx −DNb(t).

(A6)

xn,min

By comparing Eqs. (A4) and (A6), we see that in order for the cell population balance equation to satisfy Eq. (A4), the following condition should hold:

&

the physics of the problem: (a) cell division does not change the population mass; (b) cell growth does not change the population number. References

xmax

Taking the first moment of the cell population balance Eq. (6), we obtain: dNb(t) + dt

1439

( [xr(x,S)N(x,t)]dx =0 xn, min(x or xn, minr(xn, min,S)N(xn, min,t) =xmaxr(xmax,S)N(xmax,t). (A7) xmax

In order for Eqs. (A3) and (A7) to hold simultaneously and since xn,min " xmax, the containment conditions in Eq. (8) for the single-variable case: r(xn, min,S)N(xn, min,t) = r(xmax,S)N(xmax,t) =0 should be satisfied. Thus, we see that the containment conditions force the solution of the cell population balance Eq. (6) to satisfy two simple facts imposed by

Bailey, J. E., & Ollis, D. F. (1986). Biochemical engineering fundamentals (second ed.). New York: McGraw-Hill. Behnken, D. W., Horowitz, J., & Katz, S. (1963). Ind. Eng. Chem. Fundam., 2, 212. Block, D. E., Eitzman, P. D., Wangensteen, J. D., & Srienc, F. (1990). Slit scanning of Saccharomyces cere6isiae cells: quantification of asymmetric cell division and cell cycle progression in asynchronous culture. Biotechnol. Prog., 6, 504 – 512. Canuto, C., Hussaini, M. Y., Quarteroni, A., & Zang, T. A. (1988). Spectral methods in fluid dynamics. Berlin: Springer. Chiu, C. (1990). A numerical method for nonlinear age dependent population models. Differential and Integral Equations, 3 (4), 767 – 782. Christofides, P. D., & Daoutidis, P. (1997). Control of nonlinear distributed parameter systems: recent results and future research directions. In J. C. Kantor, C. E. Garcia, & B. Carnaham, Chemical process control — V, AIChE symposium series, vol. 93 (pp. 316, 302 – 306). Diekmann, O., Heijmans, H. J. A. M., & Thieme, H. R. (1984). On the stability of the cell size distribution. J. Math. Biol., 19, 227–248. Eakman, J. M., Fredrickson, A. G., & Tsuchiya, H. M. (1966). Statistics and dynamics of microbial cell populations. Chem. Eng. Prog., 62, 37 – 49. Fredrickson, A. G., Ramkrishna, D., & Tsuchiya, H. M. (1967). Statistics and dynamics of procaryotic cell populations. Math. Biosci., 1, 327 – 374. Fredrickson, A. G. (1999). Some unsolved problems in population balance modeling. In AIChE annual meeting, Dallas, Texas, Oct. 31 – No6. 5 (p. 1999). Fredrickson, A. G. & Mantzaris, N. V. (submitted for publication). A new set of population balance equations for microbial and cell cultures. Chem. Eng. Sci. Gerald, C. F., & Wheatley, P. O. (1989). Applied numerical analysis (fourth ed.). Reading, MA: Addison-Wesley. Godin, F. B., Cooper, D. G., & Rey, A. D. (1999a). Numerical methods for a population-balance model of a periodic fermentation process. AIChE J., 45 (6), 1359 – 1364. Godin, F. B., Cooper, D. G., & Rey, A. D. (1999b). Development and solution of a cell mass population balance model applied to the SCF process. Chem. Eng. Sci., 54 (5), 565 – 578. Greiner, G., & Nagel, R. (1985). Growth of cell populations via one-parameter semigroups of positive operators. In Mathematics applied to science (pp. 79 – 105). New York: Academic Press. Gustaffson, B., Kreiss, H. O., & Sundstro¨ m, A. (1972). Stability theory of difference approximations for mixed initial-boundary value problems, II. Math. Comput., 26, 649 – 686. Gyllenberg, M., & Webb, G. F. (1990). A nonlinear structured population balance model of tumor growth with quiescence. J. Math. Biol., 28, 671 – 694. Hatzis, C., Fredrickson, A. G., & Srienc, F. (1997). Cell-cycle analysis in phagotrophic microorganisms from flow cytometric histograms. J. Theor. Biol., 186 (2), 131 – 144. Hatzis, C., Srienc, F., & Fredrickson, A. G. (1995). Multistaged corpuscular models of microbial growth: Monte Carlo simulations. Biosystems, 36, 19 – 35. Hjortso, M. A., & Bailey, J. E. (1983). Transient responses of budding yeast populations. Math. Biosci., 63, 121 – 148. Hjortso, M. A., & Nielsen, J. (1994). A conceptual model of autonomous oscillations in microbial cultures. Chem. Eng. Sci., 49 (8), 1083 –1095.

1440

N.V. Mantzaris et al. / Computers and Chemical Engineering 25 (2001) 1411–1440

Hjortso, M. A., & Nielsen, J. (1995). Population balance models of autonomous microbial oscillations. J. Biotechnol., 42, 255 – 269. Kim, M.-Y. (1996). Galerkin methods for a model of population dynamics with nonlinear diffusion. Numerical Methods for Partial Differential Equations, 12, 59–73. Kim, M.-Y., & Park, E. J. (1995a). Mixed approximation of a population diffusion equation. Comput. Math. Applic., 30 (12), 23 – 33. Kim, M.-Y., & Park, E. J. (1995b). An upwind scheme for a nonlinear model in age-structured population dynamics. Comput. Math. Applic., 30 (8), 5 –17. Kostova, T. V. (1988). Numerical solutions of a hyperbolic differential-integral equation. Comput. Math. Applic, 15 (6-8), 427 – 436. Kostova, T. V. (1990). Numerical solutions to equations modelling nonlinearly interacting age-dependent populations. Comput. Math. Applic., 19 (8 – 9), 95 –103. Kostova, T., & Marcheva, M. (1989). Numerical solutions to the Gurtin– MacCamy equation. Math. Balkanica, 3, 265–277. Kromenaker, S. J., & Srienc, F. (1991). Cell-cycle-dependent protein accumulation by producer and nonproducer murine hybridoma cell lines: a population analysis. Biotechnol. Bioeng., 38, 665 – 677. Lax, P. D., & Wendroff, B. (1960). Systems of conservation laws. Communications on Pure and Applied Mathematics, 13, 217 – 237. LeVeque, R. J. (1992). Numerical methods for conser6ation laws (second ed.). Basel: Birkha¨ user. Liou, J. J., Srienc, F., & Fredrickson, A. G. (1997). Solutions of population balance models based on a successive generation approach. Chem. Eng. Sci., 9, 1529–1540. Mantzaris, N. V., Liou, J.-J., Daoutidis, P., & Srienc, F. (1999). Numerical solution of a mass structured cell population balance model in an environment of changing substrate concentration. J. Biotechnol., 71 (1 – 3), 157 –174. Monod, J. (1942). Recherches sur la croissance des cultures bacteriennes. Paris: Hermann. Natarajan, A., & Srienc, F. (1999). Dynamics of glucose uptake by single Escherichia coli cells. Metabol. Eng., 1, 320–333. Osher, S. (1969a). Stability of difference approximations of dissipative type for mixed initial-boundary value problems, I. Math. Comput., 23, 335 – 340.

Osher, S. (1969b). Systems of difference equations with general homogeneous boundary conditions. Trans. Am. Math. Soc., 137, 177 – 201. Ramkrishna, D. (1979). Statistical models of cell populations. Ad6. Biochem. Eng., 11, 1 – 47. Ramkrishna, D. (1985). The status of population balances. Re6. Chem. Eng., 3, 49 – 95. Ramkrishna, D. (1994). Toward a self-similar theory of microbial populations. Biotechnol. Bioeng., 43, 138 – 148. Ricthmyer, R. D., & Morton, K. W. (1967). Difference methods for initial-6alue problems. New York: Wiley. Roels, J. A. (1983). Energetics and kinetics in biotechnology. Amsterdam: Elsevier. Srienc, F. (1993). Flow cytometry in biotechnology: potential and limitations. ECB6. In Proceedings of the 6th European congress on biotechnology, Florence, Italy, June 13 – 17. Srienc, F., & Dien, B. S. (1992). Kinetics of the cell cycle of Saccharomyces cere6isiae. Biochemical Engineering VII: Annals of the New York Academy of Sciences, 665, 59 – 71. Strikwerda, J. C. (1989). Finite difference schemes and partial differential equations. London: Chapman & Hall. Subramanian, G., & Ramkrishna, D. (1971). On the solution of statistical models of cell populations. Math. Biosci., 10, 1–23. Sulsky, D. (1994). Numerical solution of structured population balance models. II Mass structure. J. Math. Biol., 32, 491 – 514. Sweeney, P. J., Srienc, F., & Fredrickson, A. G. (1994). Measurement of Unequal DNA partitioning in Tetrahymena pyriformis using slit-scanning flow cytometry. Biotechnol. Prog., 10, 19 – 25. Trucco, E. (1965a). Bull. Math. Biophys., 27, 285. Trucco, E. (1965b). Bull. Math. Biophys., 27, 449. Tsuchiya, H. M., Fredrickson, A. G., & Aris, R. (1966). Dynamics of microbial cell populations. Ad6. Chem. Eng., 6, 125 – 206. Williams, F. M. (1967). A model of cell growth dynamics. J. Theor. Biol., 15, 190. Young, D. M., & Gregory, R. T. (1988). A sur6ey of numerical mathematics, vol. 1. New York: Dover Publications, Inc. Zachmanoglou, E.C., & Thoe, D.W. (1986). Introduction to Partial Difference Equations with Applications, New York: Dover Publications, Inc.