Monte Carlo simulation of microbial population growth

Monte Carlo simulation of microbial population growth

Monte Carlo Simulation of Microbial Population B. H. SHAH, Growth J. D. BORWANKER AND D. RAMKRISHNA* Indian Instituie of Technology, Kanpur, ...

1MB Sizes 0 Downloads 28 Views

Monte

Carlo

Simulation

of Microbial

Population

B. H. SHAH,

Growth

J. D. BORWANKER

AND

D. RAMKRISHNA* Indian Instituie of Technology, Kanpur, U.P. 208016, India Communicated

by Charles

J. Mode

ABSTRACT Mathematical models individual cells differing

of cell populations from one another

which recognize them to be segregated into in state and behavior lead to complicated

integro-differential

equations.

Analysis

even more complex

situation

of coupled

integro-differential

the behavior

of populations

by Kendall

for simulating

model has been extended cellular state and behavior. that the random

variables

of fluctuations

in this work to include The resulting algorithm to be generated

in small populations equations. which

obey

leads to the

An early technique the birth-and-death

the factor of variation in individual is probabilistically exact in the sense

have exactly

derived

distribution

functions,

and

involves no arbitrary discretization of the time interval, characteristic of simulation techniques in general. The simulated systems include populations distributed according to their age in a birth-and-death process and a batch culture of cells distributed according to their mass which grow and reproduce by binary fission.

INTRODUCTION Microbial populations consist of individual cells that differ from one another in their physiological states and consequently in their growth and reproductive processes in any particular environment. Mathematical models have attempted to account for this complexity by assigning a simple identifying structure such as cell mass or age to the individual cell and accounting for variation among individual cells through a statistical distribution function for which equations are obtained based on the concept of a

*Current address: School Lafayette, Indiana 47907. MATHEMATICAL 0 American

Elsevier

of

BIOSCIENCES Publishing

Chemical

Engineering,

31, l-23 (1976)

Company,

Inc., 1976

Purdue

University,

West

1

2

B. H. SHAH ET AL.

number balance. The penalty for this sophistication, even with such gross approximations of the physiological state of the cell, has been the solution of integro-differential equations which result from the model. Small populations add to the complexity by the presence of fluctuations whose calculation would require the solution of coupled integro-differential equations [9]. Solutions have been obtained for large populations, which are distributed according to their age [2, 141 or mass [l, 121. A mathematical model of the type just referred to is essentially a set of stipulations regarding the behavior of the individual cell. Thus a characteristic state such as age or mass may be defined for the individual cell, together with the kinetics of its growth and a transition probability function associated with the birth of new cells from the parent cell. The term “growth” in the present context refers to the rate of change of state, which for cell age is of course unity. The growth rate and the transition probability of cell division will in general depend on the cell state, time and the “environment”, which may be defined by concentrations of one or more essential nutrient substances. Also other events such as death of an individual may be accounted for through a transition-probability function for them. These stipulations, together with a knowledge of the initial state of the population and the environment, should suffice for the prediction of all future states. Kendall [4] considered an “artificial realization” of a simple birth-anddeath process in the following framework. The process involves the random appearance of new individuals and the disappearance of existing individuals with respective transition probabilities or “frequency functions” which continually change the strength of the total population, indicated by the total number. Thus a “birth” event will add one to the population while a “death” event will delete one from it. As is customary, during a time interval dt, multiple events will occur with probability of O(dt2). Kendall defines a “time interval of quiescence”, the period between successive events during which the number of individuals would remain constant. Obviously, this interval of quiescence is a random quantity because of the random nature of the events themselves. Kendall shows that the interval of quiescence has an exponential distribution with a coefficient parameter which depends on the number of individuals at the beginning of the interval, although this latter qualification depends upon the model for the transition probabilities. Thus if the size of the population is known at the beginning of the interval of quiescence, the size at the end of the interval will depend on whether the last event is a birth or a death. Given that a birth or a death has occurred, the probability of either event is readily obtained as the ratio of the corresponding transition probability to the sum of the two frequencies. An artificial realization of the process is now made possible by the generation of random numbers which represent successive

SIMULATION

OF MICROBIAL GROWTH

3

quiescence intervals and identify whether the event at the termination of each interval is a birth or a death. The random numbers to be generated should obey their governing probabilistic laws. Several such realizations can be obtained from the instant of reckoning of the process, from which all information about the process, probabilistic or average, must be calculable. The objective of this paper is to show that Kendall’s approach is easily extended to obtain artificial realizations of mathematical models of populations which incorporate the additional detail of individual variation in some characteristic state. The resulting algorithm is a powerful tool for the simulation of particulate systems in general. Consequently, we derive the algorithm for an abstract particulate system and subsequently discuss its application to examples of biological populations. THE SIMULATION

ALGORITHM

We consider a system of particles in which the individual particle is distinguished by a single scalar property X which can take on values between 0 and co. The extension to a vector-valued particle state will become evident in the course of the following discussion. We assume that the rate of change of particle state, denoted X, is a function of the instantaneous state. The situation when X may depend on some environmental variables is, in principle, not difficult to accommodate, although computationally burdensome. The procedure for this generalization is indicated at a subsequent stage, but the derivation of the algorithm and the examples of application will exclude this feature. The particles in the system are subject to events which we will call particulate events, the occurrence of any of which will result in a change in the number of particles. Thus the breakage of a particle into smaller particles or fission of a cell into daughter cells are particulate events which increase the number of particles, while the agglomeration of two particles to form a single particle diminishes the number. The latter type, which necessarily involves the participation of more than one particle (i.e., biparticulate events), will not be treated here. Thus, this paper will deal with monoparticulate events such as cell division or cell death. Random particulate events are modeled through frequencies or transition probabilities. We suppose for the sake of generality that there are m monoparticulate events, the ith event possessing a frequency bi, which is assumed to depend on the particle state. Since environmental variables are excluded for the present, their influence, if any, on the frequencies b, would require the generalization outlined subsequently. At this point we recognize that the model may require additional stipulations pertaining to the states of newborn particles. The state of a newborn particle may depend on the state of the parent particle and may have a probabilistic distribution. Thus a daughter cell formed from the

B.H.SHAH

4

ETAI..

division of a mother cell of mass (say) m’ may have any mass between 0 and m’. If, however, age is the characteristic state of a cell, newborn cells will have age zero. Following KendalI 141,we first seek to determine the probability distribution of the interval of quiescence following an instant at which the state of the system is known, i.e.. the number of particles along with the states of individual particles is known. DISi-RIEUTION

OF QUIESCENCE

TIME

INTERVALS

The probability distribution to be determined is conditional on a known state of the system a? time t. Suppose there are N particles at time I with states _x,.x2,. . ,xN. We will refer to the particle of state x, as thejth particle. The state of thejth particle at any time 7 subsequent to I is given by the solution of the differential equation

with the initial condition that x=x, at T=O; the solution is most conveniently represented by x(! + T.x,,I). We denote the state of the system at time I by A, and define the following conditional probability: P (TIA,) = Pr{ interval of quiescence > TIA,) The interval of quiescence will be greater than 7 + do, if it is greater than 7 and during the time interval (r? 7 + dr) none of the N particles then existing underwent any of the m particulate events. The mathematical equivalent of the preceding statement is

1

P(r+dTlA,)=P(TjA,)

I-

2 5 4(X(t+T,Xj.t))dT].

!=I

Rearranging the terms in Eq. (2). dividing obtains the differential equation

i P(rjAt)

must

satisfy

by d7 and letting

dT+0,

one

2 5 bi(x(t+7,xj,t))

-$P(TiA,)= -P(rlA,) Further

(2)

J-1

i=l ,‘I

the obvious

(3) 1

initial

condition

P(O]A,)=

I,

SIMULATION

5

OF MICROBIAL GROWTH

subject to which the solution

of Eq. (3) is readily found. Thus

which is an exponential distribution with a variable parameter. The cumulative distribution function F(rlA,), which is the probability that the quiescence interval is less than r, is then given by

F(rlA,)=

l-exp

1

- 2 2 i=i

JTb,(,(r+~‘,~l,l))d7r

j=l

O

.1

(5)

Next we compute the probability that, given that a particulate event has occurred after a quiescence interval T following time t, thejth particle was involved in the ith event. This is conveniently done by defining a twodimensional indicator variable u which takes on the value (ij) for the event just mentioned to be true. It is readily seen that

Pr{ u = (W% T } = m

bi(x(t+ T,Xj,t)) Ncrj

(6)

Z: C bi(X(t+ T,Xj,t)) . i=l

j=l

When the quiescence interval T is specified along with the identity of the particle and the identity of the event it has undergone, the state of the system at t+ T is exactly calculable from that at time t, if the states of any newborns are precisely known. SIMULATION

OF SYSTEM

EVOLUTION

The particulate system under study evolves in time starting from a given initial condition. The preceding discussion has established that any particular realization of the system is represented by a sequence of random variables

(T,,u,)>

(Tz>dr...>(T,,,u,,)

(7)

which represents the state of the system until time t =Cy=,Tj with changes in the number of particles occurring at T,, T, + T2, T, + T,+ T3,. . . , T, + T2+ . . . + T,,. Of course, the tacit assumption has been made that the states of newborns are precisely known. The situation when this assumption does

6

B. H. SHAH

ET AL.

not hold will be treated through an example. Since the probability distributions of the random variables in the sequence are known, artificial realizations are possible by generating the random variables on a digital computer. The random variables are sampled sequentially until a specified value of the final time is just exceeded. A sufficient number of repetitions starting from the initial condition would then allow us to estimate the expected values of various quantities of interest and fluctuations about them. If the initial condition is known only in a probabilistic sense, then this would require the generation of specific realizations of the random initial conditions. We will not concern outselves here with the method of generating random variables, since it is available elsewhere. The simulation strategy should now be evident. The state of the system at any time t, represented by A,, starting from that at t = 0, is obtained as follows. From A,, generate the interval of quiescence T,. Update the state of the jth particle to i= 1,2 )...) N,. x = x ( T,, x/7O), Compute the event frequencies at t = T, and simulate the event by generating the indicator variable u. At this point one may need to generate additional random variables, especially in assigning the “birth states” of newborns. Since samples of a random variable with an arbitrary cumulative distribution function are obtainable by the methods in the preceding section, the states of newborns can be generated when a probabilistic model is available for their distribution. ESTIMATION

OF SYSTEM

PROPERTIES

Essentially, the number of simulations must be large enough to provide consistent estimates of the various quantities of interest which describe the system. Some of the quantities evaluated in this work are as follows. If n(x, t)dx represents the number of particles at time t with states between x and x+ dx, then N (x, t), the number of particles with states less than x, is given by N(x,t)=fn(x’,t)dx’.

Since N (x, t) is in general a random

(8)

quantity,

,8, (xpt) = E

one may obtain

[N (x, 41.

(9)

One may also write PI as P, (x,t)=

Lxf, (x’,t)dx’.

(10)

SIMULATION

7

OF MICROBIAL GROWTH

wheref,(x, t) is a product density of order 1, a concept first developed by A. Ramakrishnan [7] in the analysis of cosmic-ray showers; f,(x,t)dx is the expected number of particles between x and x + dx at time t. One may also obtain a function /?*(x,y, t) (x
which is alternatively

expressed in terms of the second product density

P2(KY, t>= J’,yyf2 (X’>Y’> t)dx'dy'.

fi as

(12)

,B2is thus associated with the calculation of fluctuations about mean values The particulate process defined in the preceding sections leads to integrodifferential equations in f,, fi etc. [9] which may be solved subject to appropriate initial conditions to obtain the values of /It and & However, the present aim is to compute their values from the results of direct simulation. Cumulative properties of the system derived from individual properties may also be evaluated. Let W(X) be the property associated with a particle of state x. Then define W(t) by

W(t)=

~mw(x)n(x,l)dx,

which is obviously a stochastic process. nature of the particle population,

Equivalently,

(13) due to the discrete

N(r)

W(t)=

z

w(x,).

(14)

j=l The expected value of W(t) and the variance

may be obtained

EW(1)=Jmw(x>f, (x,t)b

(15)

0

+

Yv2(x)f, I

0

(x,t)dx-

as in [9]:

V

mw(x)fi

0

2

(-Gt)dx

1 .

B. H. SHAH ET AL.

8

The formulae for estimating the quantities (9) (1 l), (15) and (16) from the simulation results are easily identified. Representing the estimated quantities by a circumflex on the symbol, we obtain for a total of S simulation runs the following results.

(17)

(18)

EW(t)= +

E Wj(l),

J=l

PW(t)=

&

,g [ u:(+~~wl*.

(21)

The confidence intervals of estimates were also determined using a confidence limit of 95%. The confidence interval for EW, for example, is given by l?W-k
where k = c[ PW/S]‘/* with c a constant which depends on the number of degrees of freedom S- 1 and can be found from standard tables of the Student’s t-distribution (for example, see Owen [6]). A general flow chart for a fixed number of simulations S, which is self-explanatory, is shown in Fig. 1. We consider the application of the algorithm to two processes: the first, a simple birth-and-death process including the detailed age distribution of individuals, and the second, a cell mass-distribution model due to Eakman, Fredrickson and Tsuchiya [l]. The first is distinguished from the second by the existence of analytic solutions with which the simulation results may be compared. The second problem has been solved earlier by the method of weighted residuals [ 121.

SIMULATION

OF MICROBIAL

FIG. I. A general flow chart for simulation of a random growing particles and estimation of system properties.

THE BIRTH-AND-DEATH

9

GROWTH

particulate

system

with

PROCESS

We consider a birth-and-death process for a population distributed according to age. We refer to Kendall [3] for details and directly proceed to the differential equations in terms of the product density functions f,(x,r) and f>(,x,y, r), where x and y now refer to age. If X(x) and p(x) are respectively the birth and death frequencies for an individual of age x, the differential equations in .f, and fi are easily derived:

(23)

B. H. SHAH

10

ET AL.

The differential equations contain no birth frequency, because birth has been assumed to leave the parent unaltered while adding an individual of age zero to the population. The possibility of multiplets requires some modifications, and A. Ramakrishnan and Srinivasan [8] have accounted for twins. The simulation procedure in this paper needs no special consideration to account for this effect, which is done naturally and readily. Equations (22) and (23) must be subject to the boundary conditions (24)

(25) If the initial population consists of No individuals with ages xo,,xo2, then the initial conditions for f, and fiare given by f,

(x,0)=

2 stx--oj)

XONo7

(26)

j=l

NO

f2 Equations arbitrary on letting simulation dence of increased.

NO

c 8(Y - XOj>, j=l

(XYY> 0) = Ix f3tx- XOi) i=l

i#j.

(27)

(22) and (23) are readily solved subject to (24) through (27) for frequencies h(x) and p(x). The algebra is considerably simplified h and p be independent of age, i.e., constants. Of course, the procedure is in no way complicated by allowing for age depenA and CL, only the number of simulations required may be For h and p constant we obtain r

Noe(X-“)‘-ti 3

O
(284

Pb) No2h2,2(h-Cl)f--h(x+Y)+N

[(A+

_~e(“-“)‘{e-“X-r?v+e-~-AY

N,&(“-P)‘-h

f2(XYY,t)=<

}],

$$(y 6

i=

-

xoi

-

+A-p j=*

x
(294

x<

(29b)

t)e-N

I

h -e-P[e(X-P)(t-X)-l]

2 2 _ i=r

p)e2(“-‘)‘-X(X+Y)

NO 2 ~(y-xoi-t)e-W i-l

S(x-xx,i-t)6(y-xoj-t)[email protected]


t
(29~)

SIMULATION

OF MICROBIAL GROWTH

11

In (29a) and (29b) the first term dominates that one may verify that

for large enough values of N,, so

f2 (xJJ, t)=:fr (x, t)f, (v, t).

(30)

The quantities ,B,(x,t) and &(x,y,t) can be obtained by substituting (29a), (29b) and (29~) into Eqs. (10) and (12). The result is

PI (x, t) =

(28)

Noe(h-P)r(l - emh),

x
(31a)

Noe(h-C)‘(l-ee-“‘)+e-P’N(x-tt,O),

x>t

Plb)

[1+$(Z)

P, (-Gt)P,

(YJ)] x
-~[~,(x,t)(l-e-F)+P,(y,t)(l-e-r)l,

Wa) [

P2(XXY~ t>= [

l+~(~)]a,cx,t)~,(t,t)-~[~,(x.i)(l-e-~)

PI (x>t)+

1 x_ll.

1

EP,(x,t)-h(l-e-‘“) N,

1

1,

x
Wb)

me-a[N(x-t,O)+N(Y-t,O)l

I

PI (t,t)+

kemZp’N(x-

We thus the model model. The discussion,

t,O)[N(y-

G

t,O)-

1 ( 11,

~P,(t,t)-A(l-e-Pr N,

t
4 (32~)

have exact theoretical values of the quantities of interest from with which could be compared the simulated results of the adaptation of the algorithm to the model is presented without since it is readily understood.

12

B. H. SHAH

The cumulative

distribution

function

for quiescence

intervals

F(r)A,)=l-exp[-(X+p)~N(f)l. A realization

of T in terms of the uniform

ET AL.

is (33)

random

variable

R is given by

The indicator variable u is two-dimensional, but a slight variation here is desirable from the point of view of convenience. First, we simulate the type of event by generating a random indicator variable I, such that I= 1 if a birth takes place and Z= - 1 if a death takes place. Clearly, Pr{Z=l[A,,T}=&,

(35)

Pr{Z=-liA,,T}=&.

(36)

The addition of Z to N automatically updates N, while the ages of the individuals present are updated by the addition of T. Birth adds an individual of age zero. When death occurs, we generate another indicator random variable u such that U= i if the ith individual dies. Clearly, Pr{u=iiA,,T,Z=-l}=&. The representative

RESULTS

The defined X Inax, at selected assumed Thus’

sequence

OF SIMULATION;

(37)

for the process is then

COMPARISON

WITH

THEORY

numerical results are presented in terms of dimensionless age as the ratio of the actual age x to the age of the oldest individual, the initial time. The numerical values of h and p were arbitrarily respectively. The initial condition was to be l/x,,, and 1/2x,,, to be N, particles spaced uniformly on the age interval (O,xma.J xai = Lx,,,, No

‘Note that in general

for a uniform

i=1,2,...,N,,.

initial distribution

[i.e., G(x) = x], one must space

the ith particle on xoi after generating the ith uniform random variable R, (between 0 and 1) and obtain xoi = Rix,,. The choice Ri = i/N, is essentially a “deterministic analog” of the uniform random variable. A similar procedure is indeed conceivable for arbitrary G(x).

SIMULATION

OF MICROBIAL

13

GROWTH

In the figures, the initial expected cumulative age distributions are shown as steps wherever convenient. For large N,, however, the cumbersome steps have been replaced by straight lines which enclose them. In all computations, the simulated results contained the analytical curves within the bands of 95% confidence limits. With increased number of simulations these confidence-limit bands narrowed considerably, giving very close agreement with the analytical result. Increasing the initial number of particles, N,, decreases the number of simulations required for convergence of the estimated results for the expected cumulative age distribution, ,L?,,to the analytical values p,. Figure 2 shows the simulation results along side the analytical values for N,= 10, obtained with 60 simulation runs. The agreement between the two is evident. Figure 3 shows the results for N,=25 with 40 simulation runs.

I6 SIMPLE No= S

16

BIRTH-AN&DEATH

iROCESS



IO

I

,=I

‘60

A = I/X,,, ,o, /xmax

,p = 0 5/x,,, = I/N,

14

12

6

05

IO x/x

FIG. 2. Simulation mated and theoretical dimensionless age.

max,

15

20

2.5

DIMENSIONLESS AGE

of a simple birth-and-death process values of the expected cumulative

with age distribution. Estinumber as a function of

14

B. H. SHAH 6( l-

I

I

!

SIMPLE BIRTH-AND-DEATH )r* I/Xmax xoi Ix,,,&

I

I

I

ET AL.

-

PROCESS

9r=0.5/x,,, i/N,,o=I, ,N,

No’25 5c ,-

s =40

40

s: x 30

20 0

IO

1

0 025

05

/

075 x/x max

FIG. 3. Simulation mated and theoretical dimensionless age.

I

1.0

1

1.25

DIMENSIONLESS

I

1

I.5

1.75

,

20

AGE

of a simple birth-and-death process values of the expected cumulative

with age distribution. Estinumber as a function of

Figure 4 shows the results for ,8, along side the analytical values & at two different times, plotted against x,/x,,, with x~/x,,,=~ as a parameter. The analytical values were always within the 95% confidence-limit bands about the estimated values, as shown in Fig. 4. The age-distribution model considered above does not pertain to microbial populations for which binary division of the mother cell leads to replacement by two daughter cells, both of age zero. The mathematical model and the simulation procedure are both altered slightly. The latter changes in regard to the updating of the population state when a birth event takes place. The modification is, however, easily made. Since our eventual interest in the paper is the simulation of microbial populations, we consider the following mass-distribution model for cell populations, due to Eakman et al. [l].

SIMULATION

OF MICROBIAL

32oc

I

I

SIMRE

I

I

BIRTH-AND-DEATH

x,_Jx~~~I/N~,

28OC

I= I,

1.5

GROWTH /

I

I

I

PROCESS A= I/xmm,r

, N,

=o5/x,o,

No=25 2400 /x,,,~~ =

; x. x-


I’

s =60 + =2

2000

Anal.

05

S,m.

-

0

1600 G L. x

n

I200

< Q 600

400

r 1

I

I

I

02

04

06

08

0 0

/

,

1.0

1.2

xl’%lax~DIMENSIONLESS FIG. 4.

Simulation

mated and theoretical dimensionless age.

of a simple values

birth-and-death

of the second-order

/

I

1

1.4

16

1.8

AGE

process cumulative

with age distribution. number

Esti-

as a function

of

A CELL MASS DISTRIBUTION MODEL FOR MICROBIAL POPULATIONS In the interest of conserving space, we refer the reader for details to the works of Eakman et al. [l] and of Subramanian et al. [13]. Briefly, the model considers a population of microorganisms distributed according to their mass, in which individual cells grow by appropriating substrate from the environment and multiply by binary division. The cell growth follows Michaelis-Menten kinetics, which after accounting for endogenous metabolism yields

m.

--pc

(38)

I

The transition

probability

of division,

b(m),

is given by

exp[ -(?)‘I b(m)=

Km

&VG

(39) erfc( y)



B. H. SHAH ET AL.

16 The daughter-cell mass m’ is taken

mass distribution for cells formed to be the one used by Subramanian

from a mother et al. [13]:

cell of

(40)

The cumulative

distribution

function

~mp(m”.m’)dm’r=

The population mass distribution

in a batch given by

for daughter-cell

lO( z)‘-

microbial

mass

15( $)4+6(

propagator

is given by

2)‘.

evolves

(41)

from

an initial

(42)

where_/,

is the initial

number

density

of cells. The initial

EC, =

mmjo(m)dm=N,a(n+l)

biomass

concentra-

tion is given by

SIMULATION

I0

(43)

OF THE MODEL

For the simulation,

we neglect

the effect

of a changing

environment-

one that in principle is accommodated by the algorithm, without the hazards of substantially increased computations. tion balance equation is given by

ah Cm,4

at

+ $i;

(m,t)ti=2_fwb(m’)p(m,m’)l, m

although not The popula-

(m’,t)dm’

- Wm).h (m,t) with the initial

(44)

condition (45)

h (wO>=h Cm). To obtain the quiescence single cell from given initial

interval distribution conditions using

dm -=fiil;ikm, dr

we trace

m (0) = m,,

the growth

of a

(46)

SIMULATION

OF MICROBIAL GROWTH

where k is a constant, from some arbitrary in = m(t + r,m,, t) is

17

since the environment is constant. If 7 = 0 is measured time t, then the solution of (45) represented by

m =

mie kr.

(47)

Let the system at time t be A,=[m,,i=1,2 where N(t) is the total number function for the quiescence interval

F(rlA,)=l-exp I

,..., N(t)], of cells. The cumulative following t is given by

- C /‘b(m(t+r’,m,,t))dr’ i=* O

A sample value of the quiescence

interval

distribution

.

(48)

I

T is then found from

h’(r)

-ln(l-R)=

Substitution result

x J$(nl(t+7,mi,t))d7. ;=, 0

of (39) into (49) together

3

(49)

with some algebra

erfcj miek~-mc)

will produce

the

=,_R (50)

i-1



erfc( F)

which is the source from which T is evaluated given a sample value of R. The identity of the dividing cell is governed by an indicator random variable u such that u = i if the cell of mass m, divides. Clearly, b(miekr) Pr{ u = ilA,, T } = N(r) x b(m,e”‘) i=l Further,

if one of the daughter

(51) .

cells has mass M, the other will have mass

18 mi-

B. H. SHAH M,

ET AL.

where m, is the mass of the mother cell. From (41) we have

Pr{M~m~A,,T,u=i}=lO

( 111 m,)3-15($)4+6(35.

(52)

The sample values of u and M may be obtained from successive generations of the uniform random variable as indicated earlier. Some specific computational aspects that simplified the simulation procedure have been discussed by Shah [ 1 l] and will be omitted here. Table 1 presents the numerical values of the constants. TABLE

1

Numerical Values of Constants Used in the Cell Mass Distribution Model Quantity

Numerical

k

0.1787 hr-’ 0.4734 hr -

SIMULATION

m,

3x10-‘2g

E

3ti

a

2x10-‘*g

x lo-”

Value

’ g

RESULTS

The results of the simulation stand by themselves, since no analytical solutions are possible with which to compare. The results obtained, however, have the same features as the approximate solutions obtained by Subramanian and Ramkrishna [ 121 through the method of weighted residuals. Simulations were carried out with a small initial population, since the average behavior is independent of the population size, as shown by Ramkrishna and Borwanker [9, lo]. Thus for No=20 only 10 simulation runs were sufficient to obtain consistent estimates of the expected quantities ,8, and total biomass. Of course, the confidence levels narrowed considerably with increased number. Figure 5 shows the expected number of cells and the expected biomass concentration as a function of time with the first initial condition. Confidence limits are indicated for the expected number. The distribution 8, are shown in Fig. 6, plotted at different times. Figure 7 shows the asymptotic nature of the normalized cumulative mass distribution j?,/gN, which is attained with the onset of exponential growth.

SIMULATION

OF MICROBIAL I

40

I

I

= N,,+P-m’a

ho)

f

N,

19

GROWTH 1

1

a = 2 x lo-”

= 20

u s

=I0

EN (1) WITH 95% CONFIDENCE LIMITS

32

- 51

47

24 01

/

5

/


/

/

/

/

/’

z E F -43

,/’

/

/

/

B

/’

5

/’ :

39

z 0 m t,

OL 0

I

I 0.4

I

TIME, hr

1 08

I

I

1

1.2

J35 Ij 6

VARIATION OF TOTAL NUMBER AND BIOMASS CONCENTRATION

FIG. 5.

Mass-distribution

tion at various

model.

Variation

of total number

and biomass

concentra-

times.

CONCLUSIONS The simulation technique presented in this paper has been shown to be an effective method of handling models of particle populations which recognize the segregated nature of these systems, and in particular models of population growth. There are no strong constraints of a mathematical nature on obtaining answers, as there are in the exact or approximate solutions of the governing population balance equations. Besides, fluctuations characteristic of smaller populations can also be obtained from the simulations; the alternative route of obtaining them from successive product density equations presents immense difficulties. The examples in this paper did not account either for explicit time dependence of any quantity or for the effects of environment on the population. These effects are in principle readily accommodated. Thus, for example, k may depend on time and some environmental variable (say JJ),

20

B. H. SHAH

40 h,d

f

ET AL.

=+,-m/a a=ZxlO-'*gm

N,

= 20

s

2

k

= 0.1787

IO hr-’

32

24

m,

MASS, grn X IO” 2

I

MASSDISTRIBUTION FIG. 6.

Mass-distribution

mass at different

model.

3

AT DIFFERENT

Variation

of cumulative

4

TIMES.

number

as a function

of

times.

which calls for rewriting dx’ = dr

_

Eq. (1) as

B (x’,y, t + 7),

j=

1,2 )...) N(t),

(14

for thejth cell, satisfying the initial condition xj= xj at T =O, where xi is the state at time t. The reciprocal effect of the population on the environment may be described by a differential equation such as

-4 = ?(x’,x2 )...) xN(‘);y,t+7) dr

supplemented by the solution of Eqs. (la) should be substituted for the functions x(t +

(lb)

initial condition T = 0, y =y(t). The simultaneous and (lb) determines the functions xj(t + T) which into Eqs. (2) through (6) as respective replacements T,x/, t). It should then be clear that the extraction of z

SIMULATION

OF MICROBIAL

21

GROWTH

f cm r oj =Bs -m/o [email protected]

~‘2xlO“~gm

NO = 20 s = IO k = 0.1787 hi’

1.0

t,llr 0

0.2

0 .

0.4 0.8

08

II

1.2

h

16

0.6 s z ZW . = 04

.

0.2

m,MASS,gmXlO” 0

ASYMPTOTIC NORMALIZED MASS DISTRIBUTION FIG.

function

7.

Mass-distribution

model.

Asymptotic

normalized

mass

distribution

as a

of mass.

from Eq. (12) for given values of R may well involve prohibitive computation unless some ad hoc simplifying calculation schemes are concocted in each case.

NOMENCLATURE LA TIN SYMBOLS a

‘4, b C

E F

F,

Variable coefficient of exponential distribution. State of population at time t. Transition probability of monoparticulate event. Concentration of biomass c,: substrate concentration in Eq. (48). Statistical expectation. Cumulative distribution function for quiescence interval. Cumulative distribution function for random variable Z.

B. H. SHAH ET AL.

22 _h G Z k ni m m, N n P

Pr R r S T t u V W W

X R X

Y P Y Z Z

GREEK

Product density of order i (i= 1 and 2). Cumulative distribution function for the initial state of the population. Indicator random variable in Eqs. (41) and (42). Growth-rate constant for cell growth, as in Eq. (56). Also coefficient of linear birth frequency. k,: constant in Eq. (48). Growth rate of cells. Cell mass. Critical mass above which division probability is very high. Total number of particles. N,=value of N at t =O. Number density of particles. Also an integer in Eq. (52). Probability function for quiescence interval defined below Eq. (1). Probability. Uniform random variable in the interval (0,l). Specific values of R. Also radius of cell in Eq. (48). Number of simulations. Quiescence interval. Time. Indicator random variable. Statistical variance. Macroscopic property of particulate system derived from w. Particle property. Particle state. Particle growth rate. Value of particle state. Environmental variable. Rate of change of environment. Value of particle state different from x. Also value of environmental variable. Random variable. Value of random variable Z. SYMBOLS

Parameter in Eq. (52). Cumulative quantity derived from i-fold integration of ith-order Pi product density. Dirac delta function. 6 Parameter in division frequency of cell in Eq. (49). Birth frequency. : Death frequency. Also parameter in cell growth rate in Eq. (48). P pcL,,{ Parameters in cell growth-rate equation (38). Quiescence time after time t. 7

a

SIMULATION

OF MICROBIAL

23

GROWTH

REFERENCES

4

J. M. Eakman, A. G. Fredrickson and H. M. Tsuchiya, Statistics and dynamics of microbial cell populations, Chem. Eng. Prog. Symp. Ser., No. 69, 62, 3749 (1966). A. G. Fredrickson and H. M. Tsuchiya, Continuous propagation of microorganisms, Am. Inst. Chem. Eng. J. 9, 459468 (1963). D. G. Kendall, Stochastic processes and population growth, J. Roy. Stat. Sot., Ser. B 11, 230-240 (1949). D. G. Kendall, An artificial realization of a simple “birth-and-death” process, J. Roy.

5

Srat. Sot., Ser. B 12, 116126 (1950). J. Moshman, Random number generation,

1 2 3

Computers, Vol. II, (A. Ralston

in Mathematical Methods jar Digital

and H. S. Wilf, Eds.), Wiley,

New York,

1967, pp.

249-284. 6 7 8 9 10 11 12 13

D. B. Owen, Handbook of Statistical Tables, Addison-Wesley, Reading, Mass., 1962. A. Ramakrishnan, Stochastic processes relating to particles distributed in a continuous infinity of states, Proc. Camb. Philos. Sot. 46, 595-602 (1950). A. Ramakrishnan and S. K. Srinivasan, On age distribution in population growth, BuU. Math. Biophys. 20, 289-303 (1958). D. Ramkrishna and J. D. Borwanker, A puristic analysis of population balance. I, Chem. Eng. Sci. 28, 1423-1435 (1973). D. Ramkrishna and J. D. Borwanker, A puristic analysis of population Chem. Eng. Sci. 29, 1711-1721 (1974). B. H. Shah, A simulative and analytical study of particulate systems,

Ph.D.

II,

thesis,

Indian Institute of Technology, Kanpur, 1975. G. Subramanian and D. Ramkrishna, On the solution of statistical models of cell populations, Math. Biosci. 10, l-23 (1971). G. Subramanian, D. Ramkrishna, A. G. Fredrickson and H. M. Tsuchiya, On the mass

14

balance.

distribution

model

for microbial

cell populations,

Bull. Math. Biophys. 32,

521-537 (1970). E. Trucco, Mathematical models for cellular systems. The Von Foerster Parts I and II, BUN. Math. Biophys. 21, 285-304, 449471 (1965).

equation.