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
integrodifferential
equations.
Analysis
even more complex
situation
of coupled
integrodifferential
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 birthanddeath
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 birthanddeath 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, l23 (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 integrodifferential 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 integrodifferential 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 transitionprobability 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 birthanddeath 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 vectorvalued 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. DISiRIEUTION
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)
J1
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,)=
lexp
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 cosmicray 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?Wk
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 tdistribution (for example, see Owen [6]). A general flow chart for a fixed number of simulations S, which is selfexplanatory, is shown in Fig. 1. We consider the application of the algorithm to two processes: the first, a simple birthanddeath process including the detailed age distribution of individuals, and the second, a cell massdistribution 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 BIRTHANDDEATH
9
GROWTH
particulate
system
with
PROCESS
We consider a birthanddeath 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 stxoj)
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(hCl)fh(x+Y)+N
[(A+
_~e(““)‘{e“Xr?v+e~AY
N,&(“P)‘h
f2(XYY,t)=<
}],
$$(y 6
i=

xoi

+Ap j=*
x
(294
x<
(29b)
t)eN
I
h eP[e(XP)(tX)l]
2 2 _ i=r
p)e2(“‘)‘X(X+Y)
NO 2 ~(yxoit)eW il
S(xxx,it)6(yxojt)[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(hP)r(l  emh),
x
(31a)
Noe(hC)‘(lee“‘)+eP’N(xtt,O),
x>t
Plb)
[1+$(Z)
P, (Gt)P,
(YJ)] x
~[~,(x,t)(leF)+P,(y,t)(ler)l,
Wa) [
P2(XXY~ t>= [
l+~(~)]a,cx,t)~,(t,t)~[~,(x.i)(le~)
PI (x>t)+
1 x_ll.
1
EP,(x,t)h(le‘“) N,
1
1,
x
Wb)
mea[N(xt,O)+N(Yt,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(lePr 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,)=lexp[(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 twodimensional, 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 confidencelimit 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
BIRTHAN&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 birthanddeath process values of the expected cumulative
with age distribution. Estinumber as a function of
14
B. H. SHAH 6( l
I
I
!
SIMPLE BIRTHANDDEATH )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 birthanddeath 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% confidencelimit bands about the estimated values, as shown in Fig. 4. The agedistribution 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 massdistribution model for cell populations, due to Eakman et al. [l].
SIMULATION
OF MICROBIAL
32oc
I
I
SIMRE
I
I
BIRTHANDDEATH
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
birthanddeath
of the secondorder
/
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 MichaelisMenten 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 daughtercell 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 daughtercell
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,)=lexp 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(lR)=
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)
i1
’
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,)315($)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,,+Pm’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.
Massdistribution
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.
Massdistribution
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.
Massdistribution
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). Growthrate 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 ifold integration of ithorder 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 growthrate 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, 230240 (1949). D. G. Kendall, An artificial realization of a simple “birthanddeath” 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.
249284. 6 7 8 9 10 11 12 13
D. B. Owen, Handbook of Statistical Tables, AddisonWesley, Reading, Mass., 1962. A. Ramakrishnan, Stochastic processes relating to particles distributed in a continuous infinity of states, Proc. Camb. Philos. Sot. 46, 595602 (1950). A. Ramakrishnan and S. K. Srinivasan, On age distribution in population growth, BuU. Math. Biophys. 20, 289303 (1958). D. Ramkrishna and J. D. Borwanker, A puristic analysis of population balance. I, Chem. Eng. Sci. 28, 14231435 (1973). D. Ramkrishna and J. D. Borwanker, A puristic analysis of population Chem. Eng. Sci. 29, 17111721 (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, l23 (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,
521537 (1970). E. Trucco, Mathematical models for cellular systems. The Von Foerster Parts I and II, BUN. Math. Biophys. 21, 285304, 449471 (1965).
equation.