- Email: [email protected]
Lag time in microbe growth as an age-dependent branching process with two phase types
B U L L E T I N OF
MATHEMATICALBIOLOCY VOLUME 40, 1978
L A G T I M E IN M I C R O B E G R O W T H AS AN A G E - D E P E N D E N T BRANCHING PROCESS WITH TWO PHASE TYPES
. M A R K J. CHRISTENSENa n d R. SHONKWlLER School of Mathematics, G e o r g i a Institute of Technology, Atlanta, Georgia 30332, U.S.A.
We study a two-type, age-dependent branching process in which the branching probabilities of one of the types may vary with time.Specifically this modification of the Bellman Harris process starts with a Type I particle which may either die or change to a Type 11 particle depending upon a time varying probability. A Type 11 particle may either die or reproduce with fixed probabilities but may not return to a particle of Type 1. In this way the process models the lag phenomenon observed in microbe growth subsequent to transfer to a new culture medium while the organism is adapting to its new environment. We show that if the mean reproduction rate of Type 1I particles exceeds 1, then the population size grows exponentially. Further the extinction probability for this process is related to that of the Bellman Harris process. Finally the governing equations are solved for several choices of the growth parameters and the solutions are graphically displayed showing that a wide variety of behavior can be modeled by this process.
Introduction.
A familiar biological phenomenon is the appearance of a growth lag when microorganisms are transferred to a new nutrient medium. Accordingly a mathematical descriptio n of the population size as a stochastic process does not conform to the well-known Age-Dependent Branching Process due to Bellman and Harris (1952). In the following we offer a biological explanation of this growth lag and, based upon it, a mathematical model for this stochastic process. Our result is a Waugh modified Two-type Age-Dependent Branching Process which can be viewed as a process with evolution or as a process with immigration. Let Gi(t) be the life length distribution function for the phase type i (i = 1, 2) and Z(t) be the total number of organisms or particles at time t given 1 newly 347
348
MARK J. CHRISTENSEN AND R. SHONKWILER
born particle of phase type 1 at time 0. Denote by q,(t) (respectively p,(t)) the probability that if a particle in phase 1 (respectively 2) ends its life, it will be replaced by n particles in phase 2. Finally let M (t) denote the expected number of particles in phase 2 at time t given l newly born particle in phase 2 at time 0. Then the expected population size is given by the pair of equations
E(Z(t))=I-GI(t)+
/0
ml(U)M(t-u)dGl(U)
(0.1)
M(s)=I-G2(s)+
m2(w)M(s-w)dG2(w ),
(0.2)
and
where
tnl(b/)= ~ n=O
rlqn(bl)
and
m2 (w) = ~ kpk (w). k=O
Some features of this system are worthy of comment. The assumption that particles in both phase types are replaced by ones of type 2 only leads to the uncoupling of the system. Thus, having solved the second equation for M, the final solution E (Z (t)) is expressed as a definite integral of M. Another important feature is that the probability functions qn and Pk may be time-dependent so this model can be used to extend that of Waugh (1955). From the point of view of the biological phenomenon this capability is necessary in phase type 1. Finally a fixed minimum lag time 6 can be implemented in either phase, say phase 1, by setting G1 (t) = 0 for t<6 (cf. Mode, 1969).
1. Biochemical Background. In one of the classic experiments supporting the Operon Theory of gene regulation it was discovered that wild-type E. coli cells cultured on glucose have about 5 molecules of fl-galactosidase per cell while those cultured on lactose have about 5,000 per celt. This enzyme is not needed for the catabolism of glucose and in the presence of glucose the enzymes'
LAG TIME IN MICROBE G R O W T H
349
production is repressed by a chemical feedback mechanism with the result that energy and raw materials are thereby conserved. However, upon transfer to a nutrient medium containing only lactose this enzyme is now vitally needed to hydrolyze the glucosidic bond of the metabolite whereupon the hydrolyzed products may enter the usual glycolysis chain. The time required by E. coli to produce the necessary level of this enzyme at its reaction site is 1 to 2 minutes at physiological temperatures. This same time lag phenomenon is observed to a greater or lesser extent, ranging up to several hours, by most bacterial cells (certain mutants excepted) cultured on various nutrient pairs. It is observed by those population counting techniques which can distinguish viable cells from expired ones that an unusually large number of cells die during the "tooling-up" period. This is to be expected as the cells are no less vulnerable to their environment in the lag period and simultaneously are temporarily without a renewable energy source.
2. The Mathematical Model. We envision the population of the branching process to be composed of organisms of two different types or alternatively as a single organism subject to two regimes of growth. The first regime corresponds to the adaptation period and no reproduction is assumed here. The second corresponds to the normal growth of those organisms which have successfully adapted. No particles are assumed to be in this regime at time t =0, and no particles from this regime ever return to regime 1. In the derivation which follows, both regimes are treated as a Waugh modification of the BellmanHarris process in that the branching probabilities can vary with time. The second regime is further altered to allow for immigration from regime 1. However, in the balance of the paper only the branching probabilities in regime 1 are time varying corresponding to the fact that if an event occurs soon in regime 1, that event is most likely expiration, qo(t). In the same manner reproduction at the moment of adaptation, q,(t)> 0 (n > 2), is included in our model, but in the applications we assume adaptation occurs without simultaneous reproduction, qo(t)+ql (t)= 1 (t >0). Let X(t). and Y(t) be the number of particles in regimes 1 and 2 respectively at time t, and let F(s, t) and H(s, t) be the probability generating functions for the random variables Z(t)=X (t)+ Y(t) and Y(t) respectively,
F(s,t)= ~ s~P(Z(t)=z) z=O
and
H(s,t)= ~ s~P(T(t)-z).
(2.1)
350
M A R K J. C H R I S T E N S E N
A N D R. S H O N K W I L E R
Setting Qz (t) = P (Y (t) = z), the probability that there are z particles in regime 2 at time t given 1 particle is newly born in regime 2 at time 0, we see by the Cauchy product of power series that
Hk(s, t ) = L s~
~
Qi,(t) ... Qik (t).
il+...+ik=z
z=O
We recognize the second summation to be exactly the probability P(Y(t) = z 1Y(0-) = 0, Y (0 + ) = k ) that k particles enter regime 2 at time 0 and lead to z particles at time t. In fact each of the k particles, zh,..., Zrk,leading independent and identically distributed lives, results in ii..... ik particles respectively at time t and this number is z if il +... + ik = Z. Therefore
Hk(s't) = L s:P(Y(t)=z[Y(O )=O,Y(O+)=k).
(2.2)
--=0
Further, from Waugh (1955, p. 294), H satisfies the integral equation
H(s,t)=
/o
h(H(s,t-u),u)dG2(u)+s(1-G2(t)
),
(2.3)
where
h(s, u ) = L PJ(u)sj"
(2.4)
j=0
To find F(s, t) we compute the probabilities appearing in the sum. First consider z = 0, the probability of extinction at time t. This can occur either in regime 1 at any time u up to time t (qo(u)dGl(U)) or, if adaptation with k offspring has taken place at some time u less than t (qk (U) d G 1(u)) then extinction can occur in regime 2 during the time t - u ( P ( Y ( t ) = O [Y(u-)=0, Y(u + ) = k) = Q ~ ( t - u ) ) . Thus
P(Z(t)=O)=
I
t
qo(u)dGl(u) 0
+
Q~(t-u)qk(u)dGl(U) k=l
=
0
qk(u)Qko ( t - u ) d G l ( u ) . k=O
0
(2.5)
LAG TIME IN MICROBE GROWTH Next let z > l . E i t h e r X ( t ) = 1 o r X ( t ) = 0 s o
351
that
P(Z(t)=z)=P(X(t)--1, Y ( t ) = z - 1 ) + P ( X ( t ) = O , Y(t)=z).
(2.6)
Regarding the first term, no particle can appear in regime 2 ifX(t) = 1 and hence
• (P(X(t)=I)=I=GI(t), P(X(t)=I,Y(t)=z-1)=~O, for x > l .
for
z=l (2.7)
For the second term, the event X ( t ) = 0 , Y(t)= z can be decomposed into the time u of migration from regime l whereupon k ( k = l , 2 .... ) newly born particles appear in regime 2 and lead to z particles during time t - u. Therefore
P(X (t)=0, r ( t ) = z ) = It k~"1 P(Y(t)=ztY(u- )=0, 0
=
Y(u+)=k)qk(u)dGl(u).
(2.8)
Substituting equations (2.5) through (2.8) into (2.1) we obtain
F(s,t)=s(1 - G l ( t ) ) +
+
z=l
k=O
s:
0
Qko(t--u)qk(u)dal(u )
k=l
P(Y(t)=zlY(u-)=O, Y (u + ) = k )qk(u )dG 1(u ).
Or using (2.2),
F(s, t) = s(1 - G1 (t)) +
fi
Ok=O
qk(u)Hk(s, t -- u)dGl (u).
(2.9)
Equation (2.9) in conjunction with (2.3) constitutes the mathematical model. The various moments of Z(t) are easily derived from these equations. For example
E(Z(t))
OF(s,t) as s=l = l - G l ( t ) + it ~ qk(l,)kHk-l(1, t--u) ~H(Scs--U)'tA dOk=~
s~ 1dG1 (u).
352
MAR K J. C H R I S T E N S E N A N D R. S H O N K W I L E R
However H k - ~ ( 1 , t - u ) - - 1 and, denoting the partial derivative appearing under the integral sign by M ( t - u), this becomes equatio'n (0.1). At the same time, from (2.3), and (2.4),
M(t) =~H(s, t) (~S
= 1 - Gz(t) s= l
+
f
t ~
0 k=O
kPk(u)H k 1 ( 1 , t - u )
CH(s,t_u)
~S
=1
dGl(u),
which is equation (0.2).
3. The Behavior q[" the Moments oJ Z ( t ) J o r Large t. There are standard results for the asymptotic behavior of the Bellman-Harris process. In this section we obtain the modification of these results to the system (0.1) and (0.2). As might be expected, the original result for the unmodified Bellman Harris process will be used in this analysis. Therefore, we now state that result; LEMMA 1.
I f M is the solution of the integral equation t
I
M(t)= 1 -G(t)+m
M(t-s)dG(s), 0
where m > 1, d G ( s ) > 0 , G ( 0 + ) = 0 , G ( , ~ ) = 1 and G is not a step Junction all of whose steps occur at regularly spaced points, then M ~ n 1 exp(at) as t--* ::~, where a satisfies the equation 1= m
f oc~exp( - as)dG(s) 0
and na=(m-1)/(am z
sexp(-as)dG(s)). 0
F o r the p r o o f of this lemma we refer the reader to Feller (1966), BharuchaReid (1960) or Athreya and Ney (1972). In this last reference the reader will also find a discussion of the various cases that can occur if m 2 =< 1. We note, however, that their discussion will be transferable directly to our (0.2) since this is exactly the Bellman Harris case. We can now state our result.
LAG
TIME
IN MICROBE
GROWTH
353
Letm 2 > 1 andsupposeGzsatisfiesalltheconditioninLemmal. Then the solution of(0.1) will obey,for large t THEOREM 1.
E(Z(t)) ~ n'1 exp(at), where a is the root o[the equation
fo
1 =m2
exp(-as)dG2(s)
while 77'1is given by,
j
" :7.~,
tTrl ~-111
t?l 1 (s)exp( -
as)dG1 (s).
0
Proof Let c > 0 be given and then choose Tso that 1 - G I ( T ) < e approximation M(t) ~ t71 exp(at) is valid for t > T. Then, since
and the
t
E(Z(t))=I-G1(t)+
I
ml(s)M(t-s)dGl(s), 0
we see that for t > T t-- T
E(Z(t))=I-Ga(t)+
I
ml(S)M(t-s)dGl(s)
do t
+
I
ml(s)M(t-s)dGl(S).
Now since t is large the first term will be less than e, while in the first integral we see that the argument of M is t - s > t - ( t - T ) = T a n d hence M may be replaced by the asymptotic form nl exp(at). Doing this we see that,
f
t- T 0
ml(s)M(t-s)dGlls)~
ft
7"
ml(s)nlexp(a(t-s))dGl(s)
0 t
7'
ml(s)exp(-as)dGllS).
= n 1 exp(at) (}
354
MARK
J. C H R I S T E N S E N
AND
R. S H O N K W I L E R
As for the second integral, since ml ( s ) < K and t > T we'see that
f
'
ml(S)M(t-s)dGl
(s) <2Kna
exp(aT)(Gl(t)-Gl(T))
t-T
< 2 K n l exp(aT)e. Thus we see that, for t > T,
E(Z(t)),,-nl exp(d! )
I
t - 7"
ml(s)exp(-as)dG1 (s) 0
+ 2Knl exp(aTle ml(s)exp(-as)dGl(s)+2Ke, exp(a(T-t)) .
= n l exp(at) 0
Therefore, as t--, oo, the expression in braces becomes,
f
,- r m,
(s ) e x p ( - a s )d G ~(s ) + 2Ke exp(a( T - t ) )
0
f
oo
---,
ml(s)exp(-as)dGa(s)+O. 0
Thus the proof is completed.
4. Extinction Probability. In this section the branching probabilities q,(t) and p,(t) (n=0, 1.... ) may vary with time. Let p denote the probability that starting from 1 newly born particle in regime 1 at time 0, extinction of the colony occurs at some time t < ~ . p = lira P(Z(t)=O). t -+ 30
I f a is the extinction probability for the process H(s, t) describing regime 2 (cf. Athreya and Ney, 1972, p. 143), then
P= ~ ~rk k=0
qk(u)dGl(U). 0
(4.1)
L A G T I M E IN M I C R O B E G R O W T H
355
This is easily seen, for denoting by I t the indicator function of the interval [0, t], we may write equation (2.5) as P(Z(t)=O)=
I,qk(u)Q~(t-u)dGl(u). k=O
0
The integrand of this expression as a function of t is bounded by the constant 1 which is a dG1 integrable function of u. Hence Lebesque's Dominated Convergence Theorem applies and since Q~ (t - u) ~ ak as t ~ ov we obtain (4.1). 5. Sample Growth Curves. In this section we present a number of the population growth curves that result from the system (0.1) and (0.2) for certain choices of the distribution functions G1 and Gz and of the branching probabilities qi and pj. (ij = 0, 1, 2 .... ). Our principal interest here is to study the lag or initial region of the growth curve. Since this region is primarily affected by the choice of parameters for regime 1, we shall assume, as is often done, that the process in regime 2 is Markovian and thereby obtain a simple closed form solution for the population growth function M ( t ) in regime 2. Specifically the life-length distribution function G2 must be exponentially distributed
Gz(t)=l-e
~a+u)t, t > 0 ,
if it is assumed that future states depend only on the present state of the process in regime 2 and if the branching probabilities # Po-2+it
and
2 P2-2+it
are independent of time. But for this exponential distribution the integral equation (0, 2) has the solution M (t) = e I~- ,~t. Hence the overall population size is given by t
E ( Z ( t ) ) = 1 - G1 (t) +
S
e ~ - ~'~"- U)ql (u)dG1 (u).
0
In regime 1 we take the life length to be gamma distributed dGl(t) dt
fl~+l F(a+l)
t ~ e -/~t
(5.1)
356
MARK J. CHRISTENSEN AND R. SHONKWILER
with parameters e and ft. This is a very flexible distribution over the range of parameters capable of including, for example the Markovian case by choosing =0.
Moreover we incorporate the seldomly used but critically important authenticity feature of allowing the branching probabilities in regime 1 to vary with time. Specifically we put ql(t)=ct~e
-bt
and
qo(t)=l-q~(t)
with a, b, and c independent parameters satisfying a>0,
b>0,
and c ( a / b )a e a = m a x ql ( . ) =< l.
This permits the simulation of an "adaptation window" centered around the maximum point for ql namely t = a / b . TABLE I
Graph number
Density c~ fl
a
Adaptation probability b
c
1
0.1
0.367
10
160
*
2
0
1
2
32
*
3
0
0.3
2
8
*
4
2
1.0
0
0.9
*
5
0.1
0.733
0.1
0.8
*
6
0
4
0.1
0.4
*
*c is chosen so that the maximum value of ql {I ) is 1.
Description G a m m a life length, wide plateau at T; thin window at T/48 Exponential life length averaging T/3: m e d i u m window at T/48 Exponential life length averaging T; m e d i u m window at T/12 G a m m a life length : m e d i u m at T; m e d i u m window at 0 G a m m a life length, wide plateau at T/2; wide window at T/24 Exponential life length averaging T/12; wide window at 7"/13
LAG TIME IN MICROBE GROWTH
357
The graphs which follow were obtained by a computer graphing program evaluating the equation (5.1) at 40 abscissa points. We used the same parameter values in regime 2 for each curve, namely 2=0.316
and
#=0.017,
which corresponds to an expected life-length of T--3 hours and a doubling probability of 95 per cent. This data is based upon the literature (cf. Shonkwiler, 1977). 3--
2:--
6
5
o 13.
I O
I
I
I
2
I
I
I
I
3
4
5
6
Time
Figure 1.
As our objective is to display the spectrum of behavior possible in the initial region of the growth curve, we choose a variety of values for the parameters c~,/~, a, b and c toward that end. The possibilities can be descriptively enumerated as follows. The distribution is taken either as exponentially distributed (~ =0) with expected life length times, 1//~, taken to be T, 1/3T, or 1/12T, or as gamma distributed with a medium (c~=2) or wide (~=0.1) plateau and expected life length times, (c~+ 1 )//~, taken to be 1/2 T or T. The adaptation probability q 1(t) is taken to have an adaptation window (maximum) at 1/12T, 1/2T, 1/48Tor 0 which is either wide (a = 0.1 if not at 0) medium (a = 2 if not at 0) or thin (a = 10 if not at 0). Table 1 is a list of the parameter values used for the six graphs. Inspection of these growth curves shows that this model is capable of describing a wide variety of behavior depending upon the profile and relative
358
MARK J. CHRISTENSEN AND R. SHONKW1LER
p l a c e m e n t of the a d a p t i v e w i n d o w . In p a r t i c u l a r c u r v e 5 s h o w s a p r o n o u n c e d d e l a y b e f o r e the o n s e t of e x p o n e n t i a l g r o w t h .
LITERATURE Athreya, K. B. and P. E. Ney. 1972. Branching Processes. New York: Springer-Verlag. Bellman, R. and T. Harris. 1952. "On Age-Dependent Binary Branching Processes." Ann. Math., 55, 280-295. Bharucha-Reid, A. T. 1960. Elements of the Theory oJ Markot: Processes and 7 heir Applications. New York: McGraw-Hill. Feller, W. 1966. An Introduction to Probability Theory and its Applications, Vol. II. New York: John Wiley & Sons, Inc. Mode, C. J. 1969. "'Lag Time in Cell Division from the Point of View of the Bellman-Harris Process." Math. Biosci., 5,341-345. Shonkwiler, R. 1977. "Consequences of the Growth and Division Model for Mitochondria Assuming No Cellular Control." Bull. Math. Biol., 39, 613-618. Waugh, W. A. O'N. 1955. "An Age-Dependent Birth and Death Process." Biometrica, 42,291 306. RECEIVED 11-22-76
MATHEMATICALBIOLOCY VOLUME 40, 1978
L A G T I M E IN M I C R O B E G R O W T H AS AN A G E - D E P E N D E N T BRANCHING PROCESS WITH TWO PHASE TYPES
. M A R K J. CHRISTENSENa n d R. SHONKWlLER School of Mathematics, G e o r g i a Institute of Technology, Atlanta, Georgia 30332, U.S.A.
We study a two-type, age-dependent branching process in which the branching probabilities of one of the types may vary with time.Specifically this modification of the Bellman Harris process starts with a Type I particle which may either die or change to a Type 11 particle depending upon a time varying probability. A Type 11 particle may either die or reproduce with fixed probabilities but may not return to a particle of Type 1. In this way the process models the lag phenomenon observed in microbe growth subsequent to transfer to a new culture medium while the organism is adapting to its new environment. We show that if the mean reproduction rate of Type 1I particles exceeds 1, then the population size grows exponentially. Further the extinction probability for this process is related to that of the Bellman Harris process. Finally the governing equations are solved for several choices of the growth parameters and the solutions are graphically displayed showing that a wide variety of behavior can be modeled by this process.
Introduction.
A familiar biological phenomenon is the appearance of a growth lag when microorganisms are transferred to a new nutrient medium. Accordingly a mathematical descriptio n of the population size as a stochastic process does not conform to the well-known Age-Dependent Branching Process due to Bellman and Harris (1952). In the following we offer a biological explanation of this growth lag and, based upon it, a mathematical model for this stochastic process. Our result is a Waugh modified Two-type Age-Dependent Branching Process which can be viewed as a process with evolution or as a process with immigration. Let Gi(t) be the life length distribution function for the phase type i (i = 1, 2) and Z(t) be the total number of organisms or particles at time t given 1 newly 347
348
MARK J. CHRISTENSEN AND R. SHONKWILER
born particle of phase type 1 at time 0. Denote by q,(t) (respectively p,(t)) the probability that if a particle in phase 1 (respectively 2) ends its life, it will be replaced by n particles in phase 2. Finally let M (t) denote the expected number of particles in phase 2 at time t given l newly born particle in phase 2 at time 0. Then the expected population size is given by the pair of equations
E(Z(t))=I-GI(t)+
/0
ml(U)M(t-u)dGl(U)
(0.1)
M(s)=I-G2(s)+
m2(w)M(s-w)dG2(w ),
(0.2)
and
where
tnl(b/)= ~ n=O
rlqn(bl)
and
m2 (w) = ~ kpk (w). k=O
Some features of this system are worthy of comment. The assumption that particles in both phase types are replaced by ones of type 2 only leads to the uncoupling of the system. Thus, having solved the second equation for M, the final solution E (Z (t)) is expressed as a definite integral of M. Another important feature is that the probability functions qn and Pk may be time-dependent so this model can be used to extend that of Waugh (1955). From the point of view of the biological phenomenon this capability is necessary in phase type 1. Finally a fixed minimum lag time 6 can be implemented in either phase, say phase 1, by setting G1 (t) = 0 for t<6 (cf. Mode, 1969).
1. Biochemical Background. In one of the classic experiments supporting the Operon Theory of gene regulation it was discovered that wild-type E. coli cells cultured on glucose have about 5 molecules of fl-galactosidase per cell while those cultured on lactose have about 5,000 per celt. This enzyme is not needed for the catabolism of glucose and in the presence of glucose the enzymes'
LAG TIME IN MICROBE G R O W T H
349
production is repressed by a chemical feedback mechanism with the result that energy and raw materials are thereby conserved. However, upon transfer to a nutrient medium containing only lactose this enzyme is now vitally needed to hydrolyze the glucosidic bond of the metabolite whereupon the hydrolyzed products may enter the usual glycolysis chain. The time required by E. coli to produce the necessary level of this enzyme at its reaction site is 1 to 2 minutes at physiological temperatures. This same time lag phenomenon is observed to a greater or lesser extent, ranging up to several hours, by most bacterial cells (certain mutants excepted) cultured on various nutrient pairs. It is observed by those population counting techniques which can distinguish viable cells from expired ones that an unusually large number of cells die during the "tooling-up" period. This is to be expected as the cells are no less vulnerable to their environment in the lag period and simultaneously are temporarily without a renewable energy source.
2. The Mathematical Model. We envision the population of the branching process to be composed of organisms of two different types or alternatively as a single organism subject to two regimes of growth. The first regime corresponds to the adaptation period and no reproduction is assumed here. The second corresponds to the normal growth of those organisms which have successfully adapted. No particles are assumed to be in this regime at time t =0, and no particles from this regime ever return to regime 1. In the derivation which follows, both regimes are treated as a Waugh modification of the BellmanHarris process in that the branching probabilities can vary with time. The second regime is further altered to allow for immigration from regime 1. However, in the balance of the paper only the branching probabilities in regime 1 are time varying corresponding to the fact that if an event occurs soon in regime 1, that event is most likely expiration, qo(t). In the same manner reproduction at the moment of adaptation, q,(t)> 0 (n > 2), is included in our model, but in the applications we assume adaptation occurs without simultaneous reproduction, qo(t)+ql (t)= 1 (t >0). Let X(t). and Y(t) be the number of particles in regimes 1 and 2 respectively at time t, and let F(s, t) and H(s, t) be the probability generating functions for the random variables Z(t)=X (t)+ Y(t) and Y(t) respectively,
F(s,t)= ~ s~P(Z(t)=z) z=O
and
H(s,t)= ~ s~P(T(t)-z).
(2.1)
350
M A R K J. C H R I S T E N S E N
A N D R. S H O N K W I L E R
Setting Qz (t) = P (Y (t) = z), the probability that there are z particles in regime 2 at time t given 1 particle is newly born in regime 2 at time 0, we see by the Cauchy product of power series that
Hk(s, t ) = L s~
~
Qi,(t) ... Qik (t).
il+...+ik=z
z=O
We recognize the second summation to be exactly the probability P(Y(t) = z 1Y(0-) = 0, Y (0 + ) = k ) that k particles enter regime 2 at time 0 and lead to z particles at time t. In fact each of the k particles, zh,..., Zrk,leading independent and identically distributed lives, results in ii..... ik particles respectively at time t and this number is z if il +... + ik = Z. Therefore
Hk(s't) = L s:P(Y(t)=z[Y(O )=O,Y(O+)=k).
(2.2)
--=0
Further, from Waugh (1955, p. 294), H satisfies the integral equation
H(s,t)=
/o
h(H(s,t-u),u)dG2(u)+s(1-G2(t)
),
(2.3)
where
h(s, u ) = L PJ(u)sj"
(2.4)
j=0
To find F(s, t) we compute the probabilities appearing in the sum. First consider z = 0, the probability of extinction at time t. This can occur either in regime 1 at any time u up to time t (qo(u)dGl(U)) or, if adaptation with k offspring has taken place at some time u less than t (qk (U) d G 1(u)) then extinction can occur in regime 2 during the time t - u ( P ( Y ( t ) = O [Y(u-)=0, Y(u + ) = k) = Q ~ ( t - u ) ) . Thus
P(Z(t)=O)=
I
t
qo(u)dGl(u) 0
+
Q~(t-u)qk(u)dGl(U) k=l
=
0
qk(u)Qko ( t - u ) d G l ( u ) . k=O
0
(2.5)
LAG TIME IN MICROBE GROWTH Next let z > l . E i t h e r X ( t ) = 1 o r X ( t ) = 0 s o
351
that
P(Z(t)=z)=P(X(t)--1, Y ( t ) = z - 1 ) + P ( X ( t ) = O , Y(t)=z).
(2.6)
Regarding the first term, no particle can appear in regime 2 ifX(t) = 1 and hence
• (P(X(t)=I)=I=GI(t), P(X(t)=I,Y(t)=z-1)=~O, for x > l .
for
z=l (2.7)
For the second term, the event X ( t ) = 0 , Y(t)= z can be decomposed into the time u of migration from regime l whereupon k ( k = l , 2 .... ) newly born particles appear in regime 2 and lead to z particles during time t - u. Therefore
P(X (t)=0, r ( t ) = z ) = It k~"1 P(Y(t)=ztY(u- )=0, 0
=
Y(u+)=k)qk(u)dGl(u).
(2.8)
Substituting equations (2.5) through (2.8) into (2.1) we obtain
F(s,t)=s(1 - G l ( t ) ) +
+
z=l
k=O
s:
0
Qko(t--u)qk(u)dal(u )
k=l
P(Y(t)=zlY(u-)=O, Y (u + ) = k )qk(u )dG 1(u ).
Or using (2.2),
F(s, t) = s(1 - G1 (t)) +
fi
Ok=O
qk(u)Hk(s, t -- u)dGl (u).
(2.9)
Equation (2.9) in conjunction with (2.3) constitutes the mathematical model. The various moments of Z(t) are easily derived from these equations. For example
E(Z(t))
OF(s,t) as s=l = l - G l ( t ) + it ~ qk(l,)kHk-l(1, t--u) ~H(Scs--U)'tA dOk=~
s~ 1dG1 (u).
352
MAR K J. C H R I S T E N S E N A N D R. S H O N K W I L E R
However H k - ~ ( 1 , t - u ) - - 1 and, denoting the partial derivative appearing under the integral sign by M ( t - u), this becomes equatio'n (0.1). At the same time, from (2.3), and (2.4),
M(t) =~H(s, t) (~S
= 1 - Gz(t) s= l
+
f
t ~
0 k=O
kPk(u)H k 1 ( 1 , t - u )
CH(s,t_u)
~S
=1
dGl(u),
which is equation (0.2).
3. The Behavior q[" the Moments oJ Z ( t ) J o r Large t. There are standard results for the asymptotic behavior of the Bellman-Harris process. In this section we obtain the modification of these results to the system (0.1) and (0.2). As might be expected, the original result for the unmodified Bellman Harris process will be used in this analysis. Therefore, we now state that result; LEMMA 1.
I f M is the solution of the integral equation t
I
M(t)= 1 -G(t)+m
M(t-s)dG(s), 0
where m > 1, d G ( s ) > 0 , G ( 0 + ) = 0 , G ( , ~ ) = 1 and G is not a step Junction all of whose steps occur at regularly spaced points, then M ~ n 1 exp(at) as t--* ::~, where a satisfies the equation 1= m
f oc~exp( - as)dG(s) 0
and na=(m-1)/(am z
sexp(-as)dG(s)). 0
F o r the p r o o f of this lemma we refer the reader to Feller (1966), BharuchaReid (1960) or Athreya and Ney (1972). In this last reference the reader will also find a discussion of the various cases that can occur if m 2 =< 1. We note, however, that their discussion will be transferable directly to our (0.2) since this is exactly the Bellman Harris case. We can now state our result.
LAG
TIME
IN MICROBE
GROWTH
353
Letm 2 > 1 andsupposeGzsatisfiesalltheconditioninLemmal. Then the solution of(0.1) will obey,for large t THEOREM 1.
E(Z(t)) ~ n'1 exp(at), where a is the root o[the equation
fo
1 =m2
exp(-as)dG2(s)
while 77'1is given by,
j
" :7.~,
tTrl ~-111
t?l 1 (s)exp( -
as)dG1 (s).
0
Proof Let c > 0 be given and then choose Tso that 1 - G I ( T ) < e approximation M(t) ~ t71 exp(at) is valid for t > T. Then, since
and the
t
E(Z(t))=I-G1(t)+
I
ml(s)M(t-s)dGl(s), 0
we see that for t > T t-- T
E(Z(t))=I-Ga(t)+
I
ml(S)M(t-s)dGl(s)
do t
+
I
ml(s)M(t-s)dGl(S).
Now since t is large the first term will be less than e, while in the first integral we see that the argument of M is t - s > t - ( t - T ) = T a n d hence M may be replaced by the asymptotic form nl exp(at). Doing this we see that,
f
t- T 0
ml(s)M(t-s)dGlls)~
ft
7"
ml(s)nlexp(a(t-s))dGl(s)
0 t
7'
ml(s)exp(-as)dGllS).
= n 1 exp(at) (}
354
MARK
J. C H R I S T E N S E N
AND
R. S H O N K W I L E R
As for the second integral, since ml ( s ) < K and t > T we'see that
f
'
ml(S)M(t-s)dGl
(s) <2Kna
exp(aT)(Gl(t)-Gl(T))
t-T
< 2 K n l exp(aT)e. Thus we see that, for t > T,
E(Z(t)),,-nl exp(d! )
I
t - 7"
ml(s)exp(-as)dG1 (s) 0
+ 2Knl exp(aTle ml(s)exp(-as)dGl(s)+2Ke, exp(a(T-t)) .
= n l exp(at) 0
Therefore, as t--, oo, the expression in braces becomes,
f
,- r m,
(s ) e x p ( - a s )d G ~(s ) + 2Ke exp(a( T - t ) )
0
f
oo
---,
ml(s)exp(-as)dGa(s)+O. 0
Thus the proof is completed.
4. Extinction Probability. In this section the branching probabilities q,(t) and p,(t) (n=0, 1.... ) may vary with time. Let p denote the probability that starting from 1 newly born particle in regime 1 at time 0, extinction of the colony occurs at some time t < ~ . p = lira P(Z(t)=O). t -+ 30
I f a is the extinction probability for the process H(s, t) describing regime 2 (cf. Athreya and Ney, 1972, p. 143), then
P= ~ ~rk k=0
qk(u)dGl(U). 0
(4.1)
L A G T I M E IN M I C R O B E G R O W T H
355
This is easily seen, for denoting by I t the indicator function of the interval [0, t], we may write equation (2.5) as P(Z(t)=O)=
I,qk(u)Q~(t-u)dGl(u). k=O
0
The integrand of this expression as a function of t is bounded by the constant 1 which is a dG1 integrable function of u. Hence Lebesque's Dominated Convergence Theorem applies and since Q~ (t - u) ~ ak as t ~ ov we obtain (4.1). 5. Sample Growth Curves. In this section we present a number of the population growth curves that result from the system (0.1) and (0.2) for certain choices of the distribution functions G1 and Gz and of the branching probabilities qi and pj. (ij = 0, 1, 2 .... ). Our principal interest here is to study the lag or initial region of the growth curve. Since this region is primarily affected by the choice of parameters for regime 1, we shall assume, as is often done, that the process in regime 2 is Markovian and thereby obtain a simple closed form solution for the population growth function M ( t ) in regime 2. Specifically the life-length distribution function G2 must be exponentially distributed
Gz(t)=l-e
~a+u)t, t > 0 ,
if it is assumed that future states depend only on the present state of the process in regime 2 and if the branching probabilities # Po-2+it
and
2 P2-2+it
are independent of time. But for this exponential distribution the integral equation (0, 2) has the solution M (t) = e I~- ,~t. Hence the overall population size is given by t
E ( Z ( t ) ) = 1 - G1 (t) +
S
e ~ - ~'~"- U)ql (u)dG1 (u).
0
In regime 1 we take the life length to be gamma distributed dGl(t) dt
fl~+l F(a+l)
t ~ e -/~t
(5.1)
356
MARK J. CHRISTENSEN AND R. SHONKWILER
with parameters e and ft. This is a very flexible distribution over the range of parameters capable of including, for example the Markovian case by choosing =0.
Moreover we incorporate the seldomly used but critically important authenticity feature of allowing the branching probabilities in regime 1 to vary with time. Specifically we put ql(t)=ct~e
-bt
and
qo(t)=l-q~(t)
with a, b, and c independent parameters satisfying a>0,
b>0,
and c ( a / b )a e a = m a x ql ( . ) =< l.
This permits the simulation of an "adaptation window" centered around the maximum point for ql namely t = a / b . TABLE I
Graph number
Density c~ fl
a
Adaptation probability b
c
1
0.1
0.367
10
160
*
2
0
1
2
32
*
3
0
0.3
2
8
*
4
2
1.0
0
0.9
*
5
0.1
0.733
0.1
0.8
*
6
0
4
0.1
0.4
*
*c is chosen so that the maximum value of ql {I ) is 1.
Description G a m m a life length, wide plateau at T; thin window at T/48 Exponential life length averaging T/3: m e d i u m window at T/48 Exponential life length averaging T; m e d i u m window at T/12 G a m m a life length : m e d i u m at T; m e d i u m window at 0 G a m m a life length, wide plateau at T/2; wide window at T/24 Exponential life length averaging T/12; wide window at 7"/13
LAG TIME IN MICROBE GROWTH
357
The graphs which follow were obtained by a computer graphing program evaluating the equation (5.1) at 40 abscissa points. We used the same parameter values in regime 2 for each curve, namely 2=0.316
and
#=0.017,
which corresponds to an expected life-length of T--3 hours and a doubling probability of 95 per cent. This data is based upon the literature (cf. Shonkwiler, 1977). 3--
2:--
6
5
o 13.
I O
I
I
I
2
I
I
I
I
3
4
5
6
Time
Figure 1.
As our objective is to display the spectrum of behavior possible in the initial region of the growth curve, we choose a variety of values for the parameters c~,/~, a, b and c toward that end. The possibilities can be descriptively enumerated as follows. The distribution is taken either as exponentially distributed (~ =0) with expected life length times, 1//~, taken to be T, 1/3T, or 1/12T, or as gamma distributed with a medium (c~=2) or wide (~=0.1) plateau and expected life length times, (c~+ 1 )//~, taken to be 1/2 T or T. The adaptation probability q 1(t) is taken to have an adaptation window (maximum) at 1/12T, 1/2T, 1/48Tor 0 which is either wide (a = 0.1 if not at 0) medium (a = 2 if not at 0) or thin (a = 10 if not at 0). Table 1 is a list of the parameter values used for the six graphs. Inspection of these growth curves shows that this model is capable of describing a wide variety of behavior depending upon the profile and relative
358
MARK J. CHRISTENSEN AND R. SHONKW1LER
p l a c e m e n t of the a d a p t i v e w i n d o w . In p a r t i c u l a r c u r v e 5 s h o w s a p r o n o u n c e d d e l a y b e f o r e the o n s e t of e x p o n e n t i a l g r o w t h .
LITERATURE Athreya, K. B. and P. E. Ney. 1972. Branching Processes. New York: Springer-Verlag. Bellman, R. and T. Harris. 1952. "On Age-Dependent Binary Branching Processes." Ann. Math., 55, 280-295. Bharucha-Reid, A. T. 1960. Elements of the Theory oJ Markot: Processes and 7 heir Applications. New York: McGraw-Hill. Feller, W. 1966. An Introduction to Probability Theory and its Applications, Vol. II. New York: John Wiley & Sons, Inc. Mode, C. J. 1969. "'Lag Time in Cell Division from the Point of View of the Bellman-Harris Process." Math. Biosci., 5,341-345. Shonkwiler, R. 1977. "Consequences of the Growth and Division Model for Mitochondria Assuming No Cellular Control." Bull. Math. Biol., 39, 613-618. Waugh, W. A. O'N. 1955. "An Age-Dependent Birth and Death Process." Biometrica, 42,291 306. RECEIVED 11-22-76