SOME STGCHASUC MODELS FOR MECROGRGGNESGS DEATH mazes Hussis fer the Degree of Ph 0. Wm M STATE UNVERSITY PHZSZR CK MERGE GEYER 198? Ff? )G588773 '! M WE us I Jlll'l' mu m 1121: WWW ”WU! 2/ 3 93 00766 7243 This is to certify that the thesis entitled SOME STOCHASTIC MODELS FOR MICROORGANISM DEATH KINETICS presented by Frederick Pierce Geyer has been accepted towards fulfillment of the requirements for _Eh._D..__ degree in Agricultural Engineering @MQAQ Major professor Date_Ju.1_y_Z_._l.9.61_—_ 0-169 ABSTRACT SOME STOCHASTIC MODELS FOR MICROORGANISM DEATH KINETICS by Frederick Pierce Geyer Mathematical models currently used for microbial death kinetics are deterministic. Little attention has been given to population fluctuations arising from random aspects inherent in death processes. Modern probability theory was used to show how quantitative values can be assigned to the influence of these factors on the popu- lation size during a reduction process. Consequently, real death processes were not studied, but instead a method of mathematical modeling was derived and illustrated. Death processes were considered as Markov processes with a continuous time parameter. Chapman-Kolmogorov equations were derived and solved by the use of probability generating functions. Models for organisms with one and two viable states were considered for both time dependent and constant death rates. The stochastic models obtained gave theoretical probability distributions for the discrete levels of possible population size during a reduction process. The mean of the probability distribution derived was equivalent to the prediction of deterministic models. The latter was Frederick Pierce Geyer found to give a good approximation of the stochastic model for determining the lethal treatment required to sterilize a population of microorganisms. From the theoretical probability distributions de- rived, several methods were given to simulate a population reduction process. These techniques were computer pro- grammed to generate simulated death processes. Experimental evidence of the predicted probability distribution was considered for a homogeneous population with a constant death rate. The results were inconclusive because the predicted variation was usually smaller than the expected variation due to errors of measurement and observation. _For a model with two viable states, statistical esti- mation of transition parameters was considered. The least squares estimators required the simultaneously solution of a complex system of non-linear equations. The methods of successive substitutions and the generalized Newton method were developed. Their application was successful for data simulated according to the derived probability distributions. But these techniques did not give convergence for data with deviations larger than predicted by the stochastic models. Approved é ZZU %’ aJor Professor and Department Chairman Date 890*Z; ¢:/ . . Simulated Population (State 2) for Model with Two Viable States—-Method Using Equations 5.1A—l6. Simulated Population (State 2) for Model with Two Viable States--Method of Generating Time Intervals Between Each Transition . . . . . Simulated Death Process for Model with Two Viable States--Method Using Equations 5.14—16. . . . viii Page 31 37 39 H2 45 52‘ 54 67 75 79 82 Figure . Page 5.5. Simulated Death Process for Model with Two Viable States-—Method of Generating Time Intervals Between Each Transition . 83 6.1. Simulated Data for Test of Numerical Procedures . . . . . . . . . . 93 ix LIST OF APPENDICES Appendix A. Stationary Transition Probabilities B. Lagrange' 3 Method for Solution of a Partial Differential Equation of First Order . . . . C. Solution of Differential Equation for Stochastic Model with Intermediate State . . . . . . . . . Page 106 108 113 c, cl, etc. CV exp(a) g1 and gl(t) h h hl and hl(t) h2 and h2(t) NOMENCLATURE arbitrary constants. coefficient of variation event a e transition rate, per unit time, from state I to state 3 (at time t) death rate, per unit time death rate, per unit time transition rate, per unit time, from state 1 to state 2 (at time t) transition rate, per unit time, from state 2 to state 3 (at time t) integer number integer number integer number representing total number in states I and 2 initial number in state 1 number in state 1 at time t initial number in state 2 number in state 2 at time t initial size of population size of population at time t probability transition matrix probability xi pJ(t) p£,m(t) R and RN s(t) t u V X(t) A W 0, then the following two properties hold: Property 1 in(t’s) 1 0 (1.5) N Property 2 Z i=1 pji(t.8> = l (1.6) The following special case of the Chapman-Kolmogorov equations can then be obtained: "MZ pk1(u.8) ka(u,t) in(t,s) (1.7) J 1 0 i s < t < u Matrix notation may be used where P(t,s) is a N by N transition matrix. Equation 1.7 can be written P(u,s) = P(u,t)*P(t,s) (1.8) If the transition probabilities depend only on the difference t-s and not on the initial value, 8, the transition probabilities and matrix are stationary. Equation 1.7 can be written pki(8+t) = 321 ka(t) in(S) (1.9) s, t > 0 For matrix notation, this becomes P(s+t) = P(t) * P(s) (1.10) A detailed derivation is given in Appendix A. The fundamental stochastic differential equations can then be derived from Equation 1.7. If u is replaced by t +tAt and s by t in this equation, it can be differenti- ated (exact differential for a stationary model) with re- spect to t by taking the limit as At + 0 after a few algebraic manipulations. This yields the Kolmogorov "for-. ward system" of differential equations. 0n the other hand, if u is replaced by s and s-by s - As, the differential may. be obtained with respect to the variable 8 yielding the "backward system" of equations. Kolmogorov (1931) was the first to derive these two systems of equations. A more' detailed derivation may be found in most stochastic texts such as Doob (1953, p. 235) and Feller (1957, p. 423). Stochastic Population Models 3 The classical theories of population growth treat the size of the population as a continuous variable that proceeds deterministically throughout the whole process. The fundamental assumption is that the future development of the population can be exactly predicted once the state at some initial point is completely specified. Feller (1939) was the first to use the methods of Kolmogorov to treat population change as a Markov process in continuous time. His study led to a birth and death process with a discrete number of states. Much earlier, Yule (1924) had given a stochastic birth process in con- nection with the mathematical theory of evolution. He considered the creation of new species by mutation as a random event. Feller's work was further developed by Arley (1943). He used a simple birth and death process in the stochastic theory of the "cascade showers" initiated by cosmic ray. particles. In a series of papers, Kendall (1948a, 1948b, 1949) treated both homogeneous and non-homogeneous birth and death models. He also illustrated the use of gener- ating functions for these processes. Following Kendall, many other authors have contributed to this field. Among the special cases considered is the inclusion of immi- gration into the population. ~r 3See Lotka, 1945, for a review of deterministic models of population size. 10 As a partial development of the birth and death models, the simple stochastic death model can be easily obtained. While early population studies were interested in processes that involved growth or mutations, death events were only included to make the model more realistic. Consequently,the study of population reduction with stochastic death models did not-receive much consider- ation. However, the concept that human death is a random phenomenon has a long history. Medieval artists often pictured death as something that "sooner or later" en- slaved the individual. During the Great Plagues the mysterious ways in which death took the lives of many and left others untouched created a marked.impression that death was an unpredictable event. Similarly;the concept of chance was identified with death as that which obeys no rule and defies all measure and prediction. But with the development of the scientific theory of probability, chance took on meaning as a measurable quantity. Pearson (1897) argued that human death statistics could be identi- fied with probability distributions which could be defined in mathematical terms. From studies of human mortality statistics, he_found five periods of human life that showed regular chance distribution of mortality. In this section, individual deaths have been con- sidered in a population where life was normal. That is;~ some individuals may die, but others live and perform the ll usual functions of life. But in this study, the whole population is assumed under the stress of a lethal con— dition which inhibits the normal processes of growth and reproduction. Ideas about the action of radiation on organisms have stimulated a large number of deterministic models” for the survival of organisms. These can usually be classed into two general types, "hit" and "target" models. "Hit" theory defines an event (death) as taking place when the organisms have received a determinable number of "hits" or quantities of radiation. "Target" theory extends this concept by theorizing that there are two or more targets, Peach of which must receive one or more hits for an event to occur. Bharucha-Reid and Landan (1951) suggested a proba- bility model for radiation damage. They theorized a. chain of states with the ends being the absorbing states of death and complete recovery with immunity to further destruction. For hypothetical transition rates from one state to the next, the time dependent probabilities of reaching the two absorbing states are derived. Hoffman (1957) postulated that death of individual cells was a random event, but he did not give any mathe- matical models. Recently, Fredrickson (l966a)has “See Zimmer (1960) for an extensive review with references for death models in radiation biology. 12 suggested the use of stochastic models to describe the killing of microorganisms. He gives a probability model of organism viability for three different cases. The first model is the simple death process for a homogeneous population. Geyer (1966) also illustrates this model for the death of microorganisms. The second model gives the time-dependent probability that a clump of organisms has one or more viable organisms remaining. In the third model, Fredrickson derives the stochastic equations from a model suggested by Johnson (1963). This model assumes that the spore contains at least one each of several different types of subcellular structures. The proba- bility of a viable organism remaining is then obtained in terms of the destruction of all the different substructures. The basic stochastic death model for microorganisms was also developed by Terui (1966). He used this model to predict the most probable time to kill a population of microorganisms. Some additional details of the work of Terui and Fredrickson will be given in Chapter 3 where the basic stochastic model is analyzed. CHAPTER 2 ELEMENTARY STOCHASTIC MODEL Fundamental Observations Microorganism survival is currently considered to be reproducible according to deterministic laws. Fluctu- ation is ascribed, sometimes correctly and sometimes in- correctly, to experimental error. This study proposes that fluctuation in part is due to the random processes basic to the cause of death. And death kinetics for microorganisms are irreproducible processes. The exact cause of microorganism death1 is not known, but a number of rational explanations have been proposed for some lethal agents. If death is induced by moist heat, Rahn (1945, p. 39) has reasoned that death results from the denaturation of a single molecule. While other theories do not support this position, most agree that some unimolecular or complex molecular action occurs caus- ing the loss of reproduction ability. If this is the case, the hypothesis of Bartholomay (1957, 1958) may be utilized. He- argued that molecular processes are random processes. lDeath is defined in the usual sense of nonviability when placed in a favorable environment. 13 14 To substantiate this proposition, Bartholomay considered the random implications of modern chemical reaction theory. From this viewpoint, randomness may be found in the Brown- ian—like motions of molecules, in the random intermolecular collisions, and in the accompanying intramolecular "random walks" from one discrete quantum energy level to another. Bartholomay (1962) has extended this line of reasoning to enzyme kinetics and concluded that stochastic models should be used for enzyme reactions. This conclusion can be rele- vant to the theory of Isaacs (1935) that explains cell death from disinfectants as a result of enzyme inactivation. The action of chemical agents on microorganisms can also be considered as a molecular process. In some circum- stances, the individual chemical particles execute Brownian motion with small, rapid steps in a random manner to pene- trate and destroy the cell membrane. In the case of irradiation, death may be caused by X-rays, gamma-rays, and alpha-rays from radioactive material, fast electrons (cathode rays and beta rays), and other fast charged particles produced by the use of accelerators. In all these, energy is transferred in discrete quantities, and the time between emissions is a random variable. Rutherford gt_al. (1931) observed this randomness early in this century with radioactive material. The waiting times between decompositions are often described with a negative exponential distribution and the number of emissions by a. Poisson distribution. In-addition, the spatial distribution 15 of the radiation in the medium containing the cells can also be considered a random variable. The absorption of radiation by organisms has been theorized to have effects such as local or point energy release, molecular trans- formations following quantum Jumps, polarization, sepa- ration of charge and production of free organic radicals (Zimmer, 1960, p. 15). These actions collectively or singly are a consequence of statistical properties of the organisms and of the basic constitution of matter. Without further consideration of the causes of micro- organism death, the evidence from the several cases con- sidered indicate that one or more random factors contribute to all death processes. The modern theory of quantum mechanics establishes a comprehensive foundation that there is a basic physical randomness of molecular motion within all organized protoplasm. But an exact or deterministic relationship may appear on the macrOSCOpic scale. This gives an illusion that the process is reproducible. Usually, the instrumentation lacks the sensitivity to measure fluctuation on the microscopic scale. On the other hand, the growth of microorganism populations can be ob- served accurately. This is possible by measuring the time intervals between cell divisions with a microscope (Kelly, 1932). But since death is not an observable event such as cell division, better instrumentation is required to ap- praise the occurrence of individual deaths. l6 Derivation of Elementary Model This section describes the complete derivation of a non-homogeneous death process. It illustrates the mathe- matical techniques used in this study and presents the derivation for the time dependent case. Derivations for the homogeneous process may be found elsewhere. For example, Bailey (1964, p. 90) considers the homogeneous death process. A lethal environment will be assumed and not speci- fied as to whether it results from heat, chemical poisons, irradiation, etc. The lethal condition is applied at time zero with uniform intensity throughout the initial population of organisms. In order to derive the mathe- matical model for the population during the process, the following axioms are accepted. 1. A probability parameter h(t) is defined so that the probability of death for any organism dur- ing a short interval t, t+At is h(t+¢At). The function h(t) represents the death rate of the organism at time t. And ¢ is chosen so that: t+At g? f h(T)dT = h(t+¢At) (2-1) t 17 The death rate may be defined in terms of physical and chemical properties of the organism and the environment, but this re- lationship need not be established to derive the general stochastic model. 2. The probability of more than one death during the interval t, t+At is o(At) where o(At) is the zero order of At. That is, o(At) is some function of At such that the limit of 2%%El as At approaches zero is zero. 3. The Joint occurrence of events occurring in non-overlapping time intervals is statistically independent. Thus, the probability of two or more of these events is calculated by multiply— ing together the probabilities for each. Axiom 1 implies that the deaths of cells are inde— pendent events. Starting with an initial population No’ let n(t) be the random variable representing the number of viable cells at time t. Note that n(t) is a time de- pendent discrete random variable with a finite number of states. The probability of having n living organsism at time t is designated pn(t). This probability is not directly obtainable, but it can be derived from the stochastic differential-difference equations (the Kolmogorov equations). The "forward system" of these equations may be obtained from consideration of pn(t+At). This probability can be obtained by applying Equation 1.7 of Chapter 1. 18 First, consider all events leading to n organisms at time t+At. Three mutually exclusive events (El’ E2, E3) may produce this condition. They are described as follows: E : The Joint occurrence of events E and E . l 11 12 E11: NO - n + 1 deaths occur in time t(o,t). By definition this probability is pn+l(t). E12: One death occurs among the n + 1 organ- isms in time At. According to axiom 1, this probability is (n+l)h(t+¢At)At Therefore, according to axiom 3 the probability of E1 is given by the equation p{El} = p{Ell}p{El2} = pn+l(t) (n+1)h(t+¢At)At (2.2) E2: The Joint occurrence of events E21 and E22. E21: NO - n + i (i l 2) deaths occur during time t. This probability is pn+i(t)° E22: During t, i (i 1 2) deaths occur among p{E2} = p{E21}p{E22} = i the n + i organisms. According to axiom 2, this probability is o(At). Thus, 1 2 pn+i(t) 0(At) (2.3) 19 E : The Joint occurrence of events E and E 31 32 E31: NO - n deaths take place during time t. This probability is pn(t). E32: No deaths occur for the n remaining organisms during At. Since the proba- bility of one or more deaths is nh(t+¢At)At + o(At), the probability of no deaths is 1 - nh(t+At)At - o(At). According to axiom 3, p{E3} = p{E3l}p{E32} -- pn(t)|_l - nh(t+¢At)At - o(At)] (2.14) Since events El’ E2 and E3 are mutually exclusive ways in which pn(t+At) may occur, the probability for each of the three events is summed, pn(t+At) = p{El} + p{E2} + p{E3} (2.5) Substituting the expressions obtained for the terms on the right side of this equation, pn(t+At) = pn+l(t)(n+l)h(t+¢At) + 1:2 pn+i(t)o(At) + pn(t) 1 — nh(t+¢At)At - o(At)] (2.6) By subtracting pn(t) from both sides of this equation and dividing by At, it becomes: 2O Pn(t+At) - pn(t) At = pn+l(t)(n+l)h(t+¢At) o(At) pn+i(t) At + Z i 2 - pn(t) n h(t+¢At) (2.7) Now consider the limit as At tends to zero, the left side d (t of 2.7 becomes the derivative pdt ). By axiom 2, the value of 9L%%1 goes to zero, and h(t+¢At) becomes h(t). Consequently, Equation 2.7 can be written d pn(t) “_‘6577 = pn+l(t)(n + l)h(t) - pn(t) n h(t) (2.8) where n = 0,1,...,N 0 Equation 2.8 represents a system of NO + 1 equations. Each will have two terms on the right side of 2.8 except for the case where n is 0 or No’ Since the state NO + 1 can not exist, pNO+1(t) has value zero. Therefore Equation 2.8 is reduced to the following equation for n = No' d p (t) ———§9——— = (t) N h(t) (2 9) dt pNO o ' And if n = 0, Equation 2.8 yields d po(t) “EF_—’ = pl(t) h(t) (2.10) 21 The whole system of equations represented by 2.8 can be efficiently represented by the matrix notation: dP(t) _ ‘51?“ — h(t) A P(t) <2-ll> where P(t) is the vector pN (t), pN _1(t), pl(t), pO(t)] o o ' and A is the No+l by No+l matrix' P .No 1 NO -(NO-l) (No-l) —(NO—2) Solution of the System of Differential Equations The system of differential equations derived in the previous section may be solved in several ways. Provided h(t) is defined, these equations can be integrated suc- cessively starting with n = No' This process will yield the general solution: ’6 t pn n < > < E x —————- = E n+1 x h t p t) n=0 dt n=0 n+1 N O n - 20 n x h(t)pn(t) (2.17) n: This equation can be rewritten in a more suitable form by shift of axis to discard meaningless terms. No n dpn(t) No n-l Z x dt 2 n x h(t)pn(t) n=0 n=0 N O n-l - x z nx h(t)pn(t) (2.18) 24 Equations 2.15 and 2.16-can then be substituted into 2.18 and a partial differential equation involving the generat-. ing function is obtained. 2.91%.21 = h(t) (1 — NEE—@5313)— (2.19) Using Lagrange's method of auxiliary equations as described in Appendix B, the following ordinary differen- tials have the same solution as 2.19. dQ(x t) _ dx 3 dt *‘L-o ‘ h(tm - x) "If (2'20) Two independent solutions of 2.20 are Q(x,t) = c (2.21) l, and (x - 1) exp(- gt h(r)dr) = c (2.22) 2 with c1 and c2 arbitrary constants. Therefore,0(x,t) is t some function of (x - l)exp(- f h(r)dr). For the initial 0 conditions, t(x,o) = xNO, the general solution is t t No Q(x,t) = x exp(- f h(r)dr) + 1 - exp(- I h(r)dr) o o (2.23) 25 By a series expansion of 2.23, the coefficients of xn are obtained. From the definition of the probability generating function, these coefficients are the values of pn(t). Thus, NO . n N-n pn(t) = s(t) [l - s(t)] (2.24) n where s(t) = exp(- fth(r)dr) (2.25) O and NO — NO! - (NO-n)! n! This result could be obtained directly from 2.23 by observing that it is a generating function for a binomial distribution with parameter s(t). Since s(t) is the sur- vival probability for any organism in the population and each organism was assumed independent, the binomial distri- bution of Equation 2.24-can be obtained from Equation 2.9 for a No of one. However, the method of solution given illustrates the techniques used for more complex models (Chapter 4) that can not be solved by simple methods. The distribution of the random variable n(t) can be used to determine common statistics such as the mean H and variance 02. These two statistics can be obtained directly from the probability generating function by applying the following formulas. 26 n(t) = w'(l,t) = NO exp(— ft h(T)dT) (2.26) O 2 _ 2 o (t) - w'Nl,t) + w'(1,t) - Ew'(l.t)l = NO exp(— gt h(T)dT)[l - exp(- gt h(r)dr)] (2.27) Next, the distribution function of the arrival time of an event may be obtained. Starting with NO organisms, the probability that at time t no event has occurred is given by 2.24 for n = No' Accordingly, t pN (t) = exp(- N f h(T)dT) (2.28) o O 0 This is also the probability that the first event happens at some instant greater than t. Therefore, the distribution function of the arrival time u of the first event is given by u F(u) = l - exp(- NO 5 h(T)dT) (2.29) and the corresponding density function is f(u) = F'(u) = Noh(u)exp(- NO xu h(t)dt) (2.30) 0 where h(T) i 0 for t(0,w) and h(t) = 0 for only a finite interval. This equation can be used to give the distri- bution of the time interval between any two successive events if the value of No is adJusted to the level of the 27 population and the time scale and h(T) are shifted to the point the last event occurred. The deterministic model of this process may be com- pared to the stochastic in several ways. Deterministic kinetics are based on the Law of Mass Action. Accordingly, any change in population is proportional to the size of the population. The ordinary differential equation —— = -n h(t) (2.31) describes the model. In this case, the population is treated as a continuous function of time,although n is discontinuous for a real death process. Equation 2.31 may be integrated to give t . n(t) = Noexp(- é h(T)dT) (2.32) This equation has the same form as the equation for the mean of the stochastic model (2.26),therefore n(t) and h(t) are the deterministic equivalent of fi(t) and h(t) in the stochastic model. Thus, the stochastic model is "consistent in the mean" with the deterministic model. Consequently, the deterministic model may be considered a special case of the stochastic model. The stochastic model would yield the same result as the deterministic only if a large number of cases were averaged. 28 The stochastic model not only predicts a fluctuation from the mean, but it specifies the expected size of these deviations. With this model, the reproducibility of the process can only be considered as a Joint probability distribution of two or more independent events whose proba- bilities are specified by Equation 2.24. However, it is possible to specify a range of values within which the process would be expected to lie for any level of signifi- cance desired. Because the probability distribution derived in Equation 2.24 is binomial, it may be approximated with a normal distribution according to the de Moivre-Laplace Limit theorem (Feller, 1957, Chapter 7). The error of this approximation will be small if the variance is large. If the variance is small, a Poisson distribution could be used as an approximation of the distribution because the variance will be small in the same intervals where the probability is very small or close to one. For a normal approximation, about 68 percent, 95.5 percent and 99.7 per- cent of the distribution would be expected to fall within one, two, and three standard deviations, respectively of the mean. CHAPTER 3 HOMOGENEOUS CASE OF ELEMENTARY MODEL Description If the organism death rate h(t) is independent of time t, the process is considered homogeneous and h(t) = h. Assuming a constant lethal environment, several micro- organism death processes exhibit homogeneous character- istics. Rahn (1945) and Stumbo (1965) both concluded that the death rate of spores by constant temperature heat in- activation is independent of time. Their conclusion was based on their own laboratory studies as well as those by other researchers. On the contrary, considerable evidence has been ob- tained that heat inactivation of spores is time dependent for some conditions. For example, Frank (1957) and Humphrey (1961) have documented time dependent cases. For spore irradiation death processes, neither the homogeneous nor the non-homogeneous elementary model seems appropriate. Instead the process is considered more complex, and a "target" or "hit" model (see p. 11,Chapter l) is usually given. 29 30 The homogeneous case of heat inactivation will be non-homogeneous if the temperature is not held constant. Under these conditions, the death rate is a function of temperature. Assuming that a death process is homogeneous, its probability distribution may be easily obtained from the non-homogeneous stochastic model in Chapter 2 (Equation 2.22). Consequently, N N -n pn(t) = [ 0] exp(-nht)[l — exp(-ht)] 0 (3.1) n L H(t) = NO exp(-ht) (3.2) 02(t) = NO exp(-ht)[1 — exp(-ht)] (3.3) Bailey (1964, p. 91), and Frederickson (l966a) gave this distribution for a stochastic model of microorganism death kinetics. As an example of this type of process, Figure 3.1 shows a hypothetical death process for the distribution of Equation 3.1. As is commonly done, the log of the number of survivors is plotted with a linear time scale. Thus, the mean is a straight line. To ex- tend this consideration to some of the unique features of the homogeneous model and its application, the follow- ing was developed. Number of Organism 31 K35 C —-Theoroctical Moan : ---Tnoonctical Moan - (1’3 Standard Deviations) P —-Slmulatod Response 4 Initial Population 40’ K) E- D. Ahflb -l . 0th Rate - 2° (tint-tits) p r I03 \ E" \ I \ \ h \ \ t \ \ \ 1’ \ \ \ \ \ K32g- \\\\\ : \ \ P \ \ I- \ \ t \ \ r \- \ \ )- \ \ I \\ \\ IO ;- \ \\ {I \ \ r \ \ b \ \ \ \ L \ \ \ \ "30 l 1 1 \ O 20 4O 60 so IOO IZO Ti me ' Figure 3.1.-—Survival Curve for An Elementary Death Process. 32 First, consider the variance of the process. The lines in Figure 3.1 showing the mean plus and minus three standard deviations may be deceptive. The time of maxi- mum variance or standard deviation is not evident from this figure. Since the variance (defined in 3.3) is a continuous function with value zero at the ends of the time scale (0,m), a time tm of maximum variance may be found for 0 < tm < w. The derivative of 3.3 with respect to time and solved for t in the usual manner to find a maximum yields: t = in 2 (3 4) The size of the maximum variance °m2 is then: No 0 = 3— (3-5) 2 At time tm, the mean is 1? or the expected popu- lation size is half that of the original population. Also, the time of maximum variance is solely a function of the death rate constant and the size of the maximum variance is only a function of initial population size. Considering another statistic of the process, the coefficient of variation, CV, is a continuously increas- ing function of time since CVpouQH oEHB mo coapmnococ mo mmooonm Sodom copmasefimll.m.m opswfim .Ama.m coapmsvmv spoon comm coozuom osae om o: om om OH .r .121 Jlrrr/JH/ com m C: OE III N m .m .mm HuA pa Hov OH as a m no a com u soapsasooa HsHoHaH uorietndog u: qgeq JeqmnN 40 used. Since each experimental sample is independent of the others and provides one data point, each simulated data point may be simulated independently. The problem is then reduced to simulating only one point for any time given the initial concentration and death rate. This point will have the binomial distribution given by Equation 3.1. The generation of random numbers with a known distribution is not difficult. Shreider (1964, p. 252) gives the procedure for continuous density functions. But this technique can be extended to the discrete case in the following manner. In order to obtain a number belonging to a set of random numbers, n1, having the probability mass function pn, generate a random number R with uniform distribution (0,1). Then choose the smallest ni such that pJ > R (3.17) This method may be simplified. Since pJ is bi- nomial, it may be approximated by a normal distribution. To generate a normal distribution, again assume that random numbers with uniform distribution (0,1) may be acquired without difficulty. Then a random variable, V, with a normal distribution (0,1) may be easily generated by the following equation derived by Box and Muller (1958). 41 1/2 V = (-2 9.0ge R1) cos(2nR2) (3.18) where R1 and R2 are random numbers with uniform distri— bution (0,1). The cosine function in the above equation may be interchanged with the sine function without chang- ing the distribution of V. Since V is normal (0,1), the required random number, nr, may be determined using the values of the mean and variance given in Equations 3.2 and 3.3. Thus, “s(t) n(t) + V*o(t) (3.19) An example of data generated using this procedure is given in Figure 3.4. A Comparison of Deterministic and Stochastic Models The deterministic model has at least two important faults. First, it assumes the population size is a real— valued continuous function of time rather than an integer- valued function of time. Second, it assumes that the population size at any given time will always be the same if the initial conditions are not changed. The second fault is the more obJectionable of the two because it ignores intrinsic random factors that may influence the destruction of microorganisms. The stochastic approach questions the assumptions of the deterministic model and thus its validity. By Number Left In Population 42 -—Process Mean x Simulated Data Point initial Population - IO£5 Death Rate -£32-102(flm 00m)" 5 N I "UVVV" IO' U V rT'UUI T “)0 1 1 1. n X 11 O 20 4O 60 80 IOO IZO Time Figure 3.4.-~Simulated Experimental Data For Elementary Death Process. 43 substituting probability relationships for deterministic ones, random fluctuation is expected even in the total absence of experimental irregularities. However, the amount of inherent variability may be small and unmeasur— able by experimental procedures. As shown in Figure 3.1 the expected range of most of the predicted variability is too small to be detected for the reduction of the first half of the population. On the other hand, the predicted fluctuations become fairly large relative to the mean as the pepulation grows small. The difference between the stochastic and determi- nistic equations may be compared for prediction of process times required for sterility of all organisms in a popu- lation. To obtain complete sterility,the viable popu- lation must be reduced to zero. For the deterministic model, the population only reaches zero as time approaches infinity. However, this model may be used to predict practical process times if the value of the mean is taken to represent the probability of viability for the whole population. To use this approximation,the mean must have a value less than one. Then the probability of any viable organisms remaining at time, t, is designated q in the following equation. q = NO eXp(-ht) (3.20) where q < 1 44 Solving this equation for t, the process time to obtain a probability of viable population, q, is Process Time (Deterministic Model) = - % 2n[§L] (3.21) o For the stochastic model, q is 1 - po(t) by definition. Using the value of po(t) given in Equation 3.1, the process time required may be specified as follows. N -1 Process Time (Stochastic Model) = - %-£n 1 — (l-q) O (3.22) To compare these two predictions, Figure 3.5 shows a plot of both for an initial population of 105. By using semi-logarithmic scales, the deterministic prediction is a straight line. This line does not take on any values for the initial time interval because the mean of the pro- cess must be reduced to one before the process time can be determined (Equation 3.21). This line can be considered 1 On the as a continuation of the mean line in Figure 3.1. other hand, the stochastic model gives a probability of via- bility for all points in time. It yields a line that asymptotically approaches the deterministic curve as the process transpires. By the time the probability of any 1The time scale in Figure 3.5 has been changed from that of Figure 3.1, but the initial population is the same. 0L IO _ ————————— \ I. \ t \ I \ b . P so" - .. 2 E . .3 __ 0 .. 5 .‘L , n .9. 2» '8; "3'2:- x P 5 : Initial Population - 105 § )- —- Deterministic Prediction :8 ‘ —-- Stochastic Prediction Q, )- 10-3 _ 1 1 1 1 1 1 4 1 1 1 02468l0l2l4I6I820 Process Development (Time -x- Death Rate) Figure 3.5.--Comparison of Deterministic and Stochastic Predictions of Viable Organisms Left in Population. 46 viability is down to .01, the difference between the two models is indistinguishable. Thus Equations 3.21 and 3.22. yield the same result for small q. This conclusion can also be procured by a Taylor series expansion of these two equations about the point q = 0. Since most processes designed to produce sterility would demand a probability of sterility of less than .01, the deterministic model is as good as the stochastic model for predicting the re- quired process time. This last conclusion implies that the deterministic model is sufficient for sterility appli- cations even though the stochastic model more accurately represents real homogeneous death processes. Experimental Evidence To find experimental evidence that stochastic models represent a death process better than deterministic models, the experimental variation from the mean was studied. Assuming a homogeneous process, the variation of the re- sult from the mean can be a result of two factors. First, the intrinsic randomness of the process will cause vari— ation. Second, errors of experimentation will also contri- bute to deviations from the expected population size. The latter factor can be caused by a large number of factors since microorganisms respond to many environmental factors. Hopefully, the researcher is able to control most of these variables. But this is a very difficult task. In all cases the part of the deviation from the mean caused by 47 experimental error and that part due to the true process variation can not be distinguished. As-already stated, the direct observation of mi- crobial death events is not possible with today's techno— logy. Therefore, death processes of large populations must be studied. Data of this type can be obtained by subJecting a large number of samples of an organism to a lethal condition and withdrawing samples at different time intervals and counting the number of organisms re- maining. Experimental data obtained by Deweyl working with Serratia marcescens were analyzed. He irradiated (X-rays) cells of this organism in an oxygen atmosphere. The re- sults of these tests are shown in Table 3.1. Ordinary microbiological techniques were used in counting popu- lations. This included the diluting of large populations to obtain a population small enough to count. First, the death constant, h, was determined from a linear regression of the logarithm of the percent popu- lation reduction versus the dose of irradiation received for each trial. A death constant of .585 (kilorads)-l was obtained. Using this value for h to determine the theoretical mean and standard deviation (Equations 3.2 and 3.3), the values in columns 5 and 7 of Table 3.1 were 1 See Dewey (1963) for a report of his experiments. The data reported here is not published in his article,but it was obtained through personal communication. 48~ mmma.m aa.m mmmm.mm .am.mm amm.ma ooo.oa N omaa.mm am.mom omma.omam sm.mmom moa.m ooo.am a mmmm.am ,mm.moa mmma.amaa mm.mmma oam.o ooo.am x maom.ma om.oma ambfi.mmm ao.mam mam.m ooo.am 3 mama.m mm.oa Haom.ma mm.mma oma.ma ooo.woa > msam.mm mm.ammm mmoa.mbmm mm.mmm.ma mam.m mmm.mo a maao.am ma.moo mmfla.moaa _mm.mHom oma.a mmm.mm a mmaa.aa mm.aw ammm.oam ao.amm mma.m mmm.mw m mmmm.a mm.mm aoao.mo o.om omm.ma moo.mma m amma.mm .am.ommm ammm.mamm mm.mom.aa Ham.m www.mm a mamm.om Hm.mmm mama.aama ao.mamm mam.a www.mo a mmoH.:H .oo.HHH oamm.aom mm.mam mam.m www.ma o mmoa.a mm.ma oaae.mm as.aa amo.ma mmm.mma z a:mm.am am.ammm oomm.mamm am.mam oam.m mmm.mm z mmma.am aa.aam moam.mma o.oam omm.o oma.mm a moom.m _aa.mm ammm.om “mm.ama oma.m mmmsmm a moms.m ma.HH ammm.ma mm.mm oaw.ma mmm.mm a ammm.mm mo.ammm omam.moa.aa .mm.mmm.ma oaa.m www.mm H mmem.wm ma.am: maam.mmaa o.omma o:m.: www.mm m mmfio.aa wa.am omam.amfl mm.mHH oaa.oa mmm.om a mama.a mo.mm oaoa.am _am.afl omm.ma mma.maa a a:am.a Ho.m omoa.m m.m amm.am ooo.oom.mm m mmaa.m mH.m mmmfl.mm o.mm mmm.am ooo.ooo.mmm a wmam.m oa.m mome.a mm.a Hma.mm ooo.ooo.omm.a o aamm. om. memo. o.H aao.mm oam.mam.mmm m meow. Ha. ease. mm. aoa.om oaa.oao.ooo a soapmfl>oa new: Scam and: .COfipmHowan mnemoaam coaumHSQom pudendum coaumfl>om houm< whom mmnlx HmeHcH mxm Hmoaponoose COHpmHsdom Hmoauonoone HGOfipnasaom .. .Ammmav meson no when HmpcmsHpoaxo mcoomoohme mfiumnpom mo Ammmpuxv,QOHumHemaaH no mammamc<11.a.m mqm<9 49 .coxmp mes moHdEmm oops» mo owmno>m cm omdmoon czonm one methane HMCOHpomhm H mmm:.:m nmmo.mm mmma.ma mmmm.m ma:>.mw mmmm.mm mmwo.ma omom.m mooo.am mmmm.mm mamm.ma mnea.m m:am.aw ommm.mm Hmmm.m mamm.m Nmmm.m~ mmmm.om :MHH.NH mmza.z maam.m> mmmm.mm mzmo.HH mmmm.m mmma.mm wooa.~m mm:m.=a ma.amam ma.oam mm.mH _-.:H em.amm: mm.ama mo.woa o>.mm em.m:: mm.ao: H0.NHH mw.mH :m.mmoa om.=~ am.mm ww.m mm.m:zm mm.mmm 0:.mma 3:.mm mw.amam ma.m:u Hm.mHH mm.mm .mm.mmzm 03.0mm NH.>H mama.amma oamm.mmofl ammm.aaa woma.em mmma.momm maao.HeHH Hmfle.mmfi mmoa.mm amam.maaa aomm.aaaa mmam.aaa aamm.mm mHmH.Homm aoaa.mma aamm.am mmma.aa mamm.mamm comm.aom moam.a:a ommm.mm moma.ammm ommw.aam aama.mma Homa.am mamo.oaam mmam.moaa mama.mam .Nm.mwwm o.omaa o.HmH o.mH o.oom.HH 0.0HOH as.mm 0.0H o.ooma 0.0a» 50.:0 o.e am.oamm em.mom o.mm pa.m o.mmm o.omma .wm.mmm o.mm um.ommm o.omma um.mmm .mm.:m o.mmo.aa o.omoa o.mmH soa.m mam.a mmm.m mam.H mmo.m moo.m moo.m mma.ma mmm.m ama.a ama.m mam.ma mmm.m OHH.A maa.oa omm.aa aom.m :Hz.m Hma.m www.ma amm.m :He.m Hao.oa mma.ma aom.m :H:.m amm.m ama.oa omm.mm aaa.oa omm.mm mma.oa mma.oa mme.oa mm>.oz moa.am mma.am oma.am ooH.Hm oom.oa oom.ma oom.ba oom.aa oom.oa momma: com on common mma.aa mma.aa mma.aa mom mm oom.mm mmo.om moo.om ¢ D< B< m4 m< a< m< O4 z< 2¢ Q< m< h< H¢ m¢ 0< m< m< Q< o< m¢ << 50. calculated. The difference between the theoretical mean and the population count is given in column 6. If the deviation from the mean was only a result of the stochastic variation predicted, its expected value would be the theoretical standard deviation. A chi— square test could then be used to test the deviations from the mean against their expected values. This test was tried, but most of the experimental deviations were so much larger (as much as 50 times larger) than the theo- retical ones that the test yielded a very negative result. However, this does not mean that the stochastic model was invalid because experimental errors could cause the de- viations to be much larger than the theoretical prediction. Thus, the experimental errors would dwarf the effect of the intrinsic variation on the outcome. It is also possible that the intrinsic variation causes a larger deviation than predicted by the model. By considering only the experiments for which the initial population was over one million, a good correlation was found between the theoretical and experimental devi- ations. Experiments A, B, C, D and E all had original populations over one million. They also had a population count of less than ten organisms after the irradiation treatment. A chi—square test for these five experiments showed the validity of the theoretical variance at a .05 level of significance. 51 An attempt was then made to correlate the experi- mental deviations to the theoretical standard deviations by a linear regression. A plot of these data is shown in Figure 3.6. From a least squares analysis of these data the following equation was obtained: Experimental Deviation = 10.87 Theoretical Deviation - 8.9 (3.23) Thus,the experimental deviations were about eleven times that expected. These two variables had a correlation co- efficient of .757. However, this correlation may be due to another variable such as the number of organisms re- maining after the treatment. The dilution of the popu- lation required to obtain a countable number tends to amplify the errors for large populations. Also, the theoretical variance decreases as population decreases for the range of the experiments since all treatments extended beyond the dosage of maximum variance (Equation 3.4). Therefore, both may be correlated with pOpulation size. In conclusion, the experimental results give evi— dence of the applicability of the theoretical model, but more accurate measurement techniques are needed to test the model. Experimental Deviation From Theorectical Mean. Y "MOO 900I 80ml 'WOO 60M) 50C) 4(K) 20C) WOO 52 .714J1121L.llllljllll'1 IO K) 20) 30 (40) 50) 60> 1“) BO) 90 KM) Theorectical Standard Deviation. X Figure 3.6.--Influence of Theoretical Standard Deviation on Experimental Deviation from Process Mean. CHAPTER 4 A DEATH PROCESS WITH AN INTERMEDIATE STATE In the elementary model of Chapter 2, there was only one possible transition: a change from a viable state to a non-viable state. As an extension of this model, one can theorize that the organism may exist in two different viable states. One state may be more death resistant than the other or one state may be a result of the lethal environment. The second viable state could also be a result of an adjustment to the lethal environ- ment or a partially destroyed state. To construct a stochastic model for this case, let the two viable states be denoted by 1 and 2, and the non— viable state as state 3. Assume organisms in states 1 and 2 both may become non—viable by transitions to state 3. In addition, assume organisms in state 1 may make a transition to state 2. All other possible transitions between the three states is assumed non-existent. Visually, this model could be represented by Figure 4.1. 53~ 54 State State A Figure 4.1.-—Three State Model with Transitions Shown by Arrows. To derive.a stochastic model for the semi-closed time interval 0 :_t < w, let the number of organisms in each state be represented by the discrete random variables as follows: m(t) = number in state 1 at time t. m(t) = number in state 2 at time t. n(t) = number in state 3 at time t. The Joint probability distribution pl gt; is defined as .9 9 the probability of having 2 organisms in state 1, m organisms in state 2 and n organisms in state 3 at time t. It follows that: p (t) = 1 > (4.1) for all t (0,w). 55 Further, assume that the initial concentration in each state is known; that is, 2(0) = Lo’ m(o) = M0 and n(o) . No' Under these conditions, n(t) = LO + MO + NO — s(t) - m(t) (4.2) (t) p£,m,n (t) bility p£,m since n(t) is known if 2 and m are chosen. Therefore, may be reduced to the two state proba- Equation 4.1 then reduces to: z p, (t) = 1 (11.3) The following system of 6 axioms is accepted to derive an analytical expression for pg’ét). l. The probability of a transition of any organ- ism in state 1 to state 2 during a short time interval t, t+At is represented by the transition parameter hl(t+¢lAt)At. The function i1 is chosen such that t+At 1 - Kt é hl(T)dT - hl(t+¢1At). (4.4) 2. The probability of a transition of any organism in state 2 to state 3 during a short interval t, t+At is h2(t+¢2At)At. And ¢2 is defined such that d'H :3' f% (T)dt = h2(t+¢2At). (u.5) 56. 3. The probability of a transition of any organism in state 1 to state 3 during a short interval t, t+At is given by the transition parameter gl(t+¢3At)At, where o3 is defined in the same way as ii and o2 in Equations 4.4 and 4.5. 4. The probability of more than one transition among the three possible types in time interval t, t+At is o(At). Where o(At) is the zero order of At. That is,o(At) is defined such that 2:20 W: 0. (4.6) 5. All possible transitions other than given in axioms 1, 2, 3 and 4 have probability zero. 6. The Joint occurrence of events occurring in non-overlapping time intervals is statistically independent. (t+At) m , consider the 9 To obtain an expression for p2 various transitions which, starting at time t, can lead to values A and m for states 1 and 2 respectively at time t+At. Since the axioms listed above dictate that the random variables can either increase by one, decrease by one or remain the same without having a probability of order o(At), the transitions shown below in Table 4.1 must be considered. 57 TABLE 4.1.--Transitions for model. State At Transition State At Event Time t During At t+At El 2+1, m-l From state 1 to state 2 £,m E2 2, m+l From state 2 to state 3 2,m E3 2+1, m From state 1 to state 3 2,m E4 1, m No transitions z,m The probability of the first three events listed in Table 4.1 above is the product of the probability of two independent events. The first event is the occurrence of being in the state specified at time t. Its probability is given by definition. The second independent event has the probability of a transition during time At. Axioms l, 2 and 3 specify these probabilities. Consequently, (t) p{El} = p£+l,m-l (2+1) hl(t+¢lAt)At (4.7) p{E2} = pz’éii (n+1) h2(t+¢2At)At (4.8) piE3} = p1+1,m (2+1) sl(t+¢3At)At. (4.9) If more than one event of any of the three types given above or any combination of them occurs, then o(At) must be a factor in the probability according to axiom 4. 58 Let E4 be the event that the process is in state £,m after two or more transitions during time interval At. If r and s are any combination such that after two or more transitions state r,s becomes state 2,m, p{Eu} = 2 pr,s(t) o(At) (4.10) r,s Event 5 can now be defined as occurring if no transitions take place during interval At. Because all possible events must have a total probability of one, the probability of event 5 is one less all probabilities of one or more transitions. Therefore, P{Es} = P£,m(t) (l - q) (4.11) where q = £hl(t+¢lAt)At + mh2(t+¢2At)At + 1gl(t+At)At + o(At). (4.12) Since the five events are mutually exclusive ways in which pn(t+At) may occur, the desired probability is the summation of the probabilities for each of the five events. P£,m(t+At) = p{El} + p{E2} + p{E3} + p{Eu} + p{E5} (4.13) 59 (t+At) (t) (£+l)hl(t+¢1At)At p1,m a p2+1,m-1 (t) p£,m+l + (m+l)h2(t+¢2At)At + p£+§E% (2+l)gl(t+¢3At)At (t) + p9"m 1 — £hl(t+¢lAt)At - mh2(t+¢2At)At - 2gl(t+¢3At)At - o(At)] + z pr,s(t)o(0t) r,s (4.14) This equation can then be reduced to a differential equation by-subtracting p2 m(t) from both sides, dividing by At and 9 taking limit as At approaches zero. This equation is: d p£,m(t) (t) a. = <2+1>h1 + p.131 new.) + “Afr; (2+l)gl(t) - p£,m(t)[:£[hl(t) +gl(t)] + mh2(t)] (4.15) Taking into account that p2 m(t) has value zero for D 2 > L0’ m > MO and 2 or m less than zero, Equation 4.15 60 may be reduced to a simpler form for the end points of the process. d PL M (t) o? o (4.16) d p0,g(t) dt = pl,0(t)gl(t) + p0,lh2(t) (4.17) As indicated in Chapter 2, the easiest method of solving the system of equations in 4.15 is by using a generating function. Therefore, the following bivariate probability generating function is defined. L M o o 2 m w(x,y,t) = Z 2 p2 m(t) x y (4.18) i=0 m=0 ’ From this definition, ' (t) at 2 m d p2 m —- = Z x y ’ (4.19) at 2,m dt 33 = 2-1 m 3x 22 2x y p£,m(t) (4.20) ,m ii a 2 m-l 8y 2 x y £,m(t) (4.21) 61 Multiplying each term in Equation 4.15 by xgym and summing over all possible values of 2 and m, Equation 4.15 becomes d p (t) 1 m 2 m 2 m E x y _____.L____ = 2 (£+l)x y h (t) p 2,m dt 1,m 1 £+1,m-l 2m + 2 (m+l)x y h (t) p 2,m 2 £,m+l 2 m + 2 (1+l)x y s (t) p - z lxzym hl(t) + 31(t)} p£,m(t) £,m - z mxzymh2(t) p2 (t) (4.22) 2 m ,m 9 The above can be reduced to a single partial differ- ential equation by substituting the functions defined in Equations 4.18 through 4.21 along with a few shifts of axes . 31(X.y,t) , a a a at yhl(t) 3% + h2(t) 5% + 31(t) 3% - x[hl(t) + g1(t)] %% - yh2(t) %% (4.23) 62 This may be reduced to 0= h(t)+ (t)-xh(t)+ (t) 9i y 1 gl 1 211 3x + h2(t)(l - y) 3% — 3—5— (4.24) Using LaGrange's method of auxiliary ordinary differential equations (Appendix B), Equation 4.24 has the same solution as: SE dx 0 yhl(t) + gl(t) - [hl(t) + gl(t)]x _ d1 2.0113. ' h2(t)(1 - §)* -1 (4.25) The solution for this was obtained for initial conditions that there were L0 in state 1 and Mo in state 2 at the beginning of the process. The complete derivation is given in Appendix C. The solution is: L M t(x,y,t) = [X8 + y-y + l - B - y] 0 [ya + l - {I O (4.46) where t a = exp(- f h2(T)dT) (4.27) O B = exp[- at [I'll-(T) + gl(T):|dT] (4.28) _ t t y - a g hl(1)exp[- é [Ll(r) + gl(t) - h2(T)]dT]dT (4.29) 63 The Joint probability distribution for t in state 1 and m in state 2 is found by selecting the coefficients of x2 and ym in Equation 4.26. By series expansion, L m L -2 L —1-1 M pg m(t) = O 8‘ 2 0 vi l-B-Y O ’ 1 i=0 i M -m+i O am-i(l-a) O m—i (4.30) The summation in the above equation includes all L -2 terms from zero to m, but only the terms where O M i and O are defined need be included. All others have m-i value zero. If MO is zero, L L -2 L -£-m o 2 o m o p,,m(t) = 2 e m y 1-e-y (4.31) If L0 is zero, the Joint distribution is not required because i will start at zero and remain there for the whole process. In this case the stochastic model reduces to the elementary model presented in Chapter 2. According to Equation 4.30, the probability of extinction is given by LC MO p0,o(t) = l—B—Y l-d (4.32) 64 From the definition of the bivariate generating function given in Equation 4.18, the generating functions of the marginal distributions for state 1 and state 2 ,wl(x,t) and w2(y,t) respectively, may be specified. L [x8 + l - 8] O 41(Xst) w oaa spas Hoooz nod Am opopmv coaccaaaoa ooosflasamsu.m.m oaswam .mfluaa.m acoaccsom mean: cocooz mafia OOH H m H HIA pas: oEHuV 0H Ga m ON N m .3 E 1|. OH I H HIAM#HC3 GEHPV OH Ed I Q "mopwm coaufimcmne OOH 1 m opmum ooa I a oompm "macauoasaom HmeHcH O llV'l‘TIITY‘TI[I'llIllII‘VIIIIIlrnll'll'I'V'I'Illl‘IIIV'IUII[llllll'l'llllIVYIYUlll'llllIUIVUllllllI‘V'llllm UI‘III [[TTUVIIIII'II om 0: ow om OOH omH a 94949 UT uses JeqmnN 76 The only event left to generate is the time interval v between the transitions from state 2 to state 3. Be- cause this transition has the same characteristics as the one in Chapter 3 used to derive Equation 3.15 which gives the time interval between events, Equation 3.15 may be used to find v. Therefore, -loge(RN) 1 V "' “W- (5.18) a 2 where RN is defined as before. This simulation procedure will have an error because i a transition to state 2 is not taken into account in determining the next transition from state 2 to state 3 until the next transition of the latter type occurs. This error will be insignificant if the reduction to zero of the number in state 1 is somewhat sooner than for state 2. By using the above method, the complex form of the distribution for the number in state 2 (Equation 5.7) is avoided. However, the distribution of the time interval between a transition from state 1 to state 2 or from state 2 to state 3 involves a very simple form of the Joint probability distribution given by Equation H.30. But this approach would also incur errors of the type mentioned in the proceding paragraph. A simulation was achieved with a computer program as follows: 77 Initialize the program for time zero. Generate the time of the first transition from state l by Equation 5.17 and the time of the first transition from state 2 to state 3 by Equation 5.18. Compare the time of the next transition from state 1 to the time for the next transition from state 2. If the latter type occurs first, proceed to step 6; otherwise continue to step 3. Generate a random number uniformly distributed 'lfl'. . on the interval (0,1). If this is less than the value of hl/(hl+gl), proceed to step 5; otherwise go to step 4. When 2 reaches zero, the time of the next transition from state 1 should be set at a very large number so as to avoid any more transitions of this type. Decrease 2 by one and generate the time of the next transition from state 1 by adding the inter- val u (obtained by another evaluation of Equation 5.17) to the current time in the simulated pro- cess. Then return to step 2. Decrease 2 by one and increase m by one; record the time and size of m. Again obtain a new value of the interval u from Equation 5.17 and determine the time of the next transition of this type as in step 4. Then return to step 2. 78 6. Decrease m by one, record time and size of m. Generate a new value for v (Equation 5.18) and determine the time of the next transition from state 2 to state 3 by using this value. Continue this procedure until the number in state 2 reaches zero. if An example of the results of this type of simulation is shown in Figure 5.3. To simulate laboratory data for this model, the method discussed in Chapter 3 (p. 38) may be utilized. The indi— J vidual data points were generated by assuming a normal distribution as an approximation of a binomial distribution. Because the generating function (Equation 5.6) is the pro— duct of two binomial generating functions, the distribution (Equation 5.7) is the sum of several binomial factors. Consequently, a normal approximation can also be used to approximate the distribution. A random data point mr(t) can then be generated by the following equation mrm = fn' + v*om(t) (5.19) where E(t) and om(t) are the mean and standard derivation defined in Equations 5.8 and 5.9. The variable V is a normal deviate which can be generated by Equation 3.10. Again the approximation of the distribution by a normal will be without significant error if the variance is sizable. .q _ I. .. rcin an...» t 8.2.844?- 1: .COHuHmcmmB comm smozpmm me>anCH mEHB mCprnmcmO uo wonpoz u-mmpmpm mHan> oze npfiz Hmcoz sou Hm mumpmv cOHpmHsdom cmpsHssHmnu.m.m magmas mEHB om o: om om OH F. » _pLL _rkpcr.p c__.tht_f_L.____ _ » Ilrl|l 0 rs ON 9 7. 0: OOH a H OO thmpHc: mEHuO dalmm I m ON m mC3® Illlllofl Huh pH sHuO OH ea O OO OH I H HnnmpHcs mEHpO dwlmm I n mmpmm coHuHmcmsB OOH OOH n m mumpm OOH I H mumpm ”mCOHumHsoom HmeHcH omH a 91913 UT aJeq JeqwnN 80 An Example of Change of Death Rate The three state model of Chapter U is also proposed by Terui (1966) and used by Komenushi, Takada and Terui (1966) to represent a change of the death rate constant of Bacillus pumilus spores during heat sterilization. These researchers assumed that a spore may exist in two * viable states,both which will germinate using standard culturing techniques. They theorized that spores in one state could make a transition to a state of lower or higher resistance to heat sterilization. The change is E considered to be a discrete molecular change. Thus the viable organism can exist in only one state or another; there is no continuum between the possible states of existence. Prior to the work of Terui, Komemushi gt_al. (1966), Scharer and Humphrey (1963) had proposed a similar model to account for the non-logarithmic order of death curves for Bacillus stearothermophilus. But they gave only two possible transitions: one transition from an initial viable state to a second viable state and a transi- tion from the second viable state to a non-viable state. On the other hand, Terui (1966) included direct transitions from both viable states to the non-viable state. Neither Terui nor Scharer suggests probabilistic models, but they give the deterministic equations. The stochastic model for this example would be described by the generating function given in Equation 81 H.U5. Adopting this generating function, wk(z,t), for the homogeneous case, Lo Mo wk(z,t) = [z(b+c) + l - b - c] [za + 1 - a] (5.20) This model gives the distribution for the total number in states 1 and 2 or k = 2 + m. As derived in Equations a. Hf. .41 I M.u6, u.u7 and 4.U8, the probability distribution of k is the following k L L —1 M M -k+i pk(t> = z [ °}(b+c)1(1-b-c) ° [ ° ]ck-i(l-c) 0 i=0 1 k-i (5.21) with mean E(t) and variance ok2(t) given by E(t) = Lo(b+c) + Moa (5.22) ck2(t) = Lo(b+c)(l-b-c) + Moa(l-a) (5.23) The values of a, b and c are the same as given at the beginning of the chapter in Equations 5.2 through 5.5. The simulation of this process can be accomplished with the same computer programs presented in the previous sections. The only change required in these programs is the recording of the total number in states 1 and 2 in- stead of Just the number in state 2. Figures 5.4 and 5.5 give examples of this type of simulation. The method to simulate time intervals between transitions (p. 73) was 82 h1l.’l.fr.ll w :3 ...o 4 .OHIOH.m mcoprsom wchO soaps: -ammpmpm mHOsH> oza OOH: HmOoz you mmmoowm spawn OmpmHssHmII.O.m msOmHm msHB Om O: om ON OH O _ILLL._t...L_+.L.FfltilLIttrk.rL_It..ct..._L_LLLIYITL.ILILL_LTOFL.LI. O .OIIJ.. - ON JFK... I H H -3 Hz. , H OO WOO w OOH W ONH OOH H W m a: s u w w HIH OH O HOV OH :O W ON I N W OOH HIHOOHOO meHOO OH as I O V OH I H w HIHOOHOO meHOO_mHImH I O . OOH mopwm coHpHmcmnB - .. omH OOH I m madam . OOH I H mumpm - ”mQOHpmHsaom HmeHcH - oom 2 Due I 893823 u: JeqwnN Iaaom w:.;:..,;. s .COHpHmcwna comm coo3pmm mHm>anCH mEHB MCprpmcmO mo nonumz IIOOOOOO OHOOH> oza OOH: HmOoz you mmmoopm OOOOO OOOOHOEHOII.O.O OOOme 83 08.2”. Om OO Om ON OH O . . _ r O r _ t t H _ f# t _F _ p L t r L L b LL _ h . O - (I -4 #flf .I ON // - OO _JJIF WI r... s OOH I H H w HIAmpHc: mEHuO OH Ca I w J. W ON m x H OO HIHOOHOO OOHOO OH OO n O W OH I H H OOH HIHOOHOO OOHOO OH OO I O mouwm COHpHmcmpB T ONH OOH I m manpm - OOH I H mumpm mCOHpmHsaom HOHpHcH - OOH - OOH - OOH . OON a pue I 991948 at JeqwnN tenet 84 used in Figure 5.5 and Figure 5.4 was plotted from the same program used to plot Figure 5.2. In addition, the generating function (Equation 5.20) is the product of two binomial generating functions. So the normal approximation of the distribution may be used in the same manner as in the previous section (p. 78) $3 to generate random data points for the distribution given ;A there. Analysis of Model As compared to the deterministic model, the stochastic L- model provides a more accurate description of the real pro- cess. It defines the population size only at the discrete levels of possible existence and predicts the influence of random variables on the individual organisms. The model of this chapter may be compared with the equivalent deterministic model in the same way as was done for the elementary homogeneous model of Chapter 3 (p. 41). Without giving the details of the analysis, the same con- clusion as found in Chapter 3 was reached for this model. Namely, the deterministic model is sufficient for pre- diction of the lethal dosage required for a probability of sterility greater than 0.99. To test the stochastic model of this chapter, the same procedure as was given for the elementary model (Chapter 3) is required. But experimental errors must be small enough so as not to overshadow the intrinsic 85 variation of the process. This requirement is very difficult to obtain, and it is beyond the scope of this study to search for better methods for testing micro- organism death processes. As an example of this model, the data of Shull gt_al. (1962) were considered and found to have very large variations, at least too large to verify the stochastic model. However, the data were analyzed to ascertain if they fitted the generalized characteristics of the model. This resulted in the need to estimate the transition rates from the experimental data. Chapter 6 examines this problem. CHAPTER 6 DETERMINATION OF PARAMETERS FOR TWO STATE DEATH MODEL F“; Statistical Methods to Determine Parameters Before any experimental evidence of the homogeneous model of Chapter 5 can be considered, numerical values g for hl, h2 and g1 are required. Unless these parameters can be isolated for measurement, all three must be deter- mined from the same set of data. In this chapter, the determination of these parameters from the same set of data is considered for the heat activation model in Chapter 5 where the experimental data would be expected to give the number in state 2. In-addition to the trans— ition parameters, the initial number in state 1, Lo’ may be unknown too. For the elementary homogeneous model (Chapter 3), only one parameter, the death rate constant, was required. This was easy to obtain by a linear regression (logarithm of population versus time). But the model considered here can not be reduced to the form of a polynomial. If. 86 87 it could, standard statistical methods could be used. However, the problem of finding four parameters can be reduced to a problem of finding two parameters for cer- tain conditions. If h2 < hl,the slope of a semilog plot of the mean of the process approaches -h2 as time be- comes arbitrarily large. The mean of the process was ‘ I'-‘ given in Equation 5.8 as: L h — o l m(t) = —————:—— [exp(—h t) - exp(-(h +g )t)] ; I + MO exp(-h2t) . (6.1) Therefore, L h — o l m(large t) = —————:—— + M exp(—h t) . (6.2) If t is zero in Equation 6.2, the ordinate intercept of the line asymptotic to the mean line for large t is ob- tained. Thus, Lohl h1+31'h2 Ordinate Intercept = + M (6.3) 0 Consequently, if a semilog plot of the data becomes linear for large t, h2 can be obtained from the slope of the linear portion and one of the other three parameters determined from Equation 6.3. 88 Assuming that hl and g1 are to be determined by statistical analysis, Lo will be determined from Equation 6.3. To find hl and g1, the method of moments, the method of maximum likelihood, and the least squares method were considered. The latter method was chosen because it was the least complex of the three to apply to the model. A logarithmic transformation of the data was introduced as is usually done for models involving exponentials. Therefore, the least squares equation for this model is: Sum of Squares = E[In m(ti) - In D(ti):|2 (6.4) i where D(ti) represents the experimental population size (state 2) at time ti. The parameters hl and g1 are to be determined so that Equation 6.A is a minimum. By taking partial deriva- tives of Equation 6.4 with respect to g1 and hi, the normal equations are obtained. Designating these equations NG and NH respectively, In m(ti) - In D(ti) ar-n’(ti) NG(h ,g ) = 0 = 2 (6.5) 1 1 1 E(ti) 881 and In E(ti) - 2n D(t1)— aim-(ti) NH(h1,gl) o = z _ ah (6.6) i m(ti) l ' 1B ‘Kwtfl _I..".flkq I 'vi ‘ , _J' 'v 89 where am(ti) Lohl ______ O _____——— t exp[-(h +g )t agl hl+gl-h2 [:1 1 1 1] [5Xp(-h2ti) ' exP(-(hl+gl)tifl - (6.7) 111+gl-h2 15* If ‘ and afi(ti) = LO (gl-h2)[exp(-h2t)-exp(-(hl+sl)ti)l g. ahl hl+gl—h2 hl+gl-h2 + hltiEXp[-(hl+gl)tii] (6.8) Numerical Methods to Solve Normal Equations Since Equations 6.5 and 6.6 are non-linear, explicit solutions are not possible and numerical methods are re- quired. To simultaneously solve these non-linear equations, the two general procedures given by Hildebrand (1956, p. A50) were investigated. The first of these was the method of "successive substitutions." For this method, functions Fl(hl,gl) and F2(hl,gl) were defined such that k+l 1 (h k k.) (6.9) h = F l + and glk l = F2(hl ,glk) (6.10) 90 where the superscripts indicate the number of iterations. If hlk+l and g1k+l converge to the solution of the normal equations as k increases, Equations 6.9 and 6.10 can be solved successively to reach this solution. The functions Fl(h ) and F2(h ) can be defined 1981 1981 several ways to obtain a convergent series. Two common i R A: methods, the Newton—Raphson method and the method of HIM“ I‘d false position (Regula Fulsi) were chosen. Ostrowski (1960) gives an extensive analysis of both of these pro- cedures. For the method of false position: h k'lNH(h k,g k) — h kNH(h k’1,g k"1) F (h k g k _ 1 1 1 1 1 1 ’ ‘ - —1 l l l NH where ‘2 ll L M o o (01+1) (02+1) (C.15) M L (a3+l) O(au+l) O (C.16) M0 L0 [ya+l—a] [x8+yy+1-8-Y] (0.17) t exp(-f h2(I)dr] (C.18) O t exp(-é [hl(T)+gl(T)1dr) (0.19) ft t 60 hl(t)exp(—é [hl(r)+gl(1) - h2(r)]dr] (0.20) REFERENCES 117 _l_ REFERENCES Abramowitz, M. and I. A. Stegum, ed. (1964). Handbook of Mathematical Functions. National Bureau of Standards, Applied Mathematics Series 55. Arley, N. (1943). On the Theory of Stochastic Processes and Their Applications to the Theory of Cosmic Radiation. G. E. C. Gads, Copenhagen. Bachelier, L. (1900). Théorie de la Speculation. Ann. Sci. Norm. Sup. 3:21-86. Bailey, Norman T. J. (1964). The Elements of Stochastic Processes. John Wiley & Sons, New York. Bartholomay, Anthony F. (1957). A Stochastic Approach to Chemical Reaction Kinetics. Ph.D. Thesis, Harvard University. (Unpublished). . (1958). Stochastic Models for Chemical Re- actions. Bull. Math. Biophs. 19:175. . (1962). A Stochastic Approach to Statistical K netics with Application to Enzyme Kinetics. Biochemistry. 1:223. Bharucha-Reid, A. T. and A. G. Landau. (1951). A Sug— gested Chain Process for Radiation Damage. Bull. Math. Biophys. 13:153. Box, G. E. P. and Marvin E. Muller. (1958). A Note on the Generation of Normal Deviates. An. of Math. Stat. 19:610. Brockmayer, E, A. L. Halstrom and A. Jensen. (1948). The Life and Works of A. K. Erlang. Copenhagen Telephone Co., Copenhagen. Cohen, Abraham. (1933). An Elementary Treatise on Differential Equations. D. C. Heath and Co., New York. . Dewey, D. L. (1963). The X-ray Sensitivity of Serratia marcescens. Radiation Research. 12:64. 118 119 Doob, J. L. (1953). Stochastic Processes. John Wiley & Sons, Inc., New York. Einstein, A. (1905). Uber die von der Molekular Kinetischen Theorie der Warme Geforderte Bewegung von in Ruhenden Flussigheiten Suspendierten Teilchen. Ann. d. Physik. 11:549. . (1906). Zur Theorie der Brownschen Bewegung. Ann. d. Physik. 19:371. (1956). Investigations on the Theory of the Brownian Movement. Dover, New YorkT’ (Includes translations of 1905 papers). Feller, W. (1939). Die Grundlagen der Volterraschen Theorie des Kampfes ums Dasein in Wahrscheinlichkeitstheoretischer Behandlung. Acta. Biotheoretica. 5:11. . (1957). An Introduction to Probability Theory and Its Applications. Vol. 1, 2nd ed. John Wiley & Sons, Inc., New York Fredrickson, A. G. (l966a). Stochastic Models for Stegilization. Biotechnology and Bioengineering. :1 7e ' . (l966b). Stochastic Triangular Reactions. Chemical Engineering Science. 11:687. Forsyth, Andrew R. (1885). A Treatise on Differential Equations. MacMillan and Co., New York. Frank, H. A. and L. Campbell. (1957). The Non-logarithmic Rate of Thermal Destruction of Spores of Bacillus coagulans. Appl. Microbiol. 5:243. Geyer, Frederick P. (1966). Stochastic Mathematical Models for Describing Research Phenomenon. Paper presented at the annual meeting of the American Society of Agricultural Engineers. June 1966, Amherst, Massachusetts. Gibbs, J. Willard. (1902). Elementary Principles in Statistical Mechanics. Charles Scribner’s Sons, New York. Hildebrand, F. B. (1956). Introduction to Numerical Analysis. McGraw Hill Book Co., Inc., New York. 120 Hoffman, Joseph G. (1957). The Life and Death of Cells. Hanover House Books, Garden City, New York. Humphrey, H. E. and J. T. R. Nickerson. (1961). Testing Thermal Death Data for Non-logarithmic Behavior. Appl. Microbiol. 9:282. Isaacs, M. L. (1935). Disinfection. Chapter 12, p. 217 in F. P. Gay (ed.) Agents of Disease and Host Resistance. Charles C. Thomas, Springfield, Illinois. Johnson, H. A. (1963). Redundancy and Biological Aging. Science. 141:910. Kelly, C. D. and O. Rahn. (1932). The Growth Rate of Individual Bacterial Cells. J. Bacteriol. 99:147. Kendall, D. G. (1948a). On Some Modes of Population Growth Leading to R. A. Fisher's Logarithmic Series Distribution. Biometrika. 99:6. . (1948b). On the Generalized Birth-and-Death Process. Ann. Math. Stat. 19:1. . (1949). Stochastic Processes and Population Growth. J. Roy. Statist. Soc. Ser. B. 11:230. Kolmogorov, A. N. (1931). Uber die Analytischen Methoden inudfir Wahrscheinlichkeitsrechung. Math. Ann. 10 : 15. Komemushi, S., N. Takada and G. Terui. (1966). On the Change of Death Rate Constant of Bacterial Spores in the Course of Heat Sterilization. J. of Fermentation Technology (Japanese). 99:543. Lotka, A. J. (1945). Population Analysis. Contribution: W. E. leGros Clark and P. B. Medawar, ed. Mathematical Theory of Evolution. Oxford Press, Oxford. Markoff, A. (1906). The Extension of the Law of Large Numbers to Dependent Events. Bull. Soc. Phys. Math. Rev. 99:135. Murrell, W. G. (1961). Spore Formation and Germination as a Microbial Reaction to the Environment. Contribution, M1crobial Reaction to Environment. Eleventh Sym. of Soc. for Gen. Micro. University Press, Cambridge, England. 121 Ostrowski, A. M. (1960). Solution of Equations and Systems of Equations. Academic Press, New York. Parzen, Emanuel. (1960). Modern Probability Theory and Its Applications. John Wiley & Sons, New York. Pearson, Karl. (1897). The Chances of Death. Vol. 1, Edward Arnold, New York. Rahn, O. (1932). Physiology of Bacteria. Blakistion's, Philadelphia. . (1945). Iniury and Death of Bacteria by hemical Agents. BiodynamICa, Normandy, Missouri. Rutherford, E., J. Chadwick and J. D. Ellis. (1931). Radiations from Radioactive Substances. University Press, Cambridge, England. Taylor, G. I. (1921). Diffusion by Continuous Movement. Proc. Lond. Math. Soc. Ser. 2, 99:196. Scharer, Jeno and A. E. Humphrey. (1963). A Mathematical Model for Explaining Non-Logarithmic Death Behavior of Bacterial Spores. Unpublished Report, School of Chemical Engineering, University of Pennsylvania, Philadelphia. Shreider, Y. A. (1964). Method of Statistical Testing; Monte Carlo Method. Elsevier Publishing Co., New York. Shull, J. J. and R. R. Ernst. (1962). Graphical Procedure for Comparing Death of Bacillus stearothermophilus Spores in Saturated and Superheated Steam. Appl. Microbiol. 19:452. Shull, J. J., G. T. Cargo and R. R. Ernst. (1963). Kinetics of Heat Activation and of Thermal Death of Bacterial Spores. Appl. Microbiol. 11:485. Smoluchowski, M. V. (1906). Zur Kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen. Ann. d. Physik. 91:756. Stumbo, C. R. (1965). Thermobacteriology in Food Pro- cessing. Academic Press, New York. Terui, G. (1966). Theoretical Consideration of Heat Sterilization Kinetics. J. of Fermentation Technology (Japanese). 99:534. 4......2’ 1,) 122 Watson, H. W. and F. Galton. (1874). On the Probability of the Extinction of Families. J. Anthrop. Ints. 4:138. Wiener, N. (1923). Differential Space. J. Math. Phys., Mass. Inst. of Tech. 9:131. Yule, U. (1924). A Mathematical Theory of Evolution Based on the Conclusions of Dr. J. C. Willis. Phil. Trans. of Roy. Soc., Ser. B. 213:21. Zimmer, K. G. (1960). Studies on Quantitative Radiation Biology. Translated by A. D. Griffith. Hafner Pub. Co., New York.