v e,. A STOCHASTIC APPRoACH T0 POPULATION BALANCE moms ' Thesis for the Degree of Ph D ' . MICHIGAN STATE UNIVERSITY ANILKUMAR NARAYAN KHARKAR ‘ 1973 '- * LIBRA R Y ' Michigan State University This is to certify that the thesis entitled A STOCHASTIC APPROACH TO POPULATION BALANCE MODELS presented by AN | LKUMAR NARAYAN KHARKAR has been accepted towards fulfillment of the requirements for PH.D. degreein CHEMICAL ENGINEERING Date ‘ 9'19'73 .m- 0-7639 ABSTRACT A STOCHASTIC APPROACH TO POPULATION BALANCE MODELS By Anilkumar Narayan Kharkar Population balance models are extremely important in the quantitative analysis of chemical process systems characterized by aggregates of matter, as well as ecological systems. A principal aim of this work is to cast these models into a very general stochastic framework. By and large the development of population balance models in chemical engineering thus far has been based upon purely deterministic considerations and analogies with molecular processes. In this dissertation a probabilistic approach has been presented to develop a common basis for the population balance models used in the residence time distribution analysis of process vessels, various particulate processes and a large class of problems in ecological systems. The development of the population balance models is based upon the theory of stochastic population processes. The state of each entity in a population is represented by a point in a n— dimensional Euclidean space and the dynamics of the population is characterized by birth, death and movement of the entities in a closed domain in this state space. Anilkumar Narayan Kharkar In the context of populations of reproducing entities, the process of crystal growth with secondary nucleation and the dynamics of biological populations with reproduction are of particular interest. A common feature of the crystallization process with the state of a crystal described by its characteristic length and many biological populations with individuals characterized by a measure of their maturity is that all the entities in the population are in a particular fixed state at the instant of their first appearance in the popula- tion as a result of reproduction, irrespective of the state of the parent. To account for the fact that identical entities do not necessarily grow at the same rate, the growth of each entity is char- acterized by a nonhomogeneous diffusion process on the interval [0,5], where a can possibly approach infinity. Some types of boundary behavior of practical interest are discussed. Backward diffusion equations for the probability generating functional and the first two factorial moment distributions, as well as forward diffusion equations for the first two factorial moment densities are derived for a population of reproducing entities with no external input. The diffusion equations for the moment distributions have been analyzed for some simple cases. The following results have been Obtained under mild conditions: i) a general solution for the first moment density in terms of an infinite series of eigenfunctions of the diffusion operator; ii) after a sufficiently large time the first and second factorial moment of the number of entities in any subinterval in [0,5] grow exponentially with parameters and 0‘1 2&1 respectively, where a1 is the dominant eigenvalue of the ———_—_w Anilkumar Narayan Kharkar diffusion operator; and iii) the coefficient of variation of the number of entities in any subinterval in [0,5] reaches the same constant value after a sufficiently large time. Moreover, this asymptotic value in (iii) has been shown to be inversely propor- tional to the square root of the initial number of entities in the entire population. An interesting facet of the diffusion process in the context of reproducing populations is that under certain conditions the location of the boundary a will determine whether the population will increase or decrease. A graphical solution for the determination of this critical value of a has been obtained in terms of dimensionless parameters for the case where a is an absorbing barrier. The diffusion equation for the probability generating functional has been solved numerically to obtain the extinction probability of a population. These results had not been obtained thus far for the particular cases considered. In general, a population will also have an input of entities from a source external to the population. Expressions have been derived for the probability generating functional and the first two factorial moment distributions for a population with a nonhomogeneous Poisson input. It is shown that in the case of a population of non- reproducing entities if the birth, death, and movement of entities in a n—dimensional Spatial domain can be characterized by a Markov pro- cess, and if the population receives an external input in the form A of a nonhomogeneous Poisson process, then the number of entities in any set in the spatial domain at any time is a Poisson-distributed random variable. The use of this result in validating some assump- tions behind the models developed in this work has been demonstrated by means of a simple experiment. IlII"""""f"'—————————‘—————————_"—————"—""'—‘——""""""-—""'----------I:|:IEIII=SI Anilkumar Narayan Kharkar The theoretical results obtained in this work are of con- siderable importance in a broad class of problems related to chemical engineering and ecological systems. Several examples of practical interest are also discussed along with the results. A STOCHASTIC APPROACH TO POPULATION BALANCE MODELS By Anilkumar Narayan Kharkar A THESIS. Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemical Engineering 1973 To VAISHALI, our 1 itt 1e angel . ii ACKNOWLEDGEMENTS The author would like to express his deep appreciation to his major professor, Dr. George A. Coulman and Dr. Peter J. Brockwell for their constant guidance and assistance during the preparation of this thesis. Gratitude is expressed to Dr. Robert O. Barr for his help during the course of this work. Thanks are also due to Professors Myron H. Chetrick and Donald K. Anderson for their participation in this work. Gratitude is expressed to Dr. Stanley G. Wellso for making the cereal leaf beetle oviposition data available for analysis, and to Professor Dean L. Haynes and Mr. Richard V. Connin for their guidance and active help in performing the experiment. Thanks are also due to the faculty and graduate students in the 'Ecosystem D9313“ and Management" project who provided the necessary biological inSights related to the interdisciplinary aspects of this Work. The author is indebted to the National Science Foundation (Grants GK~1366 and GI-ZO) and the Division of Engineering Research for Providing financial support during his graduate work. Sincere gratitude is due to my parents for their moral SUPPOrt and encouragement throughout my graduate studies. Special thanks are expressed to my wife Neela, for providing encouragement and for her patience while I spent countless hours thinking about Populations. TABLE OF CONTENTS Page DEDICATION .................................................... ii ACKNOWLEDGEMENTS ..... .............. iii LIST OF TABLES ..... ..u-coocoo-noo-aoo-o oooooo o ooooo henna-0.0.. Vi LIST OF FIGURES ..... ...... ........ ...... .............. . ..... .. vii CHAPTER I INTRODUCTION ................ ...... .... ........ ..... 1 II MATHEMATICAL BACKGROUND . ........ ................... 16 2.1 Point Processes ........... ..... ............... 17 2.2 Counting Measures .............. ....... ........ 19 2.3 Moment Measures and Moment Densities .......... 20 2.4 An Illustrative Example ..... ....... .... ....... 22 2.5 Generating Functionals ...... .................. 25 2.6 Stochastic Population Processes ............... 28 III THE IDENTIFICATION PROBLEM .... ........ ............. 29 3.1 Residence Time Distribution of Flow Vessels ... 30 3.2 Life—Span Distributions in Biological Populations ...... . .................. . ...... 32 3.3 A Dispersion Phenomenon in Crystallization .... 35 3.4 Generalizations 37 3-5 Selection of a Proper Model ..............-.-.. 39 IV DERIVATION OF THE DIFFUSION EQUATIONS 44 4 1 Description of the Process .................... 45 4.2 Boundary Conditions ...,......,................ 47 4.3 The Backward Diffusion Equations ....... ....... 50 4.4 The Forward Diffusion Equations ; ..... ......... 61 4 5 Comparison of the Backward and Forward Equations ...... 72 V SOLUTION OF THE DIFFUSION EQUATIONS .......... ...... 76 5 1 General Considerations .............. ........ 76 5 2 Solution of the Diffusion Equations for the First Moment ................... .............. 79 iv CHAPTER Page 5.3 Solution of the Diffusion Equations for the Second Moment ............... ................ 90 5.4 The Problem of "Critical Length" .............. 102 5.5 Solution of the Diffusion Equation for the PGF -- Computation of Extinction Probabilities 111 VI POPULATIONS WITH AN EXTERNAL INPUT ...... ........ ... 130 6.1 Description of the Process .......... ... ...... . 131 6.2 The Probability Generating FunctiOnal of the Process ....................................... 135 6.3 Moment Distributions of the Population ........ 141 6.4 External Input in Populations of Nonreproducing Entities ...................................... 144 6.5 Applications ............... .......... ......... 148 6.6 Justification of a Poisson Input ..... ......... 154 VII AN EXPERIMENT BASED ON THE THEORY .................. 161 7.1 Motivation .......... .................. ........ 161 7.2 Methods and Materials ...... . ........ .... ...... 162 7.3 Results and Discussion ........ . ...... ......... 170 VIII CONCLUSIONS ............ ........... ................. 181 8.1 Summary of the Contributions of this Dissertation ..... ....... .... .................. 181 8.2 Areas for Future Research ......... ............ 183 REFERENCES ....... ...... ...... .......... ...... ............. .... 187 APPENDIX A SIMPLIFICATION OF THE DIFFUSION MODEL FOR SOME BIOLOGICAL POPULATIONS ............................. 192 B NUMERICAL COMPUTATION OF THE "CRITICAL LENGTH" ..... 199 C COMPUTATION OF EXTINCTION PROBABILITIES FOR A POPULATION OF REPRODUCING ENTITIES ................. 204 D ANALYSIS OF EXPERIMENTAL DATA ................ ...... 208 Table 5.1 LIST OF TABLES Page Model Parameters for the Dynamics of Populations of Bluegills (Lepomis macrochirus) .................. ..... 110 Simulation Parameters for the Computation of Extinction Probability of a Population ................ 112 Results of Simulation No. l -- Extinction Probability for a Population Generated by One Ancestor .... ........ 116 Results of Simulation NO. 2 -- Extinction Probability for a Population Generated by One Ancestor ............ 118 Results of Simulation No.3 -— Extinction Probability for a Population Generated by One Ancestor ..... ... .. 120 Results of Simulation No. 4 -- Extinction Probability for a Population Generated by One Ancestor ............ 122 Results of Simulation No. 5 -- Extinction Probability for a Population Generated by One Ancestor ............ 124 Results of Simulation No. 6 -- Extinction Probability for a Population Generated by One Ancestor ............ 126 Results of Simulation No. 7 -- Extinction Probability for a Population Generated by One Ancestor ............ 128 Statistical Analysis of Cereal Leaf Beetle Oviposition Data with Mated Females ........... .................... 153 Statistical Analysis of Cereal Leaf Beetle Oviposition Data with Egg-laying Females .......................... 159 Important Steps in the Experiments with CLB Populations ........................................... 168 Cereal Leaf Beetle Population Data from Experiment A .. 171 Cereal Leaf Beetle Population Data from Experiment B .. 173 Statistical Analysis of Cereal Leaf Beetle Population Data from Experiment A ................................ 175 Statistical Analysis of Cereal Leaf Beetle Population Data from Experiment B ................................ 177 vi Figure 2.1a 2.1b LIST OF FIGURES Counting Measure for a Sample Realization of a Point Process ......................... ...... . ............... Product Counting Measure for the Sample Realization of the Point Process in Figure 2.1a ......... . ..... .... An Example of Lumping of Variables in the Individual State Space for a Biological Population ............ ... Graphical Solution of Equation (5.64) for "Critical Length” ...................... .. ....................... A Typical Sequence of the Daily Egg Input by a Single Cereal Leaf Beetle Female ... . . ........... . ... . .. Comparison of Plant Condition between Experiments A and B ... ..... .... ...... .... ............... ... ....... Iteration Scheme for Obtaining a Root of Equation (5.6 ...... ..... ..................................... Page 24 24 40 107 156 166 201 CHAPTER I INTRODUCTION The use of mathematical models in the analysis of dynamic systems has become a common practice since the advent of high-speed large—capacity computers. The mathematical models used in chemical engineering practice can be classified into three broad categories: (a) TranSport phenomena models, (b) Population balance models and (c) Empirical models. The transport phenomena models use physico~ chemical properties of matter such as diffusivity, thermal con- ductivity and viscosity to obtain a dynamic description of the and energy. Systems comprised of ensembles of discrete entities are described by the population balance models. The dynamic description of these systems involves conservation equations for ghg number of entities in the population. The empirical models are comprised of purely empirical equations or regression analysis of data to describe the system. In particular, the applications of the population balance models include the residence time dis- tribution analysis of imperfectly mixed process vessels (e.g., Levenspiel and Bischoff, 1963; Bischoff, 1966), and various particulate processes such as crystallization (Randolph and Larson, 1971), size reduction (Randolph and Durando, 1971), particle agglomeration (Hulburt and Katz, 1964), fermentation (Tsuchiya et al., 1966) and liquid-liquid extraction (Ramakrishna, 1972). 1 The polymerization processes can also be described by population balance models (e.g., Funderburk, 1969). Typically, the population balance models describe the entities in a population as being distributed in a n-dimensional Euclidean space E according to a density x(zl,zz,...,zn,t) at time t. The Spatial co-ordinates 21,...,zn represent the relevant properties of the entities such as geographical location, size, shape, activity in case of catalyst particles, etc. Thus, x(zl,...,zn,t)dzl...dz represents the number of entities in the n infinitesimal volume element dzl...dzn located at the point (zl,...,zn) in E at time t. The "movement” of the entities along the i-th co—ordinate (i = 1,...,n) of E is Characterized by a rate of change of position vi(z1,...,zn,t) in that direction. The rate of addition of entities to the infinitesimal Volume element dzl...dzn from a source internal or external to the population is expressed as the birth rate B(zl,...,zn,t) of entities per unit time per unit volume of E. Similarly, the rate of change in the number of entities in the infinitesimal volume element representing a permanent loss to the population is char- acterized as the death rate D(zl,...,zn,t). If a number balance is written for the entities in the infinitesimal volume element at (21,...,zn) at time t by equating the net rate of change in the number of entities due to birth, "movement” along any of the co- I'm}; ordinates of E5 and death to the net rate of accumulation of entities, the following conservation equation results (Himmelblau and Bischoff, 1968): ax(z1,...,zn,t) n a[vi(zl,...,zn,t)x(zl,...,zn,t)] —___——-_—-——-+ Z = i=1 521 at B(zl,...,zn,t) - D(z1,...,zn,t) . (1.1) The movement of the entities is restricted to a closed domain in E, enclosed by a boundary F. The behavior of the entities at the boundary determines the boundary conditions and the density x0(zl,...,zn,0) at time zero describes the initial condition for (1.1). In general, the quantities x, Vi’ B and D will be func- tions of other variables as well, Such as temperature, concentration of a chemical in the medium surrounding the entities, etc. In addition, B will be a function of x when the entities themselves act as sources (reproducing populations). The quantities vi and D may also depend on x. It should be noted that in the derivation of (1.1) it is assumed that all the entities at a given point in E move in any given direction at essentially the same rate. In the context of residence time distribution of a process vessel, (1.1) represents a case analogous to a plug flow vessel, where all the fluid elements at a given cross~section of the vessel move toward the outlet at essentially the same velocity. The study of nonideal vessels, i.e., vessels which cannot be characterized by perfect mixing or plug flow, commonly involves a tracer analysis of the flow system. A known input such as a pulse of a tracer is introduced in the inlet stream of the flow system or its scale model and the concentration of the tracer in the outlet stream is measured as a function of time. A mathematical model is then con- structed to give the same input—output relation as that obtained in the tracer experiment. Two distinct methods are available to achieve this: The first method characterizes the flow system as a combination of a number of perfectly mixed vessels connected in series (and parallel, if necessary). This approach is commonly known as the mixing cell approach. The second method aSSumes that the flow is essentially a modification of a plug flow caused by an effective "diffusion" of the fluid elements in the axial and trans- verse directions in the flow vessel. The models used in this method are commonly called the dispersed plug flow models. The equations for the concentration of the tracer in the outlet stream are dc 5 i _ . _ z \_' dt — k(Ci-1 - Ci)’ 1 — 1,2,...,m , t_>».~ ‘ (1.2) for the mixing cell model with only a series combination of the perfectly mixed (hypothetical) vessels of equal volume, and 3c(z ,z ,t) 02 320(2 ,2 ,t) 02 36(2 ,Z ,t) 1 2 =_L_ 1 2 +;La_ 1___1_..2__ _ at 2 52? 222 azz 22 azz ac(zl,z2,t) aC(zl,zz,t) _ r —————————~—— - r -—~-—-—-—--- - D (1.3) L 321 R 522 for the case of a diSpersed plug flow model for a cylindrical vessel. In (1.2) ci represents the concentration of the tracer in the i-th vessel, t denotes time , k is the mean residence time of the fluid in a single vessel and cois the inlet concentra- tion of the tracer for the first vessel. In the mixing cell models all the vessels in a series combination are often assumed to be of equal valume to keep the model as simple as possible. In (1.3) and z the quantities represent the axial and radial co- Z1 2 ordinates respectively, C(Zl,22,t) denotes the concentration of the tracer at a given point in the vessel at time t, 2 are the dispersion coefficients in the axial and radial directions respectively, rL and rR are the velocities in the two directions, D denotes the rate of loss of the tracer per unit volume per unit time due to chemical reaction, absorption, etc., and t denoteS' time. It should be noted that (1.1) thru (1.3) are essentially deterministic in nature. Since (1.1) represents the conservation equation for the members of the population in a very general sense, it can be seen that the equation can be used to describe any biological population as well. Indeed, the equation has been used to characterize populations of daphnia (Sinko and Streifer, 1969), 8 Species of worms which divide by fission (Sinko and Streifer, 1971), shrimp (Billups, et al., 1971) as well as an insect pest pOpulation (Barr, Kharkar and Lee, 1972). Equations similar to (1.2) and (1.3) have also been proposed (Weiss, 1968; Kendall, 1948) and used (Takahashi, 1968; Stuart and Merkle, 1965) to describe cell popula— tions. A detailed discussion of the applicability of the "mixing cell” and "dispersion” analogies in describing a general population is given later in the text. Mathematical modeling of biological populations has become all the more important today in view of the pressing need for mankind to engineer the effects of urbanization and industrialization on the ecological systems. Equations (1.1) thru (1.3)are commonly used in the modeling of chemical engineering systems and the methods of parameter estimation are quite well established. A direct analogy of these models with the models for . biological systems makes it possible to think of analogous methods 0f parameter estimation as well. Literature in the fields of ecology, applied probability and mathematical biology abounds with works on the dynamics of biological populations. The models for single—species population dynamics used in these fields can be classified into three broad categories: (i) models which treat all the members of a population as identical entities. These models usually consist of a single difference or ordinary differential equation; (ii) models which divide the population into distinct groups according to age, maturity, etc. (i.e., define a dis- crete state space to describe an entity) -- these models commonly use a system of difference or ordinary differential equations; (iii) models which allow the individuals in the population to occupy any arbitrary point in a n-dimensional Euclidean space E with the relevant descriptors such as age, location, etc. as co-ordinates (i.e., define a continuous state Space to describe an entity), as is done in the development of (1.1). These models involve the use of integral or partial differential equations. All the models can be further classified into stochastic and deterministic types. In the partial differential equation approach for a deterministic analysis, use is made of a form of (1.1) or (1.3). The integral equation approach uses the renewal theory (see e.g., Feller, 1966) quite extensively and is a very convenient alternative to models using (1.1) for the case where n = l and v1 is a constant. In more complex situations the integral equa- tions for the system lose their simple character. As far as the deterministic models are concerned, models of types (i) and (ii) can be shown (Kharkar, 1971) to be approximations of the partial differential equation approach. As regards the stochastic models, much of the work deals with characterizations of type (i), (ii), or simple cases of (iii) in which the individuals are characterized by only one descriptor, namely, the age of an individual. The use of integral equations is very convenient for obtaining a stochastic description of the age distribution of a population (in terms of moments of the age distribution), and is used in almost all cases. Relatively little work has been done in the stochastic modeling of biological populations characterized as in case (iii) with descriptors other than age. Adke and Moyal (1963) developed a stochastic model for a population diffusing on the real line to characterize the spatial migration of the population. Birth and death rates were assumed to be constant. Adke (1964a) extended the model to the case where the birth and death rates can be func- .tions, and later (Adke, 1964b) to the case where the movement of the diffusing individuals is restricted by absorbing or reflecting barriers. Sevast'yanov (1958, 1961) studied the extinction proba— bility of a population diffusing in a compact region with an absorbing boundary for the case where the individuals in the popula— tion reproduce by a branching process. Davis (1965, 1967a, 1967b) studied the asymptotic properties of a population diffusing in an abstract space, and multiplying according to a branching process, i.e., an individual was assumed to produce a random number of off- spring only at the end of its life. As in the model of Adke and Moyal (1963), Davis assumed that the offspring are in the same state as their parent. Radcliffe (1972) extended the results of Davis (1965) by allowing immigration into the population in the form of a nonhomogeneous Poisson process. Population models are also of relevance in some areas of physics. A great deal of work has been done on stochastic descrip- tions of processes involving transport of elementary particles (see e.g., Brockwell, 1966 for a class of related problems and bibliography). In a sense, these processes are analogous to the branching processes discussed in the previous paragraph: an elementary particle in motion imparts some or all of its momentum to other particles in its way. The new particles produced are in the same location as the parent particle. When a particle with energy 6 sets k new particles in motion after a collision, the sum of energies of all the new particles must be equal to 6- Thus even in the case where the parent elementary particle does not come to rest, it may be necessary to assume that the parent is replaced by a "new” particle with an appropriate energy level, when the energy distribution of the particles is being studied. While some problems of chemical engineering interest can possibly be described as branching processes and some phenomena in the biological world besides spatial migration can also be cast in this framework, this analysis is often too cumbersometo use. This is true of a large class of problems related to biological and chemical engineering systems, namely, the processes in which the parent and the offspring are in different states and the parent does not die while giving birth to an offspring. In a biological population of higher organisms this is the case when the individuals are described by their maturity in terms of their size, weight, etc. An example of a chemical process of this type is the production of nuclei in a crystallizer. Under certain conditions the crystals present in a crystallizer produce new nuclei ("offspring") without any effect on the size of the parent crystal. Theoretical studies on populations often concentrate on the asymptotic properties of the population, i.e., the behavior of the population for large values of thne, because an analytical treatment of the problem for short times is not possible in many cases with the currently avail- able methods. Also, the dynamics of the system are often described by assigning stationary transition probabilities to the dynamic changes in the population, thus resulting in a model with time- independent parameters. In contrast, the short—term behavior of a system may be of crucial importance in a problem of engineering interest when the model is used to evaluate control alternatives. In many cases the time-dependence of the phenomena underlying the dynamics of a process cannot be ignored. This is particularly true of ecological systems due to seasonal changes in temperature, solar radiation, etc. Although a general model of a population of engineering interest may be too complex to be solved analytically, a numerical solution of sufficient accuracy can usually be obtained. The need for population balance models which are cast in a Very general stochastic framework and can be directly applied to systems of engineering interest, has motivated this research. Such models are of great value in studying a broad class of problems related to the design and management of ecological systems. In this context the descriptors for the individuals in the population 10 may be size or weight in the case of fish and plankton in a lake, productive value for trees in a forest or orchard, maturity in case of insect pests and their parasites or predators in a farmland, or the level of acceptance of a new idea or economic status in human populations. To have a greater confidence in the management and control policies derived from the models of these systems one would like to know the expected behavior of the populations as well as a measure of deviations from the expected behavior due to random phenomena. A general stochastic formulation of population balance models is also of value in the analysis of systems encountered in chemical engineering practice. The importance of this approach to particulate processes has been recognized only recently (Curl, 1967; Katz and Shinnar, 1969). Probabilistic models of flow in packed beds, where the fluid elements are assumed to move in dis— crete jumps, have been used quite successfully for describing the residence time distributions for some time (see e.g., Levenspiel and Bischoff, 1963; Buffham et al., 1970; Schmalzer and Hoelscher, 1971; Srinivasan and Mehata, 1971). As indicated above, one use of a stochastic model is in pre- dicting the fluctuations in the properties of a population. For small populations the stochastic fluctuations are always important if the rate processes such as birth, growth and death are random in nature. A stochastic model enables one to determine how large a population must be for the random fluctuations to be unimportant. Data on many biological systems shows a considerable scatter and hence the observations are often recorded in the form of means and standard deviations. If the contribution of the errors due to IJIIIIIIIIIlIIn--....._........______________________i,i 11 crudeness of observations, subjective judgement of the observer, etc., to the scatter in the data is known, the variance in the repeated observations due to random phenomena alone can possibly be estimated and used as an additional parameter in testing the mathematical models. This particular facet of stochastic modeling is demonstrated in this work by a controlled experiment, which is discussed later in the text. Experimental data on particulate processes and residence time distribution analysis of process vessels often exhibits an appreciable amount of scatter. Since the population balance models currently used to characterize these processes are largely based upon purely deterministic considerations and analogies with molecular processes, these models cannot account for the scatter. In molecular processes the number of molecules involved in an experimental setup is almost always very large and consequently, the stochastic fluctuations are insignificant. On the other hand, a very much smaller number of entities (i.e., particles, fluid elements, etc.) are encountered in the experiments with particulate processes and residence time distribution analyses. This fact coupled with a stochastic nature of certain phenomena such as formation and growth of the particles or movement of the eddies of the fluid inside the experimental vessels is responsible for at least a part of the scatter in the data. Stochastic models for such systems provide a useful means for the theoretical analysis of such data. For the case of reproducing populations with no external input, the coefficient of variation Q/variance/mean) of the gmw A 12 population, which is often taken as a measure of stochastic fluctuations, is alwasy significant if the initial population is small, irrespective of whether the population becomes very large at a later time. Bartlett (1969) indicates this to be the case for age distribution in biological cell populations. A similar result is obtained in this dissertation for a more general des- cription of a population. This fact is quite important in the con- trol of many ecological systems. For example, in the biological control of insect pests, often a small number of parasites or predators is released in the infested area. Control of the pest is uSually achieved after a certain time lag required for the parasite (or predator) population to build up to a significant level. A stochastic model for such a system would enable one to estimate a priori the chance of success of such a biological control strategy. A similar situation may arise in chemical processes of industrial importance. The phenomenon of secondary nucleation, which is a significant factor in industrial crystallizers, serves as a good example to illustrate this. In secondary nucleation new nuclei are generated by the breakup of dendritic growth on the sur- faces of a growing crystal or release of micro-clusters of particles during the growth of the parent crystal. Thus, if the initial number of seed crystals is small, the fluctuations in the crystal size distribution will always remain significant. Kane (1971) studied the rates of secondary nucleation of ice crystals in brine. He used a stirred batch crystallizer with liquid isobutylene added directly to the supersaturated brine as a coolant. The experiment 13 was started with a single ice crystal and the number of crystals in the crystallizer (counted using a photographic technique) and the bulk temperature of the contents of the crystallizer were monitored throughout the experiment. The data on the induction period (i.e., the time required to form a sufficiently large number of crystals from the solution to significantly affect the concentra- tion of the brine so that the freezing point of the solution begins to change) shows a considerable scatter. Expressing the nucleation process as a linear birth process (cf. Parzen, 1962, p. 296), he obtained an estimate of the contribution of the stochastic nature of the process to the scatter in the data. His model was based on the total number of crystals. The rate of formation of new nuclei per parent crystal per unit time was assumed to be a constant. In reality, this rate may be a function of crystal size. The equations developed in this dissertation are applicable to more general situations where the rates of nucleation, growth as well as removal from the crystallizer are functions of size as well as time. By solving these equations one can obtain the mean and variance of crystals in any size range at any point in time. Other situations of industrial importance where similar considerations may apply include failure of a chemical or biochemical reaction due to the catalytic activity 0f trace quantities of an impurity and guowth of mutants of the micro—organisms in a continuous fermenter. Even in situations where one is not interested in the stochastic fluctuations in the population, a probabilistic approach to the problem offers a certain conceptual convenience in building postulates related to the mechanisms involved in the rate processes, 14 as well as the relevant properties of the population. In this context, the fraction of the total number of individuals having a certain property may sometimes be interpreted as the probability that an individual has the property. Similarly, the rate of occurrence of a certain event may be looked upon as the rate of the corresponding random phenomenon. In some cases such phenomeno— logical descriptions may not be possible in a deterministic frame- work. For example, the formulation of (1.1) is true only when the pOpulation is large enough to justify the use of the density func- tion x(z1,...,zn,t). It will be shown later in the text that this restriction can be removed if the quantity X(zl,...,zn,t) is looked upon as the density of the first moment measure of the popula- tion, i.e., when x(zl,...,zn,t)dzl...dzn is taken to be the expected number of individuals at time t in the infiniteshnal volume element dzl...dzn located at the point (21,...,zn) in the space E. The outline of this thesis is as follows: Chapter II is devoted to the mathematical background necessary for the development of population balance models in a general stochastic framework. In particular, a general probability space is defined for the population. Concepts of counting measures, moment distributions and generating functionals are developed, and the relation between the generating functionals and the moments is discussed. The analogy between different approaches to the residence time distribution studies and modeling problems related to other popula~ tions (e.g., crystals, insects, trees, cells, etc.) is discussed in Chapter III in order to cast all population balance models in a 15 common framework. Following this discussion, general considerations for the selection of a proper type of model are presented. A diffusion model for a population in which the individuals are described by one descriptor is derived in Chapter IV. The back- ward diffusion equation for the probability generating functional, as well as the backward and forward diffusion equations for the first two moments for a reproducing population with no external input are derived and the relevant boundary conditions are discussed. Solutions to these equations for some simple cases are presented in Chapter V. A general model for a population must also include an input of entities from a source external to the population, i.e., an input other than reproduction. This is dealt with in Chapter VI. Some interesting results for a population with external input and no reproduction are proved and their application to problems of practical interest is discussed. A controlled experiment performed to demonstrate the use of a result obtained in Chapter VI in testing the validity of some common assumptions behind the mathematical models developed in this thesis is described and discussed in Chapter VII. Chapter VIII is devoted to conclusions and recommendations for future research. fr. .00 891 am CHAPTER II MATHEMATICAL BACKGROUND To cast the problem of modeling a population into a stochastic framework it is necessary to consider the dynamic changes in the population as a stochastic process. A stochastic process [X(t): t E T} is defined as a family of random variables X(t) indexed by a parameter t varying in an index set T. In the pre- sent context, the index t refers to time and T usually refers to the set of nonnegative real numbers or the set of nonnegative integers. When all the entities in the population are identical, as is the case in the population models of type (i) in the previous chapter, X(t) is a nonnegative integer at any given time, represent- ing the total number of entities in the population. Similarly, for models of type (ii) where each entity in the population must belong to one of finitely many categories (such as age groups), X(t) is a finite vector of nonnegative integers defining a measure on the set of age groups. In models of type (iii) each entity can occupy any point in the n-dimensional Euclidean Space E, and an analogous description of X(t) will be in terms of measures defined on the subsets of the Euclidean space E. At each point in time the con- figuration of individuals in the space E is a random distribution. The mathematical tool for the analysis of such processes is the theory of stochastic population processes developed by Moyal (1962). 16 17 An excellent review of the work done in this theory is provided by Daley and Vere Jones (1972). 2.1 Point Processes For the most general description, the state of an entity in the population is described by a point 5 in an abstract space E, called the individual state space. In case of population balance models of type (iii) discussed in the previous chapter the individual state space is an n-dimensional Euclidean space. The population state space 6is the union of disjoint spacesfik, k = 0,1,...,m, where Ek is the k-fold cartesian product E X E X...X E, for k 2 1. By convention E0 represents a single point, corresponding to zero population. A point gk 6 Ek represents a pOpulation with k dis— tinguishable individuals. It is a k-tuple (51,...,§k) of points in E. On each Ek is defined the minimal c-field Bk containing all product sets e1 X e X...X ek, where each set e1 6 B, the 2 c-field of Borel subsets of E. The minimal o-field defined on the 00 population state space 6 containing B consists of all sub— 2 * k=1 ~k sets e of 6 whose intersections with each E are members of k B . A measure space (6,8) can thus be defined from the pair N ~k k (E,B) and its k-fold cartesian product (E ,B ), k = l,...,m. A probability measure 9 can then be constructed on the measure space (6,3) either directly, or by constructing measures P(k) on '"k m k ”k (E ,Bk) which satisfy 2 P( )(E ) = 1. k=0 k ~ Pk = P( )(ES (2.1) is the probability that the population consists of k entities. 18 The triplet (6,5360 constitutes a model of a stochastic population process, and is called a point process. Physical interpretations of 6, B and 9 are as follows: the population state space 6 is the collection of all possible states of all the entities in the population, the c-algebra 8 refers to the collection of all possible "events" (e.g., all possible combinations of individuals in a biological population being in all possible (denumerable) age groups) and 9 denotes the probability measure defined on the collection of all possible events in the population state Space 6. Similar interpretations can be given k (k) to E, B, P and Ek, B , P in the context of a single entity and a population of k entities respectively. If the population is comprised of indistinguishable entities, the population state space as consisting of unordered sets of points in E is a more appropriate basis for the description of the process. The subspace of 63 of populations with k individuals is denoted by E(k), and a point in E(k), which is an unordered 00 set g = {g1,§2,...,§k} represents a population with k indis— k tinguishable entities. The c-algebras B( ) (k) and as are comprised of Borel subsets of E k) and 68 respectively. The probability 00 measures P: and 95 on the members of B and BS respectively can be defined in the same way as for the case of distinguishable entities. There is an obvious natural mapping of 6 onto 65, and to each probability measure 98 on 63 there corresponds exactly one symmetric probability measure on 6. 19 2.2 Counting Measures An alternative to the point process description of a popula- tion is the counting process. For any subset A of the individual state space E, let N(A) be the number of entities in A. Formally, a counting measure N is any non-negative, integer valued measure defined on the class U of all Borel subsets of E. The elements of the space 6 define certain counting measures, namely, Nanak) = _ H "MW I (2.2) 1 where I(A|-) is the indicator function of the set A: 1 if x E A 1(A\s) ={ o if x r A . Note that all points in d which correspond to the same unordered set §(k) = {§1,...,§k} in 65 lead to the same counting measure (2-2)- Moyal (1962) has proved that (2.2) defines a one-to-one correspondence between 68, and the space N of all counting mea5ures on the class U. Every probability space (68,5g,9g) has an equivalent counterpart (%,5k,éh). It is useful to exhibit this last equivalence explicitly. If A1:~--,A are disjoint subsets of E whose union is E, the proba- r bility that there are k1 entities in A1, k2 entities in A2,..., and k entities in A , is r r (L) k1 kr PN{N(A1) = k1,...,N(Ar) = kr} = PS ((A1 x...x Ar )5) k k = —-———l‘—!——-—P(L)(A 1x...XA r) (2.3) k !...kr! s 1 r Re C0 wh de' 20 k. where L = k1 +...+ kr’ Ail is the ki-fold cartesian product of A1 with itself, and in the first expression on the right hand side, (C)S is the symmetrization of the set C, i.e., the union of all possible distinct cartesian products C1 X...X CL’ where k1 of the Cj's are equal to Al’ k2 are equal to A2, etc. There are L!/k1!k2!...kr! such possibilities, and thus the second expression on the right hand side is obtained. 2.3 Moment Measures and Moment Densities The expected number of entities in a measurable set A is L Q ~ _ mm=Emmn= jL :uMngNa“h=zm§NAxdlraA> E i=1 L=1 IIMB L 1 Moyal (1962) has shown that (2.4) can also be derived using (2.3). Note that M(A) is a non-negative, monotone nondecreasing and countably additive set function, and hence a measure on (E,B). M(A) is therefore called the first moment measure or the first moment distribution. When there are only a finite number of entities in every bounded set A with probability one, it is often pOSSible to define a density m(-) of the first moment measure in the case where E is a n—dimensional Euclidean space. Thus, the quantity m<§>d§ = M. (2'5) denotes the expected number of entities in the infinitesimal volume dzn located at g = (zl,zz,...,zn) in E. element dg = dzldz k k (E ,B ) by 2... Higher order moment measures are defined on taking expected values of the k-th product counting measure L Nk(A(k)l§(L)) = z 2: 1(A(k)\(gi ,...,gi )), ks n, (2.6) 1= ik=1 1 k where A(k) = A n E(k) for any A CZE and §(L) = {§1,...,§L}. The use of these higher-order moment measures is complicated by the fact that they contain ”mass” concentrations on subsets of Ek. For this reason it is more convenient to deal with a factorial moment measure M defined as the expected value of (k) (k) ~L _ (k) N(k)(A \E ) _ z z #zik I(A \gil,...,gik), L 2 k, (2.7) = N(A)(N(A)—1)...(N(A)-k). The two types of moment measures may be contrasted as follows: The sum on the right hand side of (2.6) is taken over all samples of size k from the population §(L). The sampling is done "with replacement", and (2.6) therefore contains Lk terms. In (2.7) the sum is over all samples ”without replacement", and there are (L)k = Z:%i;r terms. The convenient relationship m (k+j) ~. M (A00) = g (k+j) P (Am x EJ) (2.8) (k) j=0 k 5 holds. In the case of a n-dimensional Euclidean space, the k-th . ~k h factorial moment measure often has a den81ty on E when t e PFObability of having more than one entity in an arbitrary interval (§:§+d§) is 0(5g), i.e., this probability is of a smaller order 0f magnitude than 5g. When such densities exist, the point pro- cess is said to be orderly. illustrate where the e descriptor. entity, whi nonnegative closed or h of sizes). h-algebra bility dist entities be Similarly, ‘ in the siZe ’ the expecte W) is a £0! definin To 1 realizatiOn eHtities in that k: 4 is Recess“: Note that ti to (N(A)2) in Figure 2 N2Wild) = 2.4 An Illustrative Example At this point it will be worthwhile to consider an example to illustrate some of the concepts developed so far. Consider the case where the entities in the pOpulation are described by a single descriptor. To be specific, let the descriptor be the size of the entity, which may increase indefinitely. In this case E is the nonnegative real line and B is the o-algebra generated by open, closed or half open intervals on the real line (representing ranges (k) of sizes). The product space Ek (or E ) and the corresponding o-algebra Bk (or B(k)) are necessary to define the joint proba- bility distribution of k distinguishable (or indistinguishable) entities being in k overlapping and/or disjoint size categories. Similarly, the counting measure N(A) denotes the number of entities in the size range A and the first moment measure M(A) represents the expected number of entities in IA. Note that for a fixed A, N(A) is a random number. The product counting measure is necessary for defining product moments, as mentioned in the previous section. To illustrate a product counting measure, consider a sample realization of the process as shown in Figure 2.1a. There are four entities in the population, of which three are in the set A, so that k = 4 and N(A) = 3. To get the second moment measure, it 2 is necessary to evaluate the expectation E(N(A) ) = E(N(A) X N(A))- Note that the product counting measure on the set A X A is equal to (N(A)2). Applying (2.6) to the realization of the point process in Figure 2.1a, 4 / 4 4 4 N2(AXA\g4) = 2 2 1(AxAi(§i,§j)) = r r I(A\§91(Ai§j) i=1 j=l 1:1 J=1 The measure pairwise co sample real the range occurrence realization corresponds correSpondi illustrated in Figu in the part is in no waj f(§)d§ den range (U? Probability (§Z,§2+dg2) having a p0 of hailing a aquantity . Probability tion of the diagmal of Hg) with ._ second [Home] the em“: ie: in). The: a Slim 0f tWI 23 [E meaSure N2(A X Aigk) thus has a contribution of 1 from every iirwise combination (gi,§j) of the sizes of the entities in a imple realization of the process, when 51 and gj both lie in ie range A. The combination (gi,gj) refers to the joint :currence of an entity at gi and another at gj. The sample :alization of the process illustrated in Figure 2.1a thus irresponds to N(Az) = 9. Different pairwise combinations (gi,§j) >rresponding to this particular realization of the process are Llustrated in Figure 2.1b. Note that corresponding to every point in Figure 2.1a there is always a point (gi,§i) in Figure 2.1b. l the particular case where the existence of one entity at g1 : in no way dependent upon the existence of another at g2, if :§)d§ denotes the probability that an entity lies in the size inge (§,§+d§) for an orderly point process, then the joint robability that an entity lies in (g1,§1+d§1) and another in ;2,g2+d§2) is f(§1)f(§2)d§1d§2. Thus, while the probability of lving a point on the diagonal of A x A is f(§)d§, the probability fhaving a point (g1,§2) in A X A is only f(g1)f(g2)d§1d§2 = o(d§), quantity of a smaller order of magnitude than the corresponding 'obability for the diagonal. This results in the "mass" concentra— (2) ~- the .on of the second moment measure on a subset of E Lagonal of any set A X A. This mass concentration has a density 1:) with respect to E, while the points of increase of the :cond moment measure corresponding to a pairwise occurrence of m entities in the size range A have a density with respect to 12) The second moment measure can thus be expressed in terms of Sum of two component densities -- a density with respect to E L Figure 2.1 Figure 2.1 24 k = total number of entities = 4 N(A)= number of entities in A = 3 g1 g2 1i3 :4 er ' x x x 3 4‘ ‘> g (size) |4 A 2:: 'igure 2.1a Counting measure for a sample realization of a point process Ii a? x x x x - n” ------- M x x x X (€4,§3) <1 X X X X (52:52) (€3,52) N2(AXA) = number f oints x x x x 9 P (slap (@331) MA = :1 AXA ! w ....... --— I I I I l I ' I l I ‘ I I £ $ = - A :I gl ii igure 2.1b Product counting measure for the sample realization of the point process in Figure 2.1a representir with respec the higher the correst k-th order concentrat: of factoria that of prc densities i is not orde triplets, e di, even tt on the diag t0 express tripler, et triplet, et Srinivasan tributions Ponulation mvolVEd wa 2.5 w A Q0 pmbability (k = 0,1,. . flinct ional bility gene 25 representing the mass concentration on the diagonal and a density ~<2>. with respect to E Similar mass concentrations also exist for the higher moments of the counting measure along the diagonals of the corresponding volume elements. For an orderly process the k-th order factorial moment measures (k = 2,3,...) do not have such ~(k) concentrations on subsets of E and thus the use of densities of factorial moment distributions are much more convenient than that of product moment measures.. A good discussion of the moment densities is given by Ramakrishnan (1950). When the point process is not orderly, i.e., when the probabilities of having twins, triplets, etc. in (§,§+d§) is of the same order of magnitude as dg, even the factorial moment distributions have concentrations on the diagonal subspaces of Ek. In this case it is still possible to express the moment measures in terms of densities if each pair, triplet, etc. are treated as distinct populations with each pair, triplet, etc. considered as single entities. Ramakrishnan and Srinivasan (1958) have illustrated this for the case of age dis- tributions in a population. However, even for the case of a population containing only singlets and twins, the mathematics involved was rather elaborate. 2.5 Generating Functionals A complete description of a point process in the form of the probability distributions Pék) or the moment distributions “(10’ (k = 0,1,...,m) can be obtained in terms of a single generating functional -- the probability generating functional. The proba- bility generating functional (PGF) is a natural generalization of the probab‘ variable tr space def it transform < random var: nonent gene tional are Let is defined “Pgégiem to be equal be rePlaced has Proved Probability form (2.9). B, then ran the PGF red 0f the rand 26 the probability generating function from an integer-valued random variable to a random variable taking values in a population state space defined in Section 2.1. Analogous extension of the Laplace transform of the probability distribution of an integer—valued random variable (or the moment generating function), namely, the moment generating functional (MGF), and the characteristic func- tional are also used instead of a PCP. Let 9 be a bounded complex-valued function on 6. Define k we; ) = o...e(§k) on Ek for every k, and therefore on 6. The PGF k °° k k 6(6) = E(W(§ )) = kE-O‘rfik e(§l)e(§2)-..e(§k)P( )(de ) (2.9) is defined for all meaSurable functions 6 such that sungEle(§)\ s 1. In (2.9) the product 6(El)...e(§k) is defined (k) to be equal to one when k = O. The measures P in (2.9) can (k) 5 without changing the result. Moyal (1962) be replaced by P has proved that there is a one-to-one correspondence between probability distributions PS on (65,83) and the PGF 0f the form (2.9). If (A1,...,Ar) is a finite measurable partition of 1.. E, then taking 90%) = . l ‘31 1(A1it), (2.10) IIMP'i l—‘ the PGF reduces to a multivariate probability generating function of the random variables N(Al)"°"N(Ar) r r N(A ) z: g.s(A.i§)lN(d§)] = E[ H g. 1] . (2.11) j=1 J J i=1 1 r C(e) = E exp 2 f log{ i=1 A, l he probab tributions hariation when in (2.12) defined as The MGF (Wi acteristic replacing Thus) for l TheSummat iii) lvith The k‘th P1 where 27 The probability distributions Pék) and the factorial moment dis- tributions M(k) can be expressed as functional derivatives (variations) of G: k A ... A = ' M(k)( 1 X X k) %:T 69 G(h), (2.13) when 9(E) = 1(Ali§)I(A2\§)...I(Akig) _ (2:14) In (2.12) and (2.13) the k-th functional derivative 6:G(n) is defined as (cf. Hille and Phillips, 1957, p. 111) k k d ngifigegi 6 C(fl) = i _ (2.15) e ark g—o The MGF (when e(§) has a nonpositive real part) and the char- acteristic functional (for all 6(§)) are obtained from (2.9) by replacing e(§) by exp(e(§)) and exp(ie(§)) respectively. Thus, for example, the MGF of a point process is given by k m k (k) k ¢(e) = Etexpr 2 e(§i))] = z j~k exp< 2 e(§i))P (do > . (2.16) i=1 k=0 E ' i=1 k The summation 2 e(§i) can also be written as an integral of i=1 9(5) with respect to the counting measure N(-) on E, i.e., k 2 e<€i) = i- e(§)N(d§) . (2.17) i=1 E The k-th product moment measure Mk is related to the MGF by A x A r l k k 0 2 Mk( 1 x... k) — (- ) 66 ¢( ), ( .18) where 9(a) = I-o-I(Aki§)- (2-19> 2.6 one The bilistic d triplet ( collection at a given the dynami probab ilit cess at ti of achievi interested in a given model desc usually SL1 Probabilit Of an Enti interval Possib1e t moment den functional tiOn of SL1 1" QValuat of having example, t 1”(ho 28 6 Stochastic Population Pr The discussion in this chapter thus far dealt with a proba- listic description of a population only in a static sense. The iplet (6,8,9) thus refers to the population state space, the llection of all possible events and the probability measure only a given point in time. For a complete stochastic description of a dynamic changes in a population one needs to know the joint )bability distribution of all finite combinations :1 2 l 1 ;s at times t., where it E 6. It can be seen that the problem 1 _ 4 ,...,EC } of sample realizations fit of the pro— k . 1 achieving this is formidable. Moreover, in practice, one is often :erested in the first few (usually two) moments of the population a given subset of the individual state space E, and thus a iel describing these moments of a counting measure on E is rally sufficient for most practical purposes. If transition >babilities can be assigned to describe the changes in the state an entity during a small time interval (t,t+6t) (or during the :erval (t,t+l) when the time parameter is discrete), it is often :sible to write down equations for the moment distributions (or rent densities) of the counting measures on E. The generating Ictionals prove to be extremely convenient tools in the deriva— In of such equations. The generating functionals are also useful evaluating another quantity of interest, namely, the probability having no entities in a given set A E E at any time t- For mple, this can be obtained from the PGF by setting e(§) = - I(A[§). Bef to identif and thus d change in etc.) with the next 8 the noveme in the des ment is ch is Markovi ContinuOUS tion has t etc., depe residEnce to illustr filaments l miXin8‘ce1 elements 1 In sure 31 System the as Well as to decide CHAPTER III THE IDENTIFICATION PROBLEM Before constructing a model for a population it is necessary identify the relevant descriptors of an entity in the population 1 thus define the individual state space E. The analysis of ange in the descriptors of an entity (e.g., maturity, productivity, 2.) with time, representing the movement of entities in E, is 3 next step in the modeling process. In a deterministic model, 3 movement of the entities in E is described by rates of change the descriptors, whereas in a stochastic description this move- It is characterized by transition probabilities when the process Markovian. Whether the population can be characterized by a itinuous movement of the entities in E or whether the popula- )n has to be grouped according to location, size, productivity, 2., depends on the particular system itself. The example of ;idence time distribution analysis of flow vessels may be cited illustrate this. Equation (1.3) describes the movement of fluid zments inside the flow vessel by a continuous motion, whereas a ring-cell model described by (1.2) essentially lumps the fluid ements in the vessel into a number of hypothetical "mixing C9115”- Some situations (1.2) is a more appropriate description of the item than (1.3). Analysis of data collected in tracer experiments well as the physical structure of the flow vessel are often used decide the type of model best suited for the purpose. The 29 analogies the models make it p0 of models 31 R_§i_d The is defined leaving th quantity stream whi the vessel that the a and H‘dt, V95591 wou outlet. T pmbabilit interpreta series?“ cOmPOUEInt . C2m are indepemen be express. Since the [mist equa1 30 alogies between the models for residence time distributions and 3 models for other populations, which were mentioned in Chapter I, he it possible to cast the methods of analysis of both the classes models in a common framework. i Residence Time Distribution of Flow Vessels The residence time distribution (RTD) C(t) of a flow vessel defined as the age distribution frequency of the fluid elements iving the vessel under steady flow conditions. Thus, the antity C(t)dt denotes the fraction of the fluid in the exit :eam which has spent between t and t+dt units of time inside a vessel. The RTD can also be interpreted as the probability it the age of a fluid element leaving the vessel is between t I t+dt, or the probability that a fluid element entering the Isel would take between t and t+dt units of time to reach the :let. The RTD C(t) is thus synonymous with the residence time )bability density. Sinclair and McNaughton (1965) have used this :erpretation in a discussion on the calculation of the RTD of a 'ies-parallel network of flow vessels, when the RTD of each iponent vessel is known. If two vessels with RTD's C1(t) and It) are connected in series and if the residence times are iependent of each other, the RTD C12(t) of the combination can expressed as a convolution t = C t- d , (3'1) C12(t) £01“) 2( r) r ice the length of time spent by an element in the composite system t equal the Sum of the times spent in each of the vessels. For aperfectl exponentia density (F In case of o functic Imp between tt avessel n the outlet passing; c zones in t relative t dead space simple tra model (Hm 33 packed aCteristic fOllow a 1 Overall ef elements I. the diffus neat along the increm interval var 31 perfectly mixed vessel the exit age of a fluid element has an ponential probability distribution, i.e., the exit age probability nsity (RTD) is E§R%§EAEL’ where a is the mean residence time. case of a plug flow vessel this probability density is a Dirac function, with a delay a corresponding to the residence time. Imperfectly mixed vessels have flow patterns intermediate tween the two ideal cases described above. The mixing inside such vessel may be such that a small fraction of the fluid can reach e outlet much faster than the remaining bulk, resulting in by- ssing; or a small portion of the fluid may be caught in stagnant nes in the vessel and reach the outlet after a very long time lative to the bulk of the fluid. The magnitude of bypassing and ad space can be quantified by analyzing the data collected in nple tracer experiments performed on the flow system or its Scale 561 (Himmelblau and Bischoff, 1968, p. 71). In flow systems such packed beds with size of the packing much smaller than the char- teristic dimensions of the bed itself, the fluid elements must [low a long tortuous path before reaching the outlet, and the arall effect may be similar to that of a "diffusion" of the fluid aments relative to the mean velocity toward the outlet. Use of a diffusion equation such as (1.3) can be justified when the move- it along a direction 2 is characterized by a Markov process and 2 incremental distance éz covered along 2 in a small time :erval (t,t+6t) can be looked upon as a random variable with mean E(éz) = r(z,t)6t + 0(6t), variance V(éz) = 02(z,t)6t + 0(6t)a (3‘2) and E(ézn) = 0(6t) for n 2 3 In many Ca sentation assumed tc vessels cc of k pet conbinatic vessels ha tion is th Thus, in t distributi Pro variety of of imperfe Med and faCt that amean, 3t Population where C ( Z rate 0f ma the Var iab 32 many cases the mixing cell models give a more realistic repre- tation of the RTD. In the mixing cell model, the flow vessel is umed to be composed of a number of hypothetical perfectly mixed sels connected in series and/or parallel. In the particular case k perfectly mixed vessels connected in series the RTD of the bination is a convolution of k exponentials and when all the sels have the same mean residence time a, the RTD of the combina- n is the density of the gamma distribution k-l em = W . (3.3) s, in this case, the exit age of the fluid elements has a gamma tribution. Life-Span Distributions in Biological Populations Probability densities of life-spans of individuals in a wide iety of biological populations show a similarity with the RTD imperfectly mixed vessels, and analogous models have been pro- ed and used for some biological populations. To account for the t that the age at which cells divide is randomly distributed about ean, Stuart and Merkle (1965) described the dynamics of cell ilations by the diffusion equation 2 2 g_ a cgzzt} _ r acgzzt} = acgzztz , (3'4) 2 az2 az at 'e c(z,t)dz is the expected number of cells with physiological between (z,z+dz), the drift coefficient r denotes the mean i of maturation and the diffusion coefficient g— accounts for variability in the rates of maturation among the individual cells. Ti fixed matr acell int the time s distribute mean resit justed to observatic approach c for RTD an connected tribution ageneral where a ce Phl’Slologi Stuart and pOpulation As far the ma and Ashfor 3qu as C0 trees into All the tr. to the next t 0 be equ ii 33 ells. The cells were assumed to divide when they reached a certain ixed maturity. Takahashi (1968) divided the complete life cycle of cell into a number of hypothetical growth stages, and assumed that he time spent by the cell in each of the stages was exponentially istributed with the same mean a- The number of stages and the ean residence time a for an individual stage could then be ad- usted to match the output of the model with the corresponding bservations in a radioactive tracer experiment. Takahashi's pproach can be seen to be an exact analog of a mixing cell model or RTD analysis where vessels with equal mean residence thne are onnected in series, and thus this model results in a gamma dis- ribution for the life-spans of the cells. Weiss (1968) proposed general deterministic formulation for cell population dynamics, here a cell was characterized by its chronological as well as hysiological age. This model is equivalent to the formulations of tuart and Merkle as well as Takahashi for the first moment of the )pulation. As regards other biological populations, a gamma distribution )r the maturation periods of insect life—stages was used by Read 1d Ashford (1968). In a study of productivity of perennial crops rch as cocoa and palm, Abkin (1972) divided the life-spans of the 'ees into a number of stages according to their productive value. 1 the trees in a given stage were assumed to have the same pro- ctivity and the time taken by a tree to mature from one stage the next was assumed to have a gamma distribution. In both these ses each life stage (of insects and trees) was implicitly taken be equivalent to a number of identical sub-stages with exponentially distributr an elastiI single spa growth ral Cooner, La occurring of food a1 models has general, : at which 1 0f old age in the (:85 to repress an abstrac individua] t0 the mat describe t With an at Population genetic va the Variab nentmned Currie (19 of the Cig the mean a of the lat life~span 34 distributed residence times. Fish populations are known to exhibit : an elastic growth. The size distribution of the fish born during a single spawning season widens with time, due to differences in the growth rates among individuals (see e.g., Cooper and Latta, 1954; Cooper, Latta and Schafer, 1956). The magnitude of the spread occurring in a given length of time gets smaller with smaller amounts of food available to the fish (Hall, Cooper and Werner, 1970). No models have been proposed to account for the elastic growth. In general, it can be said that in any vertebrate population the age at which reproduction starts and the age at which an individual dies of old age are randomly distributed about some average values. As in the case of fish, size or biomass of an individual may be used to represent maturity in some cases, or it may be necessary to use an abstract maturity scale [0,1], with the maturity of a newborn individual to be 0 and the end of the life cycle corresponding to the maturity 1 (Stuart and Merkle, 1965). It seems possible to describe these populations by a diffusion equation such as (3.4) nith an additional term corresponding to death. In many biological populations the life—spans are randomly distributed due to the genetic variations among individuals. Another factor influencing :he variability in growth rates is the availability of food, as mentioned above for the case of fish populations. Leftkovich and Iurrie (1963) studied the effect of food availability on the larvae )f the cigarette bettle, Lasioderma serricorne(F.) and found that ,he mean as well as the standard deviation of the maturation period if the larval stage of the insect increase with a decreasing supply f food. It is likely that food availability has an effect on the ife-span distributions of many other populations as well. industrial size distr when there seed cryst studied tl ditions tc Starting w fully moni they rneasr length of Size distr saturatior the ratio 0f the 812 were sligh such CUrVe With zero the indivj the fact t tributiOn were Mn.” Size relat CTlStal Si acteriZed CFyS ta 13 b 35 3.3 A Dispersion Phenomenon in Crystallization It has been observed that the product crystals from an industrial sugar crystallizer exhibit a much larger spread in the size distribution than that of the seed introduced initially, even when there is no appreciable nucleation during the growth of the seed crystals. Wright and White (1969) and White and Wright (1971) studied the growth rates of sucrose crystals under different con- ditions to specifically analyze the phenomenon of size dispersion. Starting with a uniform seed size in a batch crystallizer and care- fully monitoring the process to see that no nucleation occurred, they measured the crystal size distribution (using characteristic length of a crystal as a basis) at different times. The crystal size distributions obtained under different temperatures, super- saturations and syrup purities were plotted on a normalized scale -— the ratio of the deviation from the mean with the standard deviation of the size distribution. Although the shapes of individual curves were slightly different from one another, the average of twenty six such curves was very close to the density of normal distribution with zero mean and unit variance. The difference in the shapes of the individual size distribution curves was at least partly due to :he fact that sieve analysis was used to determine the size dis- :ribution for many samples. Since the shapes of individual crystals vere non-uniform, conversion of weight—size relation to a number— ;ize relationship resulted in an inaccurate determination of the :rystal size distribution. Let the growth of crystals be char- icterized by (3.4) with c(z,t)dz denoting the expected number of :rystals between sizes -(z,z+dz), the drift coefficient r represent: efficient individua: a point 2 (cf. Kimur ”flux" acr condition Because 01 must be f: the seconr The solut: size 20 of normal and Miller to agree p Whi due to an faces of (. effect On distributj 36 'epresenting the (ensemble) mean growth rate and the diffusion co- 2 efficient g— being equal to half the value of the variance of individual growth rates. The expected number of crystals crossing lpoint 2 per unit time can be shown to be c“ - e02 egg—’51, (3.5) :cf. Kimura, 1964, p. 185). Since there is no nucleation, this 'flux” across the point z = 0 must be zero. Thus, the boundary :ondition at z = 0 is * c (z,t) = 0 . (3.6) Lecause of the size dispersion phenomenon, the function c(z,t) mst be finite for a Sufficiently large z at all times. Thus, .he second boundary condition is C(m,t) = finite . (3.7) 'he solution of (3.4), (3.6) and (3.7) with all seed crystals of .ize 20 at time zero is approximately proportional to the density i normal distribution for a sufficiently large 20 (of. Cox ,nd Miller, 1965, p. 224). Thus, the diffusion model does seem 0 agree with the experimental observations. White and Wright postulated that the size dispersion occurs ue to an uneven distribution of lattice dislocations on the sur— aces of different crystals. These imperfections have a marked ffect on growth rates of crystals, thus re5u1ting in a random istribution of individual growth rates. ea L11 I u Ell di: 37 3.4 Generalizations In the context of a general population balance model where each entity is represented as a point in the individual state Space 1 u i, the random variations among individual rates of movement of the entities along any given direction in E may occur due to three distinct reasons: (i) (ii) (iii) The entities may be present in a uniform homogeneous environment, with the entities themselves having a random distribution in their potential capability to move along any given direction. The random life-Span distribution in biological populations occurring due to genetic variability falls in this category. The entities may be present in a uniform homogeneous environment and the process of movement along any direction itself may be random in nature. The phenomenon of size diSpersion in sucrose crystals is perhaps of this type, with the number of dislocations on the surfaces of each crystal itself varying randomly as the crystal grows. The entities and the process of movement along any direction (under uniform conditions) may be uniform, but the environment is heterogeneous and the entities spend random lengths of time under different environ— mental conditions due to the heterogeneity of the environment. This will be the case when there are regions of different supersaturation in a crystallizer and crystals spend random lengths of time in these 38 regions, giving rise to a size dispersion. A similar situation exists in the RTD analysis of imperfectly mixed vessels —— for example, the RTD of a packed bed differs from that of the ideal cases of perfect mixing and plug flow because the fluid elements spend random amounts of time in the crevices between the packing elements, resulting in little movement for random periods of time. 'iously, different combinations of these cases will occur in some uations. Depending upon the inherent variability among individual ities, their rates of movement in E, or the coarseness of the erogeneity of the environment, situations analog0us to dead space bypassing may occur in any general population. In a study of ferent animal populations in fresh water ponds, Hall, Cooper and ner (1970) introduced approximately a cohort of fish in each of ponds, and the size distributions of the fish were measured at end of the experiment. The size distributions show a long tail h a slight peak in the tail, suggesting the possibility that the h which get a head—start in their growth continue to grow at a ter rate than the others. This is analogous to the phenomenon bYPassing in a nonideal flow vessel, where a small portion of fluid moves much faster relative to the bulk. Similarly, in study of economic mobility or acceptance of new ideas in SOCial tems one may encounter a situation where a small portion 0f the ulation, pOSSibly representing certain ethnic groupS, progress accept new ideas much faster (or much leWer) than bulk of the 39 ilation, representing a situation analogous to bypassing (or 1 space) in imperfectly mixed flow vessels. The nonideal flow vessels and analogous systems involving er populatiOns may be denoted by a common term -- the nonideal ulation-flow systems. In the general case, the word "flow" ers to the movement of entities in the population along any given ection in E, with the 'directions' specified by the relevant criptors of an entity, i.e., the co-ordinate axes of E. It >uld be noted that in the general context, a perfectly mixed ‘sel is analogous to a situation where all the entities contained a region in E are assumed to be identical in all respects l the time taken by an entity to leave this region has an ponential distribution. This description thus involves lumping "averaging" the individuals in the appropriate region of E. ure 3.1 illustrates this lumping for a biological population h the relevant domain in E consisting of a two-dimensional graphical area and a range of maturities corresponding to the e-Span of an individual. Selection of a Proper Model The movement of the entities along the i-th co-ordinate axis of E can be characterized as a diffusion in that direction if movement of an entity along 21 is Markovian in nature and the .ance ézi covered in a short time interval 5t satisfies (3.2). implies that the rates of movement of the entitites along 21 a short period of time should have approximately a normal dis- ition. It is necessary to study the nature of the individual 40 Z3 (maturity) A birth death immigration —> -—-» emigration .L——_—————> growth gure 3.1 An example of lumping of variables in the individual state space for a biological population migration . death K9° v e 9 .oé’ $7 | ‘\. l _—— ~—————’Z I’ 1 " a V" 2" ' birth 41 iability among the entities, the process of movement along 21, the heterogeneity of the environment to see whether this con— ion is satisfied: When the inherent variability among the itiea alone is the cause of the randomness in the movement, i.e., n each entity moves at essentially a constant rate at all times h a random distribution of such rates among different entities, overall effect of movement of all the entities over a short 'iod of time may be analogous to a diffusion if the number of ,ities as well as their individual growth rates are such that : assumption of a normal distribution of movement of entities over :hort time interval can be justified. In the case where the )CESS of growth is random in nature, it would be necessary that a movement of an individual during a short interval of time be tracterized by a number of independent random increments, or as :ontinuous stochastic process. In the situation where the erogeneity of the environment is the cause of random movement the entities, it may be necessary to put a constraint on the rseness of the heterogeneity for the diffusion model to be tified -- for example, while the diffusion equation (1.3) can isfactorily describe the RTD in a packed bed where the char- eristic dimension of the packing is much smaller than that the bed, it is a poor model for short beds with large packing :icles or for fluidized beds. It should be noted that when the 2 Fusion coefficients g— for the n-dimensional diffusion process zero, the equation for the mean reduces to the case commonly in the population balance models in chemical engineering tice, namely, (1.1), where all the entities at a given point 42 ~ E move along any given direction at essentially the same e. When the conditions for a diffusion equation are not isfied, it may be necessary to lump the entities in different ;ions in E and specify transition probabilities (or rates) for : movement of entities from one region to another, as mentioned the previous section. The mixing cell models used in the RTD Llysis and the model preposed by Takahashi (1968) for cell popula- >n dynamics are examples of this approach. Experhments similar to tracer analysis in RTD studies, along :h a careful analysis of the potential causes of random movement the entities in the individual state Space would usually be :essary to decide the type of model best suited for the purpose. 7 example, the analysis of sucrose crystallization study of -te and Wright (1969, 1971) clearly suggests the applicability of 2 diffusion equation, whereas the growth rates of fish may need lore careful study to decide whether the "bypassing" effect in growth of fish populations is of significance. A stochastic description of a population characterized by iffusion process would rely on the theory of stochastic pOpula- n processes. The discrete cases such as the perennial crOp iuctivity model of Abkin (1972) will involve the use of discrete :ontinuous parameter Markov chains when the transition proba- ;ties for the movement of the entities from one state (char- rized by the productivity, maturity, etc. of the entities in a p) to another does not depend on the past history of the ties. For the purpose of the models considered in the next two 43 Japters it will be assumed that the movement of the entities in can be described by a diffusion process. A stochastic descrip- ion of the diffusion of the entities includes the stochastic ounterpart of the process described by (1.1), i.e., the case where .11 entities at a point in E move in any given direction at the ame rate, as a special case. CHAPTER IV DERIVATION OF THE DIFFUSION EQUATIONS As a Specific example of a general population balance model, he case where the individual state space E consists of the non- egative real line will be considered in this chapter. This case mplies to populations in which the entities are characterized by Lsingle descriptor z, denoting a property of an entity such as .ocation, age, maturity or size. The model can be readily extended :0 the case where E is a n—dimensional Euclidean space for which :he properties of an entity representing the co-ordinate axes are independent of each other. The total number of entities in a mpulation can increase due to reproduction by the existing entities n the population, or by virtue of an external input. In the lerivation of the equations in this chapter it will be assumed that lhere is no external input of entities to the population. Moreover, I ‘11 the "newborn" entities are assumed to appear at the origin of at the instant of their birth. The descriptions in this chapter ill have to be modified slightly in case of populations where the offspring" is in the same state as the parent, as in the spatial igration of a biological population via a diffusion process (cf. dke and Moyal, 1963). The effect of external input will be con- idered in Chapter VI. In the following discussion the creation of entity by reproduction will be called birth of the entity and 44 45 oval of an entity will be termed as death of the entity. Description of the Process At time s the population is started with a single entity an arbitrary point y E E. The process of movement of an entity any point z E E is assumed to be a continuous Markov process, 2., as the time interval 6t becomes small, so also does the ;tance 62 covered by the entity during that interval, and the :ure locations of the entity depend only on its present location, : not the past history. The distance 62 covered by an entity a small time interval (t,t+6t) conditional on its survival :ing this interval is a random variable with a density g(6z;z,t) 1 satisfies the conditions specified by (3.2), i.e., (mean) E(6z) = r(z,t)6t + 0(6t), = foam + out), (4.1) (variance) var(6z) 3(521‘) = 0(5t) for k 2 3. :ing the interval (t,t+6t) the entity at 2 has a probability :,t)6t + 0(6t) of dying and a probability xi(z,t)ét + 0(6t), new entities at z = 0 with l f 1,...,m, of producing i \ CO 2 ii(z,t) s x < co (4.2) ‘ i=1 1 is’ l \ ‘obability that an entity will not reproduce during (t,t+6t)} = xi(z,t)5t + 0(6t). (4.3) 1 : state 2 of an entity is unaffected by giving birth to an "off— m— 46 : process of change of position of an entity as well as its re- )duction and death proceed independently of other entities in the )ulation. In view of the independence of the entities, the case are the initial population contains k entities (k > 1) can be 1di1y modeled by a Superposition of k populations, each starting :h a single ancestor in the appropriate location corresponding to 3 initial distribution. Two obvious cases where this description applies directly 3 the dynamics of biological populations of individuals char- terized by maturity, and crystallization processes involving :ondary nucleation. In biological populations the probability of rth refers to reproduction and the probability of death refers to e processes of natural death, predation, parasitism and harvesting. crystallization processes, reproduction refers to secondary cleation and death corresponds to precipitation,or loss of the ystals from the crystallizer due to mechanical removal or bulk pw of the magma. Most particulate processes in chemical engineering htain entities which do not reproduce. It can be seen that these cesses correspond to the case where the birth rates xi(z,t) = O, 1,...,m. In processes such as degradation of polymer molecules, e reduction of solids, etc., every breakage of an entity normally duces a number of entities of different sizes, thus generating andom size distribution of the pieces. The models derived in s chapter can possibly be generalized to include this case, but analysis of this general case will be much more complicated. For the quantitative description of the population, a count— measure N(A,t|y,s), t 2 s, is defined as the number of entities 47 :he set A CiE at time t, given that there was one ancestor :he point y at time 3. Obviously, N(A,S\y.5) = I(A\y) (4-4) 1 probability one. In (4.4), I(A|o) is the indicator function the set A. Starting with (4.4) the description of the temporal ages in the population can proceed in two directions. One way describing the population dynamics is to express the counting sure for a fixed arbitrary set A and a fixed time t, as a ction of s and y. Since —w s s s t, and since the dynamics the population is expressed relative to the condition specified (4.4), i.e., with s = t as the reference point in time, it be seen that this is a description of the population obtained moving ”backward" in time. This description is therefore called backward description. Another way of describing the changes he population is to express the counting measure corresponding n infinitesimal set dz located at z E E for fixed arbitrary nd 3 as a function of z and t. Since this description 1ves moving "forward" in time with respect to the reference t t = 5, this is called the forward description. Accordingly, ) iS a forward description of the population because it des- es the population in an infinitesimal interval (z,z+dz) as nction of z and t. Boundary Conditions The movement of the entities is restricted to a domain in alled the spatial domain 5. The spatial domain may be a subset 48 E, or the entire space E itself. In the present context, .5 taken as the closed interval [0,5], where a may possibly roach infinity. In the derivation of the model, the behavior of entities at boundaries must be considered. The possible modes of boundary avior of an entity diffusing on the real line were considered ?eller (1954). Some possibilities of practical interest are: reaching a boundary, the entity may be instantaneously returned the interior of D, i.e., the entity is always restricted to stay the interior of D, resulting in a reflecting barrier. Another sibility is when the entity is removed from the pOpulation as n as it reaches a boundary, representing an absorbing barrier. nird possible type of behavior at the boundary occurs when an ity reaching a boundary stays there for an exponentially dis- )uted random length of time, and then jumps to a point in the arior of D, or to the other boundary. It may be added here that 1 E is a multidimensional Euclidean space, there is a fourth ;ibi1ity, namely, the case where an entity reaching a boundary is along the boundary according to a Markov transition proba- .ty. In the context of biological populations and crystallization esses described above, the boundary at z = O is a reflecting ier, because no entity can cross this boundary and moreover, moment an entity appears at z a 0, it starts to ”grow", i.e., S to the interior of D. When D is a finite interval ], one may encounter a reflecting or an absorbing barrier at The third type of boundary behavior, where an entity can take 49 ?inite jump into the interior of D may be encountered in rare ;es. It is intuitively obvious that as 5 approaches infinity, : state of the population in any finite interval at any finite 1e Would not be dependent upon the type of boundary behavior at 5. Fish are known to exhibit stunted growth if the population [sities are too high (Beckman, 1940). This may be interpreted to in that there is an upper limit on the size that an individual :h may attain depending upon the availability of food. This iiting size can possibly be viewed as a reflecting barrier re- ‘icting further growth of the fish. However, it should be noted It since the stunted growth phenomenon occurs due to a severe i ,eraction among individuals for the available food, the stochastic "mulation of the population process may perhaps be only a crude >roximation in this case. In an intensively exploited fishery a forest resource) few individuals above a certain size (or ductive value) may survive, enabling one to imagine an absorbing rier at a finite size (or productive value). Similar considera- ns will apply to any intensively exploited biological population. some biological populations all the individuals above a certain urity can be considered identical in all respects. In such cases is possible to consider this maturity to be an absorbing barrier assume that an individual reaching this maturity is removed n the original population to become a member of another population ;isting of identical individuals. An absorbing barrier can be ified in case of particulate processes in chemical engineering .tice when the particles precipitate out of the system, or are ved mechanically by some separation technique when they reach 50 ertain size, weight, etc. The flow of fluids through packed beds ers an interesting possibility of an application of the third e of boundary behavior. The movement of fluid elements near the s of the bed cannot be characterized by a diffusion equation n in long beds, because the random velocity patterns of the eddies essary for an effective diffusion are not fully developed in the zones. The distance covered by an eddy at a boundary (i.e., inlet or exit end of the bed) in a short time interval will not e a normal distribution as indicated by (4.1) -— instead, the tribution may perhaps look more like an exponential. To account the end effects it may be worthwhile to consider an idealized w behavior where a fluid element at a boundary is capable of tantaneously jumping to a point in the interior according to a wn probability distribution. When the volume of the headers of packed vessel is rather large, it may be possible to assume the ders as perfectly mixed (hypothetical) vessels at the b0undaries, ing rise to finite exponentially distributed times of stay at the 1daries. A detailed mathematical description and verification of : boundary behavior is beyond the scope of this dissertation. The Backward Diffusion Equations A stochastic model for the process described in Section 4.1 be derived here in the form of a backward diffusion equation the probability generating functional of the population. Let ,tly,s), t 2 s, denote the number of entities in [0,2] at t, given that there was one entity at y at time 5. Define 51 :(e,tly,8) = E{[H 9(zi)],t|y,8} = E{eXP INEIOg 6(2)]N(dz,t\y,8)} (4-5) i D m the conditional probability generating functional for the popula- ;ion at time t, given that there was one entity in the population It y at time s. The product on the right side of (4.5) has ane term e(zi) corresponding to each entity at zi at time t. n View of the independence of the entities in the population a CF conditional on k entities at y1,...,yk at time s can be hown to be equal to k G(e,t‘y1,...,yk;s) = .H1G(e,t\yj,s) . (4.6) J: e proof of (4.6) follows from the definition of the conditional ounting measure with k ancestors k N(A,t\y1,.-.,yk;s) = 2N(A,t\y..s) (4.7) j=1 J nd the second expression on the right side of (4.5) (Moyal, 1964). den there are no entities in the population at time t, G(e,t\y,s) quals one (cf. Equation (2.9)). As mentioned before, the spatial )main 5 is taken as the interval [0,5]. The case where 5 : the semi-infinite interval [0,m) is obtained by letting a vproach infinity. Consider the change in G(e,t1y,s-5S) brought about by a all change in the ancestor during the time interval (s-6s,s). view of the fact that the ancestor can move from y to y+6y the time interval (s-5s,s) with probability g(6y;y,s)d(5y) ovided it survived during the interval, ,t'y,S-SS) = m 5 - E xi(y.S)6S) (1-u(y,S) és)£g(6y;y,S)G(e,t|y+6y,S)d(6y) 1=1 m ,5 1'M(Y,S)6S) E ki(Y)S)6S G(e,t‘O,S) l£g(6y;y,S)G(9,t‘y+6y,s)d(5y) i=1 (y,s)€>s - 1, (4.8) ’ T at t ancestor P[ancestor neither died nor reproduced y at time s—ss = during (s-és,s)]§{P[ancestor moved to y+5y during (s-5s,s)][PGF at t ancestor at y+6y at time 8]} P[ancestor did not die during (s-és,s)]P[ancestor produced i offspring during (s-bs,s)][PGF at t with the i offSpring as ancestors at O at time S]I{P[ancestor moved to y+5y D—‘ during (s-és,s)]}{PGF at tlancestor at y+5y at time s} -P[the ancestor died during (s-és,s) without producing any _offspring][PGF for no entities in the population] . probability of the ancestor dying and also producing any off- ng in an interval 53 is 0(6s)#, and hence the corresponding is omitted from (4.8). Expansion of G(e,t\y+6y,s) in a or series about y gives 2 2 Hr%YJO =G(e¢lyfi)_pagfing¥gaéy4.AEQL%UJQ_%p BY + o#. (4.9) : function f(t) = 0(t) if lf(t)/t s k < m as t a 0 (e.g., [uantity E(6z) in (4.1)). The function f(t) = 0(t) if It a as t a O, i.e., f(t) is of a smaller order of magnitude t, as in case of E(5zk), k 2 3 in (4.1). 53 )stituting (4.9) in (4.8), using (4.1) and omitting terms of ier higher than 6s, one obtains 3,C‘Y:S‘5S) = C(e,t\Y>s) + 63EG;G(9,t\y,S) + u(Y:S) _ z li(Y:S)G(G:t\Y:S) i=1 + z xi(y,S)G(e,t\0,S)iG(e,t\y.S)], (4.10) i=1 are the operator d9 is defined as 2 2 c7 = [9—igafil 5~§ + r(y,s) a_._ n(y,s)] . (4.11) Y ay By >tracting G(e,t‘y,s) ’from (4.10), dividing by as and letting approach zero, the backward diffusion equation for the PGF co Gg9,t y,s) as = a§G(e,t\y,s) + u(y,S) - z ii(y,s)G(e,t\y,s) i=1 + z 11(y,S)G(e,t\0,S)iG(e,tly,S) (4.12) i=1 obtained. It may appear that if E(6zk) = 0(6t) (instead of o(6t)) (4.1) for say, k = 3, then one would be able to get a third er partial differential equation instead of (4.12) by truncating Taylor series expansion of G(e,t\y+6y,s) in (4.9) after the rd derivative term. However, Pawula (1967) has shown that if Z3) = 0(5t), then E(ézk) = 0(6t) for all k > 3, thus demanding infinite Taylor series expansion in (4.9) and consequently an inite-ordered partial differential equation. When y = 5 is an absorbing barrier, in view of the condition t every entity at the absorbing barrier is immediately removed 54 -om the population, i.e., the number of entities at the barrier ; always zero with probability one, the corresponding boundary .ndition is C(e,t\3,s) = 1 . (4.13) As explained in the previous section, the boundary at y = 0 z a reflecting barrier. For a reflecting barrier at y = 0 the Ltity at the boundary is assumed to jump instantaneously to an merior point 6y. Thus, G(e.tl0,s) = G(e.t\oy,s) nce, _ a§i_i£tii;L C(e,t\0,s) - G + by 9 By S ly=o W‘ = 0 . (4,14) ay y=0 milarly, if y = 5 is also a reflecting barrier, = 0 , (4.15) nce there is only one entity at y at time t = s, [H9(Z.)],t\y,t)} = e(y). Hence, the initial condition is given by 1 G(e.tly,t) = 9(y) - (4-16) If (4.12) can be solved analytically with the initial con— tion (4.16) and boundary conditions (4.14) and (4.13) or (4.15), E SOlution G(e,t\y,s) would contain a complete stochastic . I n l .S no SCrlPtion of the population With one ancestor when there 1 55 Lput of entities to the population from an external scurce. How- 'er, because of the nonlinear nature of (4.12) and the occurrence 7the arbitrary function e(y) in the initial condition, such Lalytical treatment is not possible. Remark 4.1: If there are k entities initially at >,y2,...,yk; then the complete PGF for the population is obtained ’ solving the backward diffusion equation for G(e,t\yj,s) for = l,...,k and substituting these in (4.6). Remark 4.2: In view of the definition (2.9) for a PGF, the obability of having zero entities in the p0pulation (i.e., the :tinction probability) is given by the solution of (4.12) with :y)=0. Similarly, choosing e(y) = <4”) C(0,t\y,s) and solving (4.12) one obtains the probability nerating function ngf = E{§N < > ) } H my Em )+ge( > tly s>P(k)} H . e . ,t ,s = 2 z. Z. , , i=1 Tl zl +g 21 \y 1 8g k=0 Ek i=1 1 1 00 k . l (k) k 1m ‘i E i H ( (z )+(Q+A)e(z.),tly.S)P (dz ) *0 A k=OIEk i=1 n 1 1 m k k 1 - 2} k n +ge,t\y,s)P( )] k=0 E i=1 w k . 1 (k) k 1m ‘{ 2 ~ H (n(z.)+ge(z.),t\y,s)P (dz ) *0 A k=OIEk i=1 1 1 58 m k a— C+ 2 (k) + 13::0'r‘iilkm 5§[1:11m(zi)+( hA)e(zi) ’tiy’5>1+0(4 )JP (dzk) co k (k) k _ Z N H (fi(z.)+ 9(z.),t y,s)P (d ) k=0J~Ek i=1 1 Q 1 \ z l are 0 s h s 1, by the mean value theorem. Hence, by the dominated nvergence theorem (cf. Parzen, 1962, p. 274), k 00 k {E[ H (N(Zi)+ge(zi),t\y,8)]} = 2 j~k{1im [a—< n (“(21) i=1 k=0 E A~O 55 i=1 + e(zi>,t\y,s>> + 0(A)]P(k)(dzk)} k _ a_ — E{5g [igl(n(zi) + ge(zi),tly,s)]} . In particular, taking A and A 1 2 to be the intervals [0,21] 1 [0,22] respectively, defining * M1(zl,t\y,S) = M1([0,z1],t\y,s), (4-21) * = 4.22 M2(Zl:22,t‘y,s) M(2) ([0:21] X [0,2213tly)s) ( ) i using (4.19) and (4.20), the backward diffusion equations for a first and second factorial moment measures of the population in ’21] and [0,21] X [0,22] respectively can be obtained. The :kward diffusion equation for the first moment measure is 71‘ 5M (Z :tlYfi) 0° 3‘.“ 1 * - 23 - = + l . M (Z at 0,5) (4‘ ) as dy M1(zlytly:s) 121 Al 1 1 i :h the initial condition M:(zl,t\y,t) = (4.24) -~~_.__-~ .. 1d the boundary conditions 9: 8M1(Zl,t\y,s) ay y=0 = 0, (4.25) 7': 7': M1(zl,t\y,s)\y=21+ = M1(zl,tly,s)\y=zl_ (4.26) 7': 9‘: 5M1(z1,tly,s) 5M1(zl,t\y,s) _—_ _ = —-———— _ (4.27) By y-zl+ oy y-zl_ d * N M1(zl,t\a,s) = O (4.28) r an absorbing barrier at y = a, or M*( i ) a z ,t y,s __l__l_____._ N = 0 (4.29) By y=a r a reflecting barrier at y = 5. The continuity conditions .26) and (4.27) at y = 21 for M:(zl,t\y,s) and its derivative th respect to y are necessary in view of the fact that the itial condition (4.24) essentially divides 5 into two regions, S y 3 z1 and 21 < y s 5. Conditions (4.26) and (4.27) follow om the continuity of G(e,t\y,s) and its derivative with respect y throughout 5. Similarly, the diffusion equation for the cond factorial moment measure is 5M*(Z ,z ,tly,s) __§__l__§_______ = a7M*(z ,z ,t‘y,s) as y 2 1 2 m m * Zi(i-1)xi(y,s)M:(zl,t\O,s)M:(zz,t\0,8) + Biki(y,S)M2(zl,z2,tiy,S) i=1 i=1 m. * a * * 3 iEixi(y’s)[Ml(zl,t\y,s)M1(z2,t\0,S) + M1(zl,t\0,s)M1(zz,t\Y,S)] (4. o) with the initial condition 9: M2(z1.22,t\y,t> = o (4.31) and boundary conditions * 5M2(zl,zz,t\y,s) = O 4.32 By M , < ) and * ~ _. M2(zl,zz,t\a,s) — 0 (4.33) for an absorbing barrier, or BM2(21’ZZ,C\}’,S) ______..._——————\ _N = 0 (4.34) By Y" for a reflecting barrier at y = 5. Since the initial condition does not divide the spatial domain into different regions, the con- tinuity conditions for the points 21 and ZZ’ Which follow from (4.19) and (4.20) will always be satisfied, and thus need not be onsidered separately. Knowing Mi(zl,t\y,s), the expected number of entities in ny arbitrary interval (51,221 can be calculated from the dif- ference M:(Ez,t\y,s) — M:(El,t\y,s). Similarly, the second actorial moment for the number of entities in (51,221 is given by (2)((zl,zz]x(21,22],t\y,3) = N(z)((0,22]X(0,22].t\y,8) +N(z)((0:E1]X(O,Zl]ztly)s) _ 214(2)((0,ZI]X(0,ZZ],E\YaS) * ~ ~ 4 ~ ~ = M2(ZZ,ZZ>t\YaS) _ 2M2(Zl,22,tiy,5) it + M2(z1,zl,t\y,s). (4.35) 61 nowing the first two factorial moments, the variance of the popula- ion in any interval A can be calculated by using the relation ar(N(A)) = M(2)(A x A) +-M1(A) — (M1(A))2. Remark 4.5: When the transition probabilities related to irth and death, as well as the transition probability density for i he movement of an entity in D do not depend upon time, i.e., K hen the stochastic population process is stationary in time, the ependence of G(e,t\y,s) on t and s will enter only through he difference (t-s). In view of this, G(e.t\y.s) = G(e.t-s\y,0) _ C(Bat”S\Y) ' imilar relationships for the moments follow immediately from the GF. This makes it possible to define a new time scale othat 5—-—a— T = t-s, 65 ET. The backward diffusion equation for the PGF or the stationary process and the initial as well as boundary con- itions are given by replacing the partial derivative - as by and G(e,t\y,s) by G(e,T\y) in (4.12) thru (4.16). Similarly, e corresponding equations for the first two factorial moment asures are given by replacing ~ 2; by 2:, 7’: 7‘: 7': (z,t\y,s) and M2(zl,zz,t\y,s) by M1(z,T\y) and M2(zl,zz,m\y) as well as Spectively, in (4.23) thru (4.34). .4 The Forward Diffusion Equations The derivation of the forward equation for the generating nctionals is much more complicated than the backward equations cause in the derivation of a forward equation it is necessary to 62 ascribe the changes in the generating functional due to the changes 1 positions of all the entities in the population at time t during small time interval (t,t+6t). On the other hand, only the changes 1 a single ancestor had to be considered while deriving the backward quation. As indicated in Section 4.1 a forward description is in arms of the number of entities in an infinitesimal interval z,z+dz) C D. The forward equation for the moments is thus in terms f densities of the moments. The derivation of the forward diffusion quation therefore depends upon the existence of these densities. n the following derivations the moment generating functional (MGF) ill be used for obtaining the forward equations for the moment ensities (for the cases where these densities exist). The PGF or he characteristic functional can also be used in place of the MGF. Kendall (1949) first derived the forward differential equa- ions for the first moment (mean), variance and covariance densities r the age distribution of a pOpulation, using a generating func- ‘onal. The equation for the first moment density was the same as .1) with i = 1, v = 1, and the birth rate appearing only in a 1 undary condition. The densities were defined as (mean) E(N(dz,t)) = fi1(z,t)dz + o(dz), (variance) Var(N(dz,t)) = E2> - (E(N>>2, = fi2(z,t)dz + o(dz) d (Covariance) Cov(N(dzl,t)N(dz2,t)) = E(N(dzl,t)N(dz2,t)) ‘ E(N(dzlat))E(N(dzzat)) = [7112(z11229t) + o(dzl,dz2) : ere N(dz,t) denotes the number of entities in the interval ,Z+dz) at time t. It can be seen that in the terminology of 63 the theory of stochastic population processes the individual state space E for Kendall's problem consists of the non-negative real line denoting the age of an individual. The variance density is connected with the density of the second moment measure with respect to the diagonal of E X E and the covariance density is related to the density of the second moment measure with respect to E X E itself (cf. Sections 2.3 and 2.4). The method of derivation of the equations in this section will be similar to that of Kendall. Let ¢[9,t21N0,t1] denote the MGF for the population at time t2 given the measure N denoting the distribution of entities 0 in D at time t1. In view of (4.4), when there is only one ancestor at y at time s, ¢[9,5\N0,S] = Eiexp je(2)1(dZ\y)] = 9(y) - (4.36) B It can be seen from the properties of conditional expectations that ¢[e,t+5t\NO,s] = E[exp I5 e(z)N(dz,c+5t)\N0,s)] ll E{E[exp L 9(z)N(dz,t+6t)\N(z,t)]\1\l0,s]. (4.37) D Noting that the integral with respect to the counting measure repre— sents a sum over all the entities in the population, E[exp f~ e(z)N(dz,t+6t)\N(z,t)] = E[ H exp e(z_)\N(z,t)] D t+5t J = n Eiexp 8(2.)\N(Z.t)]. (4.38) t+5t 3 where the product has one term corresponding to each entity at Zj E D at time t+6t, and the second expression on the right follows from the independence of the entities. It is easy to see that 64 n EEeXP 6(2.)\N(z,t)]= H {(1-u(z..t)6t)E[exp(e(2.+6z))] t+6t J t J 3 + u(zj,t)6t}{1- )3 xi(zj,t)5t + z yi(zj,t)5t exp(ie(0))}. (4.39) i=1 '= 1 The quantity E[exp(e(zj+6z))] can be written as E{exp[e(z.+62)]} = j” g(6z; z..t)exp(e(z.+62))d(oz) J D J J = exp(e(zj)){1+6t[e'(zj)r(zj,t) O' (Z-3t) 2 + ————%——-(e"(zj)+e'(zj) )]+ o(st)}. (4.40) The second expression on the right side of (4.40) is obtained by expanding exp(e(zj+62)) in a Taylor series about 23. and using (4.1). The quantities e'(zj) and e"(zj) denote the first two derivatives of 9(2) evaluated at 2 = zj. Substituting (4.39) and (4.40) in (4.38) and simplifying, H EieXP(e(z.)\N(z,t)] ~ H exp(e*(2.)) t+dt J t J 7‘: = expij, e (z)N(dz.t)], (4.41) D where 2 e* = e(z)+ot{r(z.t)e'(z) + Q-éfl we) + e'2(z)] - u(z,t)[1-exp(-e(z))] + E xi(2,t)[exp(ie(0))-l]} . (4.42) i=1 Using (4.41), (4.37) can be written as * oie,t+6tiNO.s 1 = oie ,tiN0,s] - (4.43) Equation (4.43) along with the initial condition (4.36) represents a complete forward description of the population in terms of the MGF. Further simplification of (4.43) is not possible. 65 The first moment of the population in the set A CID is obtained from the MGF by using (2.18) and (2.19): ooige,t|N0,s] 0 = - ——.._____. ée¢( ) at g=0 = - 3— {Eiexp J“, ge(z)N(dz.t))\N ,sm a; D o = .vffi e(z)E[N(dz,t)\N0,s)] = IAM1(dz,t\N0,s), (4.44) when 9(2) = -I(A\2). In the derivation of (4.44) it was required to commute the operator :3 with the expectation and integral operators. This can be justified by using an argument similar to that in Remark 4.4. In the case where the number of entities in a small arbitrary interval dz is finite with probability one, the first moment density with respect to 2 exists and Ml(dz,t\NO,s) can be written as m1(z,t\N0,s)dz, where m1(z,t[N0,s) denotes the density of the first moment measure. When 9(2) in (4.37) is replaced by g9(z), (4.43) changes to o[g9,t+6thO,s) = 05[e**,t|NO,s)] , (4.45) where 2 2 e -(Z) = ge(z)+ot{r(z,t)ge'(z) + Egg g[9"(z)+z;e' (z)] - u(z,t)[l-exp(-g9(z))] + X(Z,t)[eXP(Qe(z)‘l)ll' (4‘46) USing (4-44) and taking the set A to be an arbitrary interval (51.52], (4.45) leads to 2 aml (2, thO ,s)dz at 22 4.47 %j m1(z,t\N0,s)dz , ( ) 1 21 N1% N1 66 * where the operator a; is given by a: = {a aztozeétm _ any) .1 _ Wm)“, (4.48) 52 It should be noted that although the indicator function I((21,22]\z) has jump discontinuities at 21 and 22, and thus the derivatives 9'(2) and 9"(z) in (4.46) cannot be defined at these points of discontinuity, it is possible to define an infinitely differentiable function f(z) which agrees with the indicator function with an arbitrarily small error e > 0, (defined as the integral of the absolute value of the difference between the two functions) over a closed interval [&1,&2] containing (21,22] and is zero outside (91,82) C [&1,&2] (see e.g., Indritz, 1963; p. 254). The function f(z), instead of the indicator function itself, has to be used in the derivation of (4.47). Since the limits of integration 21 and 22 in (4.47) are arbitrary, the integrands must be equal, and hence W = 4* m (z t]N s) (4 49) at z 1 ’ 0, . ' When the interval (21,22] includes one of the boundaries, i.e., z = 0 or a, the following boundary conditions result: a[oz(o,t)m1(o,t1N0,s)] 52 r(09t)m1(09t\NOSS) - % '5 co = g Eiyi(z,t)m1(z,t‘N0,S)dz, (4.50) i=1 and a[02(§:t)m1(§!thoas)] - az _ r('a,t)m1(5,th0,s) ‘ % (4'51) for a reflecting barrier, or 67 m1(5,t\N0,s) = 0 (4.52) for an absorbing barrier. The initial condition is obtained from (4.36) as the Diraceé function m1(z,s|NO,S) = 6(z-y) (4.53) when there is one ancestor at z = y at time s. When there are k ancestors located at yj, j = l,2,...,k, the initial condition will be the sum of k 6-functions k m1(z,s‘N0,s) = .2 6(2 - yj) . (4.54) j=l It can be seen that the forward diffusion equation (4.49) for the first moment density is the same as the diffusion equation (3.4) with an additional term accounting for the death of entities, and has the same form as (1.3) where a: = rR = 0 and the death term D replaced by the corresponding term in (4.48). The deriva— tion of (1.3) has been mainly based upon an analogy between mixing of eddies of a fluid in a process vessel and molecular diffusion (cf. Levenspiel and Smith, 1957). Use of the diffusion equation in population balance models in chemical engineering has been restricted to only the residence time distribution analysis of some process vessels. A major reason for not using the diffusion equation to characterize particulate processes such as sucrose crystallization thus far is perhaps a lack of a full appreciation of the stochastic nature of the diffusion process. The development of (1.1) is based upon purely deterministic considerations, with x denoting the number density of the entities 68 in the Euclidean space E. To account for the fact that an integer- valued population is characterized by a continuous density function x, it is always assumed that the population must be large in order to justify the "continuum approximation" (Randolph and Larson, 1971; p. 13). Thus, it is assumed that (1.1) is not valid for small populations. Moreover, (1.1) cannot be used to describe any popula- tion where the movement of all the entities at a given point in E at time t isznot identical in all respects. When the diffusion coefficient %_ is zero, (4.49) has the same form as (1.1), and thus (4.49) represents a generalization of (1.1). In view of this, (1.1) can be seen to be rigorously true even for small populations if x is interpreted as the first moment density instead of a number density and other parameters are given the appropriate probabilistic interpretations, and if the initial population is char- acterized as a sum of Dirac—5 functions instead of a smooth density. Evidently, the assumption of independence of entities underlying the stochastic formulation implies that the resulting partial dif— ferential equations and the boundary conditions cannot have any non- linearities with respect to the first moment density. When the initial population is large, the assumption of a smooth density function may be justified. It is only in this sense that the initial population must be large in order to justify the "continuum approxima— tion" in (1.1). Derivation of the forward equation for the second moment density is complicated by the fact that the second moment distribution has a "mass concentration" along the diagonal of the domain 5 X 5 (See Section 2.4). The second moment of the number of entities in 69 he set Al x A2 c E x 5 can be calculated using (2.18) and (2.19):' 2 2 59¢(0) = :65 [E exp ffige(z)N(dz,t\NO,s)]\g=O = Effie(zl)N(dzl,t\N0,s)ffie(z2)u(dzz,tiuo,s) = EI ~ ~9(zl)e(22)N2(dz1 X d22,t|NO,s) (4.55) DXD y Fubini's theorem on summations with respect to product measures see e.g., Feller, 1966; p. 120). As discussed in Sections 2.3 and .4, the second product moment can be expressed as a sum of a density 2(z,thO,s) with respect to the diagonal of D X D (i.e., the ine 21 = 22) and a density with respect to D X 5. Use of these ensities in (4.55) yields Ejfixfi9(zl)9(22)N2(dz1 X d22,t|N0,s) 2 = 9(2 ) m (z ,t‘N ,s)dz N N 1 2 1 0 l {zlzzzlniDXD} + l I e(z1)9(22)m12(zl,z2,c1NO,s)dzldz2 . (4.56) DXD . a in the case of the equation for the first moment density, taking z) to be the infinitely differentiable function approximating A1\2)I(A2\z) (or equivalently, taking 9(21) and 9(22) to be 1e infinitely differentiable functions approximating —I(A1|zl) 1d -I(A2l22) respectively), the second moment of the population 1 the set A1 X A2 can be obtained using (4.45), (4.55) and (4.56). us, for the density mlza gm (2 ,z ,t‘N ,s) 12 l 2 O _ 7‘: 71: i iA at dzldzz _ i £A {a21+a22}m12(z1’22’tiNo’S)dzldz2’ l 2 l 2 (4.57) * * where a; and a; denote the operator defined in (4.48). An 1 2 additional condition that m1(21,t\N0,s) = m2(21,t\N0,s) (4.58) needs to be satisfied in the course of the derivation to obtain the equation for the density m2: am2(zl,t\N0,s) i21=zzlniA1XA2l at =r‘ {49:11] {2 =2 }fl{A XA } z1 2(Z1»t\NO,S) + N(Z1,t)[m1(zl,t\NO,s) l 2 1 2 - m2(zl,t\N0,s)]}dz1 . (4.59) As mentioned in Section 2,4, the second factorial moment of the product counting measure N2(A X A) does not have a concentration along the diagonal of D X D only when the point process describing the population is orderly. Since the second product moment equals the sum of the first moment and the second factorial moment, it can be seen that (4.58) is satisfied only for an orderly process. When this is the case, m12(zl,22,t\NO,s) completely describes the second factorial moment of the population. In View of this discussion and the fact that A and A2 are arbitrary, the forward diffusion 1 equation for the second factorial moment density is om12(zl,z2.t\N0,s) at 71‘ air = {afl + agz} m12(zl,zz.t|NO.S) (4.60) 71 for the case where the point process describing the pOpulation at any given time is orderly. When A1 or A2 include the boundaries, the following boundary conditions result: 2 BEG (0,t)m12(0,22,tiN 35)] W2 (Oazzat) E r(03t)“112(03223t\N035) "' $2 J; l 621 a m =£ E lxi(zlat)m12(zlazz)t\No)S)dzl (4°61) i=1 (D + __r_. iii(zz.t)m1(z2.t\NO.s), 1-1 '5 oo v22(21,0,t) = g iElixi(22,t)m12(21,22,t\N0,s)d22 + E i)i(zl,t)m1(zl,t\NO,s) (4.62) i=1 and (,2 (5,22,t) = (,2 (z1,a,t) = O (4.63) 1 2 for a reflecting barrier, or m12(5,22,t\N0,s) = m12(zl,a,t\NO,s) = O (4.64) for an absorbing barrier. In view of (4.58) it is no longer nec- essary to solve for m2(z,t|NO,s) separately. When there is only )ne ancestor at y at time s, the initial condition will be m12(zl,22,siN0,s) = 0 . (4.65) if at time 3 there are k ancestors located at distinct points .j’ j = l,...,k, the initial condition will be 72 N .- 9,- x v'c‘k * ) n12(zl.zz.s\ 0.8) - izj6(zl-yi)5(zz-yj) + 12j6(zl -yi)5(22 ‘yj ’ irj i#j (4.66) where * * Z1,Z2 E [0:21] H [0,22] and * ~ * 21 6 {[0,21] U [0,22] ' [0,21] 0 [0:22]} (cf. Figure 2.1b). lariance of the number of entities in any interval A can be readily calculated by using the relation 7ar(N(A)) = £X£m12(21,22,th0,s)dzld22 + £m1(z,t\NO,s)dz 2 - [£m1(z,t\NO,s)dz} . (4.67) 1.5 Comparison of the Backward and Forward Equations It is not possible to obtain an explicit forward partial dif- ferential equation for a generating functional, whereas the backward equation can be readily derived. Although an analytical solution of :his equation cannot be obtained, the extinction probability as rell as the probability generating function for the number of entities -n any set A CiD can be calculated numerically. This possibility .as broad implications in the management and control of biological iopulations. For example, in the biological control of insect pests, ; small number of parasites or predators are usually released in an .nfested area. The extinction probability for the controlling popula- ion can serve as a measure of the failure probability of the control trategy. In case of a finite Spatial domain D and constant parameters, he forward diffusion equation can be solved by using separation of 73 ariables to yield an infinite series, as illustrated in the next chapter. The coefficients in the infinite series have to be valuated by using the initial condition. The initial condition is sually in the form of a sum of Dirac-6 functions, and thus cumber- ome to handle in general -- when the parameters are constant, as ell as in the general case when a numerical solution has to be attempted. In case of a numerical solution the integral in the boundary condition increases the computational effort considerably. As in the case of the corresponding backward equation, the solution of the forward equation with more than one ancestor can also be expressed as a sum of the required number of solutions, each with one ancestor in an appropriate location corresponding to the initial distribution, but this does not simplify the solution to any appreciable extent. When the initial population is large, the assumption of an initial density may be justified and the solution would be much simpler. On the other hand, the backward equation is nuch simpler to handle when the initial pOpulation is small. Although a general analytical solution of the backward equation would be quite difficult, a numerical solution is relatively easier than the forward equation because of the simpler boundary conditions and the fact that the initial condition is a step function rather than a Dirac-5 function. It should be noted that the solution of the back- ward equation for the first moment gives the value of the expected number of entities in [0,2] for a fixed 2 E D at time t, as a function of the initial location of the ancestor and the initial time s. In order to obtain the complete distribution of the expected ~ iumber of entities in D it is necessary to solve the diffusion 74 quation for a number of values of 2. In addition, when there are ore than one ancestor at the initial time s, an appropriate number f solutions each with a single ancestor in a location corresponding 0 the initial distribution have to be added to obtain the final nswer. On the other hand, the forward equation gives the entire istribution of the expected number of entities in the population t time t directly. Thus, when the initial population is large nough to justify a smooth initial density, the forward equation pould be computationally more efficient than the backward equation. The forward equation (4.60) for the second factorial moment iensity is a linear second order partial differential equation with :hree independent variables. In spite of its linearity an analytical ;olution is not possible even in case of constant parameters because )f the coupling with the first moment in the boundary conditions. i numerical solution will be cumbersomebecause of the three independent variables and the integrals in the boundary conditions. Moreover, :0 get the variance of the population in any interval, the densities lave to be integrated over a two dimensional domain (cf. Equation 14.67)). 0n the other hand, the backward equations for both the moments have only two independent variables and the moments of the :opulation in any interval can be readily computed from two (or three) ‘olutions corresponding to the boundaries of the interval, as dis— ussed earlier in this chapter. Thus, the total computational effort ill be much less than that for the forward equations when the nitial population is small. Furthermore, the backward equation is alid even if the population is not orderly and does not need the 2 ontinuity of r(z,t) and o (z,t) as well as the continuity of 75 the first derivative of 02(z,t). Thus it can be seen that the backward equation is a more convenient tool for obtaining the moments of the population when the initial population is small. When the initial population is large, the forward equation for the first moment is often much easier to handle and since the second moment will always be rather insignificant in this case, its evaluation will no longer be necessary. CHAPTER V SOLUTION OF THE DIFFUSION EQUATIONS Solutions of the diffusion equations for the first two moments and the PGF for a population with a one-dimensional individual state Space E are presented in this chapter. 5.1 General Considerations As in the last chapter, the spatial domain 5 is taken as the interval [0,5]. Analytical solution of the equations for the moments is possible when the parameters are constant. For the analytical treatment, values of the parameters are taken as follows: For the entities at z E U at time t, 02(Z,t) = C72 > 0: r(z,t) = r9 u(z,t) = u, (5-1) 0 for 0 5 z < b Ki(z)t) = xi for b s z 5 a, i = 1,2,...,m. The parameters as given by (5.1) are quite realistic for a broad class of problems. In the crystallization process with z denoting the characteristic length of a crystal, the assumption of a constant growth is very common when the supersaturation of the magma is relatively constant. A constant diffusion coefficient may also be justified under a steady flow of magma into and out 76 .7 77 of the crystallizer and constant supersaturation. A constant death rate reflects a perfect mixing in the crystallizer and a constant efflux of the magma, and the assumption of a zero birth rate for 0 s z < 6 implies that a nucleus must grow to a certain size before it can breed new nuclei by secondary nucleation. In biological populations with 2 representing the weight of an individual, the growth rate as measured in terms of the rate of increase of biomass of an individual is high when the individual is young and reduces progressively as the biomass increases. How- ever, it can be seen that by choosing z to be a proper maturity variable the dependence of the mean growth rate r(z,t) on 2 can be eliminated. In cold-blooded organisms and plants the metabolic activity is a function of the body temperature of the individual. Due to the diurnal and annual fluctuations in the atmospheric temperature and solar radiation the growth rates of such individuals in a general ecological system are functions of time. In all poikilothermic species (i.e., cold-blooded organisms and plants) the temperature-growth relation shows a characteristic behavior, which can be used to define a physiological time scale from the chronological time-temperature relationship, so as to have a constant mean growth rate with respect to the physiological time. A constant death rate for individuals in a biological population may be a reasonable first approximation in many cases. Similarly, the assumption of zero re- production rate until a certain value b of 2 (corresponding to the "age” of puberty) and a constant rate thereafter can also be seen to be a reasonable assumption for many populations. The re— lationship of physiological and physical activities Such as repro- duction and locomotion with temperature also shows a pattern similar 78 to the temperature—growth rate relationship. Hence, it may be assumed that the temperature dependence (and thus, to a large extent, the dependence on chronological time) of these parameters for poikilo- thermic species can be removed by the use of physiological time. A detailed discussion of the physiological time is given in Appendix A. ‘Use Of a constant diffusion coefficient in flow through packed beds is very common. Little work has been done so far in the applica- tion of the diffusion equations in the modeling of biological pOpula- tions. A constant diffusion coefficient is usually used in the spatial migration studies. Stuart and Merkle (1965) also used a constant diffusion coefficient in their study of cell dynamics. However, a careful study of the relationship of the diffusion co- efficient with temperature and the state 2 of an individual must be made to check whether the same transformations in z and t to yield new maturity and physiological time variables lead to con- stant values for all parameters including the diffusion coefficient. In the case where the diffusion coefficient is zero and the other parameters are given by (5.1), it can be seen that a simple transformation on 2 by defining E = z/r reduces the rean rate of change r to unity in terms of E. In this case the equations for the moment distributions are the same as those for the age dis- tribution of a biological population, and the results on the age- dependent birth and death processes can be used directly (see e.g., Kendall, 1949; Goodman, 1967; Bartlett, 1969). In the solutions presented in this chapter it is assumed that the diffusion co— efficient is always positive. 79 It should be noted that (5.1) implies stationarity of the diffusion process (i.e., that the parameters are not time—dependent) and the abbreviated notation defined in Remark 4.5 can be used. For abbreviation, the first moment density m1(z,t\N0,s) is denoted by m1(z,T), where T = t—s. The initial distribution N0 of the entities in D is expressed as a Dirac-6 function for a single ancestor. It may be possible to characterize N by a smooth 0 density function for a large initial population. Also, the m (D * quantities 2 iii and Z i(i-l)xj are abbreviated as x and i=1 i=1 ' ** 1 respectively. 5.2 Solution of the Diffusion Equations for the First Moment The forward equation is used to obtain a solution for the case where D is finite. For the case of a semi-infinite spatial domain the backward equation is used to obtain an explicit analytical solution when E = 0 and the process starts with a single ancestor at y = 0. Case 1. Finite Spatial Domain [0,3] wigh a Reflecting Barrier a; a. The forward diffusion equation and the corresponding boundary and initial conditions for this case are 2 02 a m1(z,'r) am1(z,'r) am1(z,w) 2 2 ‘ 1’ az - um1(z.'r) - —-—-—~aT , (5.2) AZ ~ 02 am1(a97) rm1(a,7) - 5‘ ~——g;*—— = 0, (S-Za) 2 am (On) '5 1 7': rm1(0,¢) - g— T = x [m1(z,'r)dz (5.2M E and 80 m1(z,0) = m0(z) . (5.2c) Let m1(z,T) be expressed as a product m1(z,T) = Z(z)T(T), (5.3) where 2(2) and T(T) are functions only of z and T respectively. Substitution of (5.3) in (5.2), (5.2a) and (5.2b) yields 2 dZZ dZ g_ (z! _ r {Z} _ 01(2) = (12(2), (5.4) 2 2 dz dz 2 dz 3 rZ(5) -Q——D—=o, (5.4a) 2 dz 2 d2 0 * E ‘ rZ(0) - g— —L)~ = x J‘Z(z)dz, (5.1m) : dz ~ b and 31—21% at, (5.5) r where a is a constant. When a > -(“—E + u), the general solution 2 0' of (5.4) is given by 2(2) = A1 cosh 32 + A2 sinh Bz, (5.6) where ”47' + 202m + 0!) B ‘ 2 ’ 0 (5-7) 2 and A1, A2 are constants. Similarly, when a $_(£_§ + u), the 20 general solution of (5.4) can be written as 2(2) = A3 cos E2 + AA sin 62, (5-8) where 2 — - 2 B = r 2 L + , (5.8a) 0' and A3, A4 are constants. Using boundary condition (5.4a) it is possible to eliminate A1 and A3 from (5.6) and (5.8) respectively. Thus, (5.6) yields f1'=‘(z)=A2[(gig—41%)COShcosh+sinh1 0' Botanhmoa) "r 2knz (2knz w W2 + kzl a ”+1exp(ak+1T)[ cos s( )+sm a)]}exp(rz/oz ), (5-17) where 80 = .jrz + 2021* . (5.18) Since the eigenfunctions are not orthogonal to each other, Remark 5.1 applies in this case as well. Remark 5.2: Note that in the infinite series solution of the diffusion equation it has been tacitly assumed that the eigen- functions form a complete set (see e.g., Indritz, 1963) and that all the eigenvalues are distinct. Also, the statement that asymptotically the expected number of entities in any set A C D increases exponentially with time is true only when the dominant eigenvalue (i.e., the eigenvalue with the largest real part) is real and distinct. It would be quite difficult to prove the complete~ ness of the set of eigenfunctions and distinctness of the eigen- values in general. Nevertheless, it is easy to see that for the case of a reflecting barrier at a with b = 0 discussed above, the eigenvalues are all real and distinct, and the eigenfunctions form a complete set over a class of functions continuous over a finite interval. Although the hypothesis that the dominant eigen— value is real and distinct may be very difficult to prove analytically, it is easy to show that there must always be at least one real dominant eigenvalue. The proof of this is as follows: If all the 85 dominant eigenvalues are distinct and have the form alR + ialj’ j = 1,---,k; alj > 0 and i = v-l, then the asymptotic solution will k be of the form exp(a1RT) Z Ajexp(ia1jT), which will be an oscillating i=1 function of time assuming negative values during certain time intervals. The same is true even for the case where the dominant eigenvalues are all complex, but not necessarily distinct. Since the first moment distribution must be a nonnegative quantity, at least one of the alj's must be zero, and the corresponding co- efficient Aj large enough as compared to the other coefficients so as to compensate for the negative contributions from all the terms involving the complex eigenvalues. Case 2. Finite Spatial Domain [0,5] with a3 Absorbing Barrier at 5. The forward diffusion equation, its initial condition and the boundary condition at z = 0 will be the same as in Case 1. The boundary condition (5.2a) will be replaced by m1(5,T) = 0 - (5.2d) Equations (5.3), (5.4), (5.4b) and (5.5) also apply in this case. Equation (5.4a) will be replaced by 2(5) = o . (5.4c) 2 When a > -[-—— + u], the following expression is obtained by apply- ing (5. 4c) to 0the general solution (5. 6): Alcoshsa + Azsinhfia = 0 Hence, f1(z) E Z(z) = A2[sinh(ez) - tanh(Ba)cosh(Bz)]exp(rz/02). (5.19) 86 2 Similarly, when a S -(£—§'+ u): 20 2(2) = A4(sin(fiz) - tan('65)cos(Ez))eXp(rz/Oz)- (5'20) When the boundary condition (5.4b) is applied to (5.19), it is found that the dominant eigenvalue is the root of (adu)[£78inhBa+BcoshBZ}-A*exp(r5/02){Easinh[a(a—S)] o o N 1 + Bcosh[a(§-b)]} + Bicexp(ra/gz) = 0, (5.21) with the largest real part. If the dominant eigenvalue is unique and real, it can be seen that as in the case of a reflecting barrier, asymptotically the expected number of entities in any set A C D varies exponentially with parameter a1. When the boundary condition (5.4b) is applied to (5.20), the following expression for the nondominant eigenvalues results: 1’ -~+ ~ * S/Z-r—‘[”~5 (0+M){-§ SlnBa B cosBa} - l exp(r c )[ 2 Sin 8(a— )] O' O' + g cos[B(5-b)]} + [*3 exp(ra/cz) = O . (5.22) It should be noted that Remark 5.2 is applicable in this case also. Similarly, in view of the fact that the eigenfunctions are mutually independent but not orthogonal, Remark 5.1 is also applicable. An explicit analytical solution for the eigenvalues is not possible even when S = 0. Case 3. Semi-Infinite Spatial Domain [0,m). When 5 is finite and a a m with an absorbing or reflecting barrier at 5, proceeding exactly as in Case 1 and Case 2, it can be shown that regardless of the initial condition, the asymptotic 87 value of the first moment density is m1(z A) ~ f1(Z) exp (alw) , where a1 is the root of the equation J2+2Q2g+2 r - r N _ 2 a )b} — o O 01 + u. - x*exp{( (5.2151) with the largest real part, f1(z) = alexp(pz), and a1 is a constant. In this case it is also possible to obtain an explicit analytical solution for the first moment distribution when E = 0 and the ancestor is initially at zero (Brockwell, 1972). The backward diffusion equation and the corresponding initial and boundary conditions for the first moment distribution are ‘k BM1(Z My) 7': 9c 7% = + , 0 , (5-23) 5T ayM1(z,T‘y) x M1(z T\ ) where the backward diffusion operator a; is given by (4.11), l for y s z 9 M;(z,0|y) = (5.23a) 0 for y > z, 7': W = (5.2%) ay fr _ * 5.23c) M1(z,rr\z+) — M1(z,T\z_), ( * 8M (Zn y) * 1““‘L‘\ = aEL-g-L'fldfl-l _ (5.23d) ay y-Z_ ay y=z+ M:(z,¢1m) = 0 . (5.23e) Taking Laplace transforms of (5.23) thru (5.23e) the partial differential equation can be transformed into the following ordinary differential equations and boundary conditions: 2— _ 02 d M1(Z’P\Y) dM1(z,ply) _ 5‘ 2 + r"“‘g§““‘ ' (M+P)M1(Z,P\Y) dy *_ l for y S 2 +‘k M1(Z,P\0) ={ (5-24) 0 for y > z, dfi __l_______ = 5.24 dy y=0 0’ ( 3) fi1(z,p\z+> = fi1, (5.241» _ $3 , dM1(z,P\y) = 1(2 p\y)\ (5.24C) dy y=z+ dy y=z_ fi1(z,p\m) = 0 , (5.24d) where p denotes the Laplace transform variable and M1(z,p\y) * n is the Laplace transform of M1(z,T\y). The general solution of (5.24) is {3‘4 (z,p\0)+1 _ ----— + alexp+a2exp u+P + a3exP(Bly) + aAEXP(BZY): y > z; (5.25b) where are constants to be determined from the boundary and a1,a2,a3,a4 conditions. In addition to these, fil(z,p\0) is to be evaluated by setting y = 0 in (5.25a): *— )\ M1(Z,P\O) + 1 5.27 2 u + p ( ) M1 = a1 + a In view of (5.24c), a3 = 0. Applying other boundary conditions, one obtains aaexp(32y) = Jig + alexp(31y) + azexp(32y), (5.28) 82a4eXp(Bzy) = BlaleXp(Bly) + eZaZeXP(Bzy) (5.29) and (5.30) B131 + Bzaz = 0 Solution of (5.27) thru (5.30) for a1,a2,a4 and M1(z,p\0) yields ._ exp(‘Blz) Ml(z,p\0) = “¥;“ ' _‘“;_‘-- p-x +u P‘A +H 2 l2 2 2 1 exp(—rz/O )exp([- r +20 (LL+p)lV/Gl_ , (5.31) = —:— - 7‘: p-x9+u P ' A + M Taking inverse transform (cf. Cox and Miller, 1965; p. 221, Equations 73 and 75), the analytical solution 9O * T _ * _ _ 2 2 M1(z,T\O) = eXp[(k*-u)T]{1 _ i g eXPE X u (z ru) /2uo ] du} (5.32) V 21’1'u3 is obtained. 5.3 Solution of the Diffusion Equations for the Second Moment As mentioned in the previous chapter, it is not possible to obtain a complete analytical solution of the forward or backward equation for the second moment distribution. However, an asymptotic solution (i.e., a solution for large values of time T) of the back— ward equation can be obtained by using Laplace transformations. In this section asymptotic solutions for the second factorial moment distribution will be obtained for the case where there is only one ancestor initially. These can be used to evaluate the second factorial moment distribution for a population with k ancestors initially at y1,...,yk, from the relation * k * M2(zl,z2,m\y1,...,yk) = _21M2(zl,z2,T\yi) l: k k 'c 3? + z M’(z ,T\y.) 2 M (z ,¢\y.) .- (5.33) . l 1 1 . l 2 J l=l J=1 j#i Equation (5.33) readily follows from (4.6) and (4.20). §§§g_l. Finite Spatial Domain [0,5] with a Reflecting Barrier g; a. The backward diffusion equation and the corresponding boundary conditions are given by * 5M2(zl’22’lly) B’T * ~ = f 0 s < b, (5-348) a§M2(z1.z2,T\y) or y * BM2(21’22’T\Y) = M* ** * * 3T a; 2(z1,z2,7\y) + h M1(zlaT\O)M1(Zz,T‘O) + *M* + k * k x 2(21,22,T\0) x [M1(21,T\0)M1(z2,lly) * * + M1(z2,T\O)Ml(z1,T\y)], for 5 s y g a (5.34b) where a; is the backward diffusion operator defined in (4.11); * BM2(21,ZZ,Y) —'—'—'_——T = 0, 5.34 by y=0 ( c) * ~ * ~ 5 3 M2(zl,z2,T\b_) — M2(z1,z2,T‘b+), ( - 4d) M* * a 2(Zl:ZZ:T\Y)\ ~ 2 6M2(zl’z29T\y) 5y y=b_ 5y y=b+ (5.346) and My? \ ) a z .z ,T y 2 1 2 (5.34f) ay ly=5 * ** . . Since A and A are step functions defined from xi(z,¢) in (5-1) with a jump discontinuity at y = b, the backward diffusion equation has to be split into two parts (5.34a) and (5.34b). The conditions (5.34d) and (5.34e) follow from the fact that the second moment distribution and its derivative with respect to y must be continuous on [0,5]. Taking Laplace transforms of (5-343) thru (5.34f) with respect to time, the following set of equations is obtained: 2 dzfi (z z p‘y) dfi (z ,Z ,p\y) o 2 1’ 2’ + 2 1 2 §~ -——————-——————— r ————————-——-—- 2 dy _ _ , . = 0 (n+p)M2(zl,72 P\y) for o s y < E, (5.35a) 92 2_ _ 02 d M2(z1,z2,p\y) + dM2(z1,z2,p\y) _——————_ r _— 2 2 - (”+p)b—’1 (Z ’2 :p Y) dy dy 2 1 2 \ *7? 7’:— + A 21(zl,zzip\0) + k M2(zlizz,p\0) * N + i [@12(21,Z2,p\0,y) + o21(z1,zz,p\y,0)] for b s y S a, (5.35b) where 7'c 7k §1(zl,zzip\0) = 1%[M1(21,T\0)M1(ZZ,T‘O)], (5-368) 7': 7\- @12(Zl!22:p\03y) = £T[Ml(Z1)T\O)M1 (ZZ,T\Y)15 (5'36b) 7k 7k @21(Z1)223P\Y’0) = £T[Ml(zl)T\y)M-l (22,T\0)] (5.36C) and 1%[f(z,T\y)] = Laplace transform of [f(z,¢\y)] with respect to T, 00 = gexp(—pt)f(z,t\y)dt . (5.37) The boundary conditions corresponding to (5.35a) and (5.35b) are dM2(Zl’Zz’p\Y) \ = 0, (5.35c) dy y=0 * N = ‘ ~ 5.35d) M2(Zl)zz)p\b_) M2(zl)22,p\b+)3 ( dM2(Zl’ZZ,p\Y) N = dM2(Zl>ZZ:P‘Y)\ N (5.35e) dy y=b_ dy y=b+ and db? ,. __2£il_i§_3l31l N = 0 , (5.35f) dy y=a The general solutions of (5.353) and (5.35b) are _ B 5.38s M2(zl,z2,p\y) = alexp(Bly) + azexp(BZy) for O S y < ( ) 93 P71é(z1,zzip\y) = a3eXP(Bly) + aAeXp(BZy) *— ** x M2(zl,zz,p\0) + x 21(zl.zz,p\0) + p + u .1" 2K {@1(Zlyzzap\y,0)eXP(Bly) — ®2(zl,z2,p\y,0)exp(32y)} + 2 0 (Bl _ 62) for b s y S a (5.38b) respectively, where ®1(zl,z2,p\y,0) = fexP(-81Y)[élz(zl,zz,p\0,y) + Q21(212223p‘y30)]dy (5.38c) @2(Z1,Z2:P\y50) = jeXP(-Bzy)[@12(z1,z2,pl0,y) + 221(zl.zz,p\y,0)]dy and 81,52 are given by (5.26). The constants a1,a2,33,a4 are to be evaluated by using the boundary conditions (5.35c) thru (5.35f). In addition, fi2(z1,z2,p\0) is evaluated by setting y = 0 in (5.383.), i.e., V H — a "i‘ a , 5.35 ) As discussed in the previous section, the asymptotic solution is independent of the initial location of the ancestor. It is there— fore sufficient to obtain an explicit expression for fi2(zl,zz,p\0) in order to evaluate the asymptotic solution. Applying the con- ditions (5.35c) thru (5.35g), the following expression is obtained: 94 M2(zl’z2’P\0) = Ph1 qt h1F(21)F(zz)ex9(2“lT)} (5.49) 7': ~k~ 7': 7': ~ 71‘ (n+2a1)[exp(61§)-exp(Bza)]—k [exp(51[a-S])-exp(32[§—b])] * a where Bl and 52 denote 81 and 52 respectively (see Equation (5.26)), evaluated at p = 2a1. The asymptotic value of the second factorial moment of the number of entities in any given set A 55(21,22] C.D can be easily obtained from (5.49) by using (4.35). Thus, 97 M(2)(Z X K.T\y) (k* M+2l )[F1(zz)-F1 (21)] 2E6X:(Bl(3-5))-EXP(BZ(3-5))]eXp(2a1T) (5.50) (2a1+u)[exp(81a)-exp(Bza)]- x *[exp(81(a-b)) -exp(52(a-b))] Similarly, the asymptotic value of the expected number of entities in the set A E [21,22) is given by M1(A,T\y) ~ [F1(52) ‘ F1(§1)]€XP(a1T) - (5.51) The variance of the number of entities in the set A can be readily obtained from the two factorial moments (5.50) and (5.51). The co- efficient of variation CV(A) of the number of entities in A gives a measure of the stochastic fluctuations in the population N(A) of the set A, and is defined as cv<7i> = [Var)f/E(N>, where Var(N(A)) is the variance of the number of entities N(A). Thus, using (5.50) and (5.51) the asymptotic value of the coefficient of variation can be seen to be cv(Z) , [<2a1+u)iexp(ega)-exp{exp(3 :[a—E]>-exp(e:ifi 43m}15 .(5. 52) l$2a1+u)[eXP(Bla)-exp(Bza)] x*{exp<31[a-b])-exp(82[3- b])} When the initial number of ancestors is k, located at y1,y2,...,yl; ( from (5.33) one obtains (A x A,T\y1,...,yk) ~ k M(2)(K x A,wiy) + k(k—1)[M1(A,m‘y)]2 (5,53) “(2) 98 where M(2)(A X A,mly) and M1(A,m‘y) are given by the expressions on the right side of (5.50) and (5.51) respectively. Since M1(A,T\y1,Y2:~--,yk) N kM1(A:°°l}7): (5'54) it follows from (5.53) and (5.54) that the asymptotic value of the variance of the population in any set A C D with k ancestors is k times the correSponding value of the variance with a single ancestor. Hence, when there are k ancestors initially, CV(A) 7"... 7‘1"” *7? 7% 7* N ~ * N ~ 13 {(zal'ili)[exl>(823)-exp(Bla)]+(l +31 ){6XP(Bl[a-b])~6XP(82[a-b])} __’ N 7'6 in at it N ~ at N ' k(2a1+u)[eXp(815)-exp(82a)]-x {exp(Blfa-b1)-eXP(62[3-b])} (5_55) When 5 = 0, (5.55) assumes a particularly simple form: *9: at + lzt CV(A) ~ A~——i;l———Ji} . (5.56) k(x - H) Case 2. Finite Spatial Domain [0,3] with 22 Absorbing Barrier a_ a. / The evaluation of the second factorial moment distribution I as well as the variance and the coefficient of variation of the number of entities in any given set A C.D proceeds exactly as in Case 1. The diffusion equation and the boundary conditions are the same except (5.34f) and (5.35f), which must be replaced by 4 M2(z1,z2,T‘a) = O (5.34g) fi,(z1.z2.p\a) = 0 (5.35b) respectively. The asymptotic solution for M(2)(A X A,T|y), where 99 A E (21,521 for this case is given by M(2)(Z X K.le) “' N ~ 2 ** * x x N ~ [F1<22>-F1<21>] exp(zle){(K +2i )[Bzexp[81(a-b)] * * N ~ * * - elexp[82(a-b)] ' 82 + 51]} (5.57) * *~ * *~ * a a N N (Zdlm)[82exp(61a) _Blexp(62a)] -A {B2exp[el(a_b)] * * ~ * - Blexp[82(5-b)]} + i*[ogexp-elexp-BZ+BIJ 2 — (241+n>[e:expl CV(A) ~ (5.58) * *~ * *~ * w * ~ ki(201+u)[BzeXP(Bla)'Blexp(823)]'K [Bzexp[61(a-5)] * k * 4 * - elexpieze—Em + A (62 - 31>} Case 3. Semi-Infinite Domain [0,w). Proceeding exactly as in Case 1 or Case 2, the following * expression is obtained for the Laplace transform of M2(z1,z2,T|0): 7': 2* x N x 21(zl,zzip\0) ;%gr ®1(zl.z2.plb,0)eXp(815)(n+P) 2 . (5.59) i (z ,2 ,P‘O) = . 2 l 2 ~ x (u + p)exp(Blb) - x Noting that a1 is the dominant eigenvalue given by the root of (5.21a), as before, F1(zl)Fl(zz) P - 2&1 §1(Z1,Z2,T‘O) N 100 Also, noting that Q ~ = U1(21’22’Plb’0) 1f£§12(21’z22P\0aY)+¢21(21,22,Ply50)]exp(‘51Y)dY}y;B where 51 2 0, so that exp(-Bly) is bounded on [0,w), similar argument as before can be used to show that ‘2F1(21)F1(zz)eXP(-Blb) (P - 2&1)Bl ®1(zl,zz,p\b.0) ~ Hence, for the set A E (51,52] Cifi with one ancestor initially, N N 2 9h? ) [F1-F1(z,>] (i + 2i*)exp<2o1i> (5.60) M(2)(AXZ,T\y) ~ + 2 *N x (n a1)eXp(Blb) - x and when there are k ancestors initially, the coefficient of variation of the population in any set A in D is given by ** * *N g (l + 3x ) - (n+2a1)exp(81b) CV(A,T\y) ~ *~ * . (5.61) k[(u+zal)eXP(Blb) - i 1 At this point it is worthwhile to summarize the important analytical results obtained thus far regarding the solutions of the diffusion equations for the first and second factorial moment dis- tributions with parameters given by (5.1) (i) A complete analytical solution (5.32) has been obtained for the first moment distribution for the case where 5 = m b = 0 and at time 0 there is one ancestor at y = 0. (ii) General solutions for the first moment density in the form of generalized Fourier series of eigenfunctions (assuming that the eigenfunctions form a complete set) 101 have been obtained for the case of a finite interval [0,3] with absorbing and reflecting barrier at 5. (iii) It has been shown that asymptotically the expected number of entities in any set A<: D changes exponentially with a parameter 0 (cf. (5.16)). 1 Here the parameter a1 is given by the dominant root of (5.11), (5.21) and (5.21a) for the cases of a reflecting and an absorbing barrier at a finite 5 and a semi—infinite spatial domain [0,m) respectively. Similarly, the asymptotic change in the second factorial moment of the number of entities in any subset of D is exponential with respect to time, with parameter 2&1 (cf. (5.50), (5.57) and (5.60)). (iv) The asymptotic behavior of the coefficient of variation with a given number k of ancestors at time 0, as given by (5.55), (5.58) and (5.61) constitutes a result of particular importance regarding the limiting behavior of the population: The coefficient of variation of the number of entities in any set A C.D attains the same (constant) value depending upon the boundary conditions and the initial number of ancestors. Furthermore, the asymptotic value of the coefficient of variation is inversely proportional to the square root of the initial number of ancestors. Of course, as indicated before, these results are based upon the hypothesis that al is real and distinct. 102 In view of the result regarding asymptotic behavior of the coefficient of variation it can be said that if the initial popula- tion is Small, the stochastic fluctuations in the population of any subset of D will always be significant even if eventually the population becomes very large. On the other hand, when the initial population is large, the stochastic fluctuations will be small and the equations for the first moment alone may be adequate to describe the dynamics of the population. For example, if in the study of the secondary nucleation process in a crystallizer the experiment is started with a small number of crystals, the number of crystals in any size range will always tend to exhibit relatively large fluctuations. This will lead to an appreciable amount of scatter in the data on the induction periods as observed by Kane (1971). Similarly, the results indicate that a biological population started with a small number of ancestors will always tend to exhibit relatively large stochastic fluctuations. 5.4 The Problem of ”Critical Length" Referring to (5.11) and (5.21) it can be seen that when the dominant eigenvalue a1 is real, for a given set of values of the parameters 02, r, k* and b, there exists a critical value acr of a, such that if a < acr, a1 < 0 and the expected number of entities in the population of any subset of D will always de— crease with time. The problem of determining this critical location of the boundary (the "critical length" of D) has important implica- tions in case of biological populations. For example, in an intensively exploited fishery or a forest resource, few individuals 103 in the population can be expected to survive beyond a certain maturity a. In this context, y and 2 correspond to the maturity of a tree (measured in terms of its height or productive value) or the length or weight of a fish at times 0 and T respectively, l* denotes the expected number of births per unit time, and u refers to the death rate due to natural causes as well as harvesting. Of course, it is assumed here that the growth of an individual can be char— acterized by a diffusion process. In these cases it is reasonable to construct a mathematical model with an absorbing barrier placed at a. As the harvesting pressure increases, the maturity correspond— ing to the absorbing barrier decreases. If the absorbing barrier corresponds to too short a life-span, the recruitment due to repro- duction would be too low to have a self-sustaining population. Evaluation of the maturity corresponding to this critical life— length would be of a great help in determining the necessary manage— ment and control policies such as the legal size-limits on fish so as to avoid extinction of the species. An opposite situation may occur in a continuous crystallizer where the crystals from the product stream are classified and all the crystals below a certain size 5 are returned to the feed stream. In the context of a crystallizer, y and 2 denote the size of a crystal (usually expressed in terms of a characteristic length of the crystal), X* characterizes the expected number of nuclei generated due to secondary nucleation and u refers to the rate of removal of crystals from the processing system. As in case of the biological populations mentioned above, the crystals are assumed to grow according to a diffusion process. For a stable operation of the 104 system the number of crystals in the crystallizer must remain con- stant. To achieve this it may be necessary to regulate the maximum size 5 of the crystals returned as the feedback. Another way of achieving the control would be to fix a and regulate the fraction of the total number of crystals of size less than or equal to s to be returned. Thus, the problem reduces to that of determining the critical size for a fixed ”death rate" or that of evaluating the "death rate" to make a given size a a critical one so as to have a stable operation of the crystallizer. In view of these possible applications, the problem of "critical length” was solved for the case of an absorbing barrier with parameters given by (5.1). In the case of a reflecting barrier at a with b = O, the criticality of the change in the population will be dependent upon X* and M alone. However, when b > O, for given values of A* and N there may be a critical location act of the boundary such that for 5 < ear, the population will (asymptotically) de- crease monotonically with time. The computation of scr in this case will be similar to that given below for an absorbing barrier. Asymptotically, the expected number of entities in the population changes exponentially with a parameter a (i.e., the 1 largest eigenvalue of the diffusion Operator), which is given by the dominant root of (5.21) for the case of an absorbing barrier at s. It can be seen from (5.21) that al is a continuous function of the location 5 of the barrier. The critical length is thus given by setting a = O in (5.21) and solving the reSulting expression for 5, subject to the condition that the root should be greater than or 105 equal to b. A solution act < b would represent a physically impossible situation in view of the fact that no entity in [0,b) is capable of reproduction (cf. (5.1)). Since (5.21) is nonlinear in 5, the solution has to be obtained numerically. Equation (5.21) can also be written as N * (mflelexpmz‘a’) - BZeXMBlaH - x {slexpis2(fi-S)] — Bzexpial(a-E)1} - x*(62'81) = o, (5.62) where _ —r-FV r2 + 2 2 + _ —r - r2 + 2 2 + 81 ’ 2 and E’2 ‘ 2 O 0 Equation (5.62) has to be solved for a after setting a1 = 0 to give acr. If the solution to the problem of critical length is to be obtained in the form of a set of families of curves, the para- meters in the graphical solution will be 02, r, K*, H and b. The total amount of computation required to evaluate the solution over a range of values of all these parameters can be substantially reduced if the diffusion equation for the first moment density is written in terms of the following dimensionless variables: E = z/a ? = rT/fi fi = ZrE/o' B = 5/5 (5.63) § ” Y/a p = ua/r 77 H 7 n3 \. H and an expression analogous to (5.62) is derived in terms of these dimensionless variables. This simplification reduces the problem of evaluation of the critical length to that of solving the equation 106 (mphlfiexmze) - vzf’exp(v1P)] + V2?exp[v1(l3-13B)] — leexp[v2(P-Pb)] + le - vzfi = o (5.64) for E. In (5.64) =-1+114-;/§ V1 2 V —'_1'—M_W 2 2 The computer program used to solve (5.64) is given in Appendix B. It should be noted that the quantities fi/X, v1, v2 and fib are all independent of E. Figure 5.1 gives the solution in the form of a family of curves for f corresponding to the critical value of 5, plotted against fi/fi, with fib and X/fi as parameters. Given the parameters %3 , r, u, l*: and b, the graph can be used to evaluate acr as follows: First, the dimensionless parameters l/fi and P5 are computed and the curve corresponding to these values of the parameters is located on the graph. The value of @/§ is calculated and the corresponding value of f is obtained2 from this curve. Knowing the critical value of § as well as g- and r, scr can be readily computed. As may be expected, the curves indicate that the value of 5 is less sensitive to the changes in the death rate when the CI‘ ratio of birth rate to the death rate X/fi is high or when the value of fb is low, which would occur if b is close to zero or if the mean rate of movement r toward the absorbing barrier is small as compared to the diffusion coefficient. 107 50 I 1' l l T l T I f j -__-__. M; = 2 :' __ a, = 10 _ w = 20 J L... r. I I 1 A | l L L 0.0 0.05 _,1/13__> 0.10 0.15 Figure 5.1 Graphical solution of Equation (5.64) for "Critical Length" 108 Although the main purpose in presenting Figure 5.1 is to demonstrate the method, care was taken to have a realistic range cf values of the parameters. The data of Cooper and Latta (1954) and Cooper, Latta and Schafer (1956) on populations of bluegills Lepomis macrochirus was used to obtain crude estimates of the diffusion and drift coefficients and the death rate. Assuming a constant death rate estimates of the yearly total mortalities were obtained by Cooper et. al., which were used to calculate the death parameter u. The data consisted of the siZe (length) distribution of the various age classes of the fish for several years. The age classes were char- acterized according to the age of the fish in years. Since spawning occurs during a relatively short period of time during each year, it was assumed that the size distribution of the fish in each age group had resulted from a large number of identical young ones born at the same instant of time. If the temporal changes in the expected size distribution of the fish can be characterized by a diffusion equation with constant parameters, the size frequency of the fish grown from a cohort should be approximately normal at all times if the effect of the reflecting boundary corresponding to zero size can be neglected. This may be the case when the drift coefficient is much larger than the diffusion coefficient. Under these conditions the changes in the mean and variance of the size frequency distribution during a time interval At will be rAt and ozAt respectively, and thus the drift coefficient r and the diffusion coefficient %_ can be readily estimated. The assumption that the size distribution was normal at all times was indeed rather crude, and since the purpose 109 here is mainly to demonstrate the technique, no further refinements in the estimates of the parameters were attempted. As discussed by Cooper et. al., horizontal as well as vertical estimates of the parameters were obtained whenever possible. A horizontal estimate is obtained by following the same age cohort year after year, e.g., by considering the two-year old fish in 1954 and the three-year old fish in 1955, assuming that the population estimates for the successive years were obtained with the same accuracy. A vertical estimate uses the different age groups counted during the same season, based on the assumption that the population is in a steady state as regards the yearly recruitment, mortality and growth patterns. The horizontal and vertical estimates were comparable to each other. Different estimates of the parameters are summarized in Table 5.1. It was assumed that the fish were capable of repro- duction after they reached a length of 2.5 inches. The values of birth rates were chosen to be simple mpltiples of the death rates. Based on these estimates the dimensionless parameters for the dif- ferent curves were chosen to include the ranges of parameters re- presented in Table 5.1. Since all the parameters in Figure 5.1 are dimensionless, the same curves can be used to calculate the critical value for any population with any arbitrary units (such as productivity, biomass, etc.) for 5. However, it must be remembered that in many cases the 2 * constant values for o , r, and l are only an approximation to the real situation and due caution should be exercised in using such charts for practical purpose. For example, it can be seen from 2 Table 5.1 that %_ and r are functions of the size of the fish. 110 TABLE 5.1: Model Parameters for the Dynamics of Populations of Bluegills (Lepomis macrochirus) A. Mean Growth Rate r in Inches/Year (Vertical Estimate) Age Class (years) 1 2 3 4 5 6 Remarks 2.5 2.2 1.0 0.8 1.0 0.6 Sugarloaf lake, Michigan 1954 2.2 2.0 1.2 1.0 0.9 0.6 Sugarloaf lake, Michigan 1955 2.5 2.3 1.43 1.35 0.78 0.26 Whitmore lake, Michigan 1955 as 2.4 2.16 1.21 1.05 0.89 0.47 Average of the estimates above 2 2 hr B. Diffusion Coefficient Q /2 in (inches) /yearj First two years Third year Fourth year Remarks i 0.0417 0.0245 0.024 Vertical estimate (1954) i --- 0.0201 --- Horizontal estimate (1954—55) *7’: for Sugarloaf lake, Michigan C. Average Death Rate u in (years)-1 Vertical estimate Horizontal estimate Remarks 1.11 0.655 Sugarloaf lake, Michigan, 1%2,1%4amil%5. 2.2 1.56 Whitmore lake, Michigan, l953~1956. 0.892 --— Fine lake, Michigan, 1955 0.415 -—- Fife lake, Michigan, 1956 111 Similarly, many fish spawn only once in a year giving rise to a "pulse" of young ones. In contrast, the model defines xiéw as the probability of producing i offspring in the time interval (T, +6T) irrespective of whether the individual had reproduced at any time in the interval (—m,¢). 5.5 Solution of the Diffusion Eguation for the PGF -- Computation of Extinction Probabilities As mentioned in Remark 4.2, G(e,T\y) gives the probability that there are no entities in the population, given that there was one ancestor initially at y, if the arbitrary function 9(z) is taken to be identically zero. To obtain this the backward diffusion equation for the PGF was solved numerically using finite differences. As in the previous section, the main purpose of these simulations was to study the nature of the solution for some simple cases and hence constant values were chosen for 9— , r and u from the range 2 of values covered in Table 5.1. Moreover, the probability of having more than one birth in a short time interval 6T was assumed to be negligible, i.e., only Al was assumed to be a significant parameter. Furthermore, k1(y,s) was aSSumed to be independent of time. Equa- tion (4.2) was written in the finite difference form as Elair+0ily)-G(e,ily) = g3 [C(e,Tiy+6y)-ZG§8,Tiy)+G(9,Tiy-6y)] 0T 2 6y + r[G(e,Tiy+6y)-G$9:T, 1 + nil - G(9,T]y)] (5.65) 6y - 11(y)G(e,T\y)[1 - C(e,"r\0)] I 112 for 6y S y s a—fiy and T 2 O. For y = 0 or a, (5.65) has to be modified to include the boundary conditions. For the stability 6 2 of the numerical computation scheme a ratio of g%— = 4 was chosen. Different time and space increments were tried and values of St and 5y were chosen so as to give results agreeing up to the first three decimal places with those for a much finer discretization. The values of the simulation parameters are summarized in Table 5.2. The results of simulations 1 thru 7 are presented in Tables 5.3 thru 5.9 respectively. A computer program used for the simulations is presented in Appendix C. It can be seen that in case of a reflecting barrier at 5 and constant birth and death rates per individual throughout its TABLE 5.2: Simulation Parameters for the Computation of Extinction Probability for a Population Simulation Death Rate Birth Rate Nature of Remarks No. H (years)'1 X1 (years)‘1 Boundary at a 1 0.6 0.6 absorbing For all simula— 2 0.6 0.6 reflecting tions 02/2 = 3 0.6 1.2 absorbing 0.025 inchz/year A 0.6 ZOsech(y-10) reflecting r = l inch/year ,5 0.6 SOsech(y-10) absorbing 5x = 0.2 inch 6 0.6 SOsech(y—8) absorbing 6T = 0.01 year 7 0.6 SOsech(y—6) absorbing a = 20 inches u____i, life-span, the changes in the total population will be the same as that for a linear birth and death process. The extinction probability for the linear birth and death process with one initial ancestor is equal to M - LL expf‘QfMt] i1 — n eXPE'Ql-nh] 113 when 11 # M and xlT/(l + xlT) when k1 = u; where x1, u and T denote the corresponding quantities defined in this section (Cox and Miller, 1965; p. 166). This result for a linear birth and death process checks with that of Simulation 2 for all initial states of the ancestor. Comparing the results of Simulations 1 and 2 it can be seen that when the initial location of the ancestor is sufficiently away from the absorbing boundary (note that this distance from the boundary will be dependent upon g3, r and n) the extinc- tion probability at any time is the same as that for a reflecting boundary. Thus, the effect of the boundary is apparent only near the boundary itself. Simulations 4 thru 7 were conducted to represent the situa- tion where each individual has a very high reproduction rate for a small portion of its life-span. In Simu1ations 4 and 5 the peak 1 reproduction was assumed to be at the midpoint of an individual's life-span. The results of both simulations show that althOugh the extinction probability immediately following the introduction of the ancestor in its peak reproductive state is quite low, it in— creases to a high value as time progresses. The results of Simula— tions 6 and 7 are more interesting. Simulation 6 shows that if the ancestor is "too young" or "too old” at the time of its introduction, the probability of extinction of the population is very high. In Simulation 7 a steady state was reached after a period of about ten years. Similar steady state solutions can be obtained for other simulations if the period of simulation is long enough. In general, it can be said from Simulations 4 thru 7 that for a population with constant diffusion and drift coefficients, a constant death rate 114 and the reproduction rate per individual having a functional form similar to that used in Simulations 4 thru 7, the extinction proba- bility at any given time is the lowest if the ancestor is in its peak reproductive state at the time of introduction. As mentioned in Section 4.5, an important application of the extinction probability is in the biological control of pest species. Simulations 4 thru 7 indicate that the maturity of the parasites or predators introduced for achieving the control may be a crucial factor in the survival of the controlling population and hence the chance of success of the biological control strategy. When the parameters in the diffusion equation are not constant the method described in the previous section cannot be used to find the ”critical length". In such cases computation of the extinction probability can be used to find the critical length by using the criterion that when 5 < acr’ the probability of ultimate extinction is one. Remark 5.3: Application of the results of this chapter to biological populations involve description of the growth of an individual by a diffusion process. The individuals in the population are often characterized by their maturity or a measure of their size such as a characteristic length, which are essentially non- decreasing quantities for each individual in the population. On the other hand, realizations of a diffusion process are not strictly monotone in nature. Thus, for example, when the growth of a fish in terms of its length is represented by a diffusion process, there will be a nonzero probability that the length of the fish as described by the model will decrease sometime during its lifetime. In this 115 sense the diffusion representation is essentially approximate in nature, and it is instructive to know the probability of a given decrease in the length, maturity, etc. of an individual if a dif- fusion model is used. Brockwell (1972b) has derived an expression for the probability of suCh deviations from monotonicity for a dif- fusion process on the real line (~m,m). In particular, let the growth of an individual in the maturity interval [0,5] be char- acterized by a diffusion process on the real line with constant dif- 2 fusion and drift coefficients 94- and r respectively, the death 2 parameter p = 0 and let z denote the location of the individual in [0,5] at time t. If Ta denotes the time taken by the individual to reach a for the first time (i.e., the first passage time) and M(t) = max z(T), then max(M(t)-z(t)) represents the 0STSt OStST . . . . a . . . . . max1mum decrease in maturity experienced by an indiv1dua1 during its lifetime. The distribution of this quantity is given by F (m) = P[ max a‘1(M-z> s w] = expE-c‘2-1)‘1>],<5.66> C OatsTa where c =-——Jg——. Thus, for the example of fish population char- (Zar)2 2 acterized by the parameters a, r, and 5L. in Table 5.2, (5.66) predicts that while the probability of a maximum decrease in length of greater than one inch during the lifetime of a fish is close to zero, the probability that a fish will decrease in length by more than one fourth inch during its lifetime is 0.036. This probability changes drastically with small changes in c. It is interesting to note that if the diffusion coefficient in the example just considered was 0.05 instead of 0.025, with the other parameters unchanged, so that the parameter c is 0.05 instead of 0.03536, the probabilities of a decrease in length greater than one inch and one fourth inch would be 8.24 x 10'7 and 0.93 respectively. 000.0 000.0 000.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 00.00 000.0 000.0 050.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 500.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 050.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.5 000.0 000.0 000.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 00.0 000.0 000.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 00.0 000.0 000.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 05.0 000.0 000.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 00.0 000.0 005.0 005.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 05.0 000.0 055.0 550.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 005.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 05.0 000.0 000.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 05.0 000.0 000.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 000.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 05.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 000.0 H0~.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 116 $6 85 8.0 8.0 $6 8.0 8.0 86 85 86 8.0 AI 8% suufim 8.8 8.3 8;: 8;: 8.2 8.2 86 8.0 8.0 80 8.0 “Lamb coco-occaooooaa-novoooao-o- AWQ£UCHV poummuc< mg“ Mo fiuwgwg HWHUHGH"?u.ouncanny-000.co-ouooooocv NEH—Ht .mmzo50 00 u m um pmflpumm wc00u0m0< .Haflpww00 0.0 n HA wumn Lupflm .aufiummkv 0.0 u 1 mum» Lumma .umm%\£ocH H u p mums susouw cmmz .umm%\0A£ocHV 000.0 M 0\00 ucwfioawmmoo coqmsmwflm .Mouwmoré mdozn woumuwcmo :oHumanaom w Mom kfiaflnmnoiam mofluocfiuxm I: H .02 :oHumHEfim mo mufismmm .00. MHLmH 000.0 500.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 500.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 500.0 500.0 500.0 000.0 000.0 000.0 000.0 00.00 000.0 500.0 500.0 000.0 000.0 500.0 500.0 500.0 500.0 500.0 500.0 00.00 000.0 000.0 500.0 000.0 500.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 117 0 0 c> 0° 00.0 00.0 00.0 00.0 00.0 00.0 00.0 00.0 00.0 00.0 00.0.A|.0MWMM 8.8 8;: 8.3 8;: 8.2 8.3 8.w 8.0 8.0 80 86 0308.0 noon.coo-unnoooooooo-ocoooo AmwfiwUCHV HOUWUUC< USU. MO ##UMCWJH HNWHHEHHW no.cocoons-ooooo-uooq-uooo- @EHH. Awwaswunoov .0.0 000mH 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.5 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 00.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 00.0 000.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 05.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 005.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 05.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 05.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 05.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 00.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 050.0 0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0 118 8.0 8.0 88 8.0 8.0 8.0 8.0 8.0 8.0 8.0 8.0+. 0”me 8.8 8.2 8.3 8;: 8.2 8.2 8.w 88 8.0 80 8.0 038.0 non.coo-onvv-uo-oc.c.ocuonn Amwsoaflv Houwwoc< GS“ we guwcwg HNHUHfiHHw ...-nvoouuoouc-uoooco-va-o- ”EH—H. .mw0050 00 u m um uw0uumm 0:0uom00mm -Auwmxv 0.0 n 00 womb LuMHm .0IAymw00 0.0 n 1 mumu :ummm .0mw0vSos0 0 u H wumu suBOMw :mwz .umm%\0ALUS0V 000.0 M 0\00 ucm0o0wwmoo COHmJMMHD .uOummoc< «so 00 Umumuocww coflum0saom m 000 0u0000mnoum EOHUoC0uxm u- 0 .oz CO0um05E0m mo wu0smwm .0.0 m0£mH 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0NJVN D 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0N.Nm 1 cmmd qmmé «N00 «N00 «No.0 qmmd «No.0 «N00 +No.0 «No.0 ¢N0.0 0N.0N 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0N.w0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 500.0 0N.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0N.q0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 0N.N0 00.0 00.0 00.0 00.0 00.0 00.0 00.0 00.0 00.0 00.0 00.04.... SMWMM 00.0w 00.00 00.00 0030 00.N0 00.00 00.0 00.0 00.8 00.N 00.0 5.3”»er 000.000.000.000...coco-o... Ammfiocflv Hoummoc\£uc0 0 u M wump Luzonm saw: .umo%\NALUC0V 0N0.0 u N\Nb ucm00000mou c00msmw0m .MOumwuc¢ mac %n vmumuwcwo CO0um0saom m 000 >u000nmnowm :00uuC0uxm In 0 .oz co0um03E0m mo mu03wwm .5.0 m0an 000.0 000.0 000.0 000.0 000.0 005.0 N05.0 N00.0 000.0 0w0.0 000.0 0N.qm 000.0 000.0 000.0 N00.0 qd0.0 005.0 055.0 000.0 500.0 000.0 «00.0 0N.NN 000.0 000.0 000.0 000.0 000.0 0N5.0 qQ5.0 050.0 000.0 000.0 000.0 0N.0N 000.0 000.0 000.0 000.0 0N0.0 000.0 005.0 N00.0 «00.0 050.0 N00.0 0N.w0 000.0 000.0 000.0 500.0 000.0 000.0 550.0 000.0 000.0 050.0 000.0 0N.00 000.0 000.0 000.0 qw0.0 000.0 000.0 «00.0 0Nw.0 0N0.0 050.0 000.0 0N.q0 000.0 000.0 500.0 000.0 500.0 0N0.0 000.0 005.0 000.0 «00.0 000.0 0N.N0 125 00.0 00.0. 0N.0 00.0 0N.M0 00.00 0N.M0 00.0 0N.0 00.0 00.0 Alllswwmm 00.0w 00.00 00.00 00.¢0 00.N0 00.00 00.0 00.0 00.0 00.N 00.0 Amkpmkv no......-...oo-nnc-o~.uo..o Ammgosflv HOUWWUC< MS“ MO £Uwcmq HWHUflCH H V nuu....‘.-..o.-no...-....- wEfipH 80:03:08 .5 .0 £08.. 000.0 000.0 000.0 0w0.0 NN0.0 000.0 050.0 N0m.0 0No.0 500.0 000.0 0N.00 000.0 000.0 000.0 000.0 500.0 000.0 m00.0 000.0 000.0 000.0 NM0.0 0N.00 000.0 000.0 000.0 500.0 000.0 000.0 000.0 00N.0 000.0 NNw.0 0N0.0 0N.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 05N.0 N00.0 N0w.0 0N0.0 0N.w 000.0 000.0 000.0 N00.0 000.0 000.0 5N0.0 w0N.0 000.0 005.0 000.0 0N.5 000.0 000.0 500.0 050.0 000.0 000.0 N00.0 wmw.0 000.0 005.0 000.0 0N.0 000.0 000.0 000.0 000.0 0Nw.0 000.0 000.0 00N.0 000.0 005.0 000.0 0N.0 000.0 000.0 000.0 mm0.0 005.0 000.0 500.0 NON.0 000.0 005.0 050.0 05.0 000.0 000.0 000.0 «00.0 005.0 000.0 050.0 500.0 000.0 0N5.0 000.0 0N.0 000.0 000.0 500.0 000.0 005.0 05N.0 000.0 050.0 N00.0 000.0 000.0 05.0 000.0 N00.0 500.0 w~w.0 500.0 0mm.0 000.0 N00.0 000.0 050.0 N0w.0 0N.m 000.0 050.0 oww.0 055.0 000.0 000.0 000.0 N00.0 050.0 N00.0 N55.0 05.N 000.0 000.0 005.0 005.0 000.0 000.0 000.0 000.0 000.0 N00.0 005.0 0N.N 000.0 N05.0 500.0 000.0 000.0 0N0.0 000.0 000.0 000.0 000.0 000.0 05.0 000.0 000.0 000.0 N00.0 000.0 000.0 5N0.0 N00.0 mwN.0 000.0 000.0 00.0 000.0 000.0 0m0.0 000.0 000.0 N00.0 0N0.0 050.0 00N.0 000.0 000.0 0N.0 000.0 050.0 000.0 0N0.0 00m.0 000.0 0N0.0 000.0 NON.0 000.0 000.0 00.0 000.0 500.0 000.0 000.0 00N.0 000.0 000.0 000.0 00N.0 000.0 000.0 05.0 w0N.0 00N.0 000.0 000.0 000.0 000.0 550.0 00N.0 w0N.0 00.0 000.0 000.0 000.0 000.0 500.0 000.0 000.0 000.0 000.0 000.0 000.0 0N.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 126 c c c P c \o N C) 00.0 00.0 00.0 00.0 00.0. 00.2 00.00 00.2 0.0.0 00.0 8.0 Iswwmm 8 00 00 2 8 3 8 S 00 S 8 S 8 0 8 0 8 .0 8.0 00.0 03093 ........................... AMQLUCMV HONWGUC< ®£U MO LUMCQQ HGHUHCH H V .........................-. WEMH .mwSUC0 0N u m um Hm0wpmm wC0LMOmn< . . -Aum000 Awumvnum000 n 0x mumu suMHm .0IApmev 0.0 n 1 mumu Lummm 0 .pmw%\LQC0 0 n u muMH £u30uw ammz .umm%\NA:u:0v 0N0.0 n N\Nb unm0u0mmwou c00m3000o .HOumquw $50 .00 wwumuwcwo :00um0saon0 m p00 03000300000000 :00uoc0uxm ..I 0 .oz SOHuM05EHm mo mu0smmm .w.0 @0900. 000.0 000.0 000.0 N00.0 000.0 «00.0 0NN.0 Nom.0 000.0 000.0 000.0 0N.0N 000.0 N00.0 000.0 000.0 NNN.0 000.0 000.0 000.0 000.0 0N.NN 000.0 000.0 000.0 000.0 000.0 000.0 m0N.0 000.0 N00.0 N00.0 000.0 0N.0N 000.0 000.0 000.0 000.0 500.0 000.0 00N.0 000.0 500.0 000.0 000.0 0N.w0 000.0 000.0 000.0 000.0 000.0 . 000.0 00N.0 000.0 000.0 000.0 000.0 0N.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 N00.0 000.0 000.0 0N.00 000.0 000.0 000.0 000.0 0N0.0 000.0 000.0 0mm.0 0N0.0 N0w.0 000.0 0N.N0 127 c c c 0 0 0 c ._I 00.0 00.0 00.0 0N.0 00.0 0N.00 00.00 0N.M0 00.0 0N.0 00.0.Allaswwmm 00.0w 00.00 00.00 00.00 00.N0 00.00 00.0 00.0 00.0 00.N 00.0 Awummmv ..o-o...oc-ootooo~...uuuo- AW$£UCHV Houmwocag mg“ “0 sumcw‘H HNWUHCH III W Icy....-.c-nuo-uouu03000un. wEfl'H. 00000000000 .0.0 00000 000.0 000.0 000.0 000.0 500.0 000.0 050.0 0N.00 500.0 000.0 000.0 000.0 500.0 000.0 050.0 0N.00 500.0 000.0 000.0 000.0 500.0 000.0 050.0 0N.0 000.0 000.0 000.0 000.0 500.0 000.0 050.0 0N.w 000.0 050.0 000.0 000.0 000.0 mom.0 N5o.0 0N.5 5N0.0 000.0 000.0 000.0 000.0 000.0 000.0 0N.0 000.0 000.0 N00.0 000.0 000.0 000.0 N00.0 0N.0 000.0 000.0 000.0 000.0 N00.0 000.0 000.0 05.0 050.0 0N0.0 000.0 000.0 000.0 000.0 000.0 0N.0 000.0 N00.0 000.0 000.0 000.0 000.0 N00.0 05.0 000.0 050.0 000.0 N00.0 000.0 000.0 0N0.0 0N.m 005.0 000.0 0N0.0 000.0 000.0 mmm.0 000.0 05.N 500.0 5w0.0 000.0 0N0.0 000.0 500.0 N00.0 0N.N 000.0 0N0.0 000.0 0N0.0 050.0 0wN.0 000.0 05.0 w00.0 000.0 N00.0 0N0.0 050.0 N5N.0 N00.0 00.0 500.0 000.0 000.0 NN0.0 000.0 50N.0 000.0 0N.0 5N0.0 000.0 050.0 0N0.0 N00.0 me.0 000.0 00.0 000.0 00N.0 500.0 000.0 500.0 000.0 mmm.0 05.0 00N.0 000.0 000.0 000.0 000.0 550.0 00N.0 00.0 000.0 500.0 000.0 000.0 000.0 000.0 000.0 0N.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.0 00.0 00.0 00.0 00.0 00.0 8.0 00.2 8.8 00.2 8.0 00.0410wwmm I n + 8.8 8 00 8 S 00.: 8:2 8.8 8.0 8.0 8.0 8.0 00.0 000003 ooc....-gc-oo-o‘n-oo...no-. Awwfiocflv POUWQUC< W3“ W0 fiumawg HNfiUflCHHwon...-0.00....uo0.0.00.0... NEW—H. .wwLUC0 0N u m um um0uumm wC0DMOmD< .0 Aumm>v 00100500m00 u 00 wumu nuu0m 0nAumw>V 0.0 H 1 mumu Lummm n.umm%\LQC0 0 n » wumu £03000 cmmz .MMm05NAzo:0V 0N0.0 u N\ b ucw0o0mmwoo :00m30000 N .pOuwmoc< mco >2 kuwprmU :00um0saom m pom 0u000£mnoum C00u050uxm an 5 .oz c00um0350m mo mu0smwm .0.0 m0QmH 129 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000 0 000 0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 000.0 00.00 00.0 00.0 00.0 00.0 00.0 00.0 00.00 00.00 00.00 00.0 00.0 .1: 0000 00000 0 00.00 00.00 00.00 00.00 00.00 00 00 00.0 00.0 00.0 00.0 00.0 0000000 ........................... Amwfiodwv Houwmodsw mSU m0 Luwfiwfl HMHuflSHHW........................... meH. 00000000000 .0.0 00000 CHAPTER V1 POPULATIONS WITH AN EXTERNAL INPUT In general, a pOpulation will also have an input other than reproduction by its members. In biological populations this external input usually takes place in the form of an input of individuals across the boundaries of a geographical domain due to locomotion in case of mobile species; or due to passive motion along with the carrying medium such as water for aquatic organisms or air in case of flying insects. The examples of external input in case of particulate processes in chemical engineering systems include an input of particles along with the carrying fluid for suspensions, nucleation of bubbles or crystals at the imperfections on the surfaces of a vessel and a homogeneous nucleation of crystals occurring at very high super— saturations of the magma. Formation of free radicals by the initiation reaction in polymerization processes can also be looked upon as an external input of molecules of zero chain length. Radcliffe (1972) considered the problem of external input (immigration) in the form of a nonhomogeneous Poisson process into a population where the individuals move within an abstract state space according to a Markov transition probability and multiply according to a branching process in which the offspring are produced at the end of life of the parent and are in the same state as the parent. Radcliffe also assumed that all the individuals entering the state 130 131 Space at any time due to migration are located at the same point in the state Space at the instant of their arrival. This description does not apply to many processes which are considered in the present work. In the processes described in the preceding paragraph the entities may be located at any point in the spatial domain at the instant of arrival and in the general process of reproduction des- cribed in Chapter IV an offSpring is in a state different from its parent. 6.1 Description of the Process For the analysis of populations with an external input, the individual state space E for the existing (or ”live") entities is taken as an n-dimensional Euclidean space. When the entities are removed from E by the process of "death" they are considered as * ~ N 7% being absorbed into a single point g é E. The union ‘E E E U g constitutes the individual state Space for the existing and dead entities together. The movement of entities is restricted to the spatial domain ‘3 E 5 U §*, where E is a closed nonempty set in E with a boundary F (cf. Section 4.2). The location of an entity in ‘3 is denoted by g. Only the entities in b may possibly move to another point in ‘3, or produce new entities by reproduction. In the context of a chemical process such as crystallization, E may consist of a four-dimensional Euclidean Space with three co- ordinates denoting the geographical location of a crystal and the fourth co-ordinate representing the size of the crystal in terms of its characteristic length. The domain 5 may refer to the pro- cessing system itself, Such as the interior of a crystallizer (with 132 the size co-ordinate possibly varying from 0 to 0). All the crystals that are removed from the system by flow or by mechanical means (i.e., by ”death”) are considered as being transferred to the absorbing state 5* E E. Similarly, in the case of a biological population E may consist of a four—dimensional Euclidean space of three geo- graphical co—ordinates and a maturity co-ordinate, 5 may be a given geographical domain of interest and "death” may be interpreted as the combined effect of natural death, emigration and harvesting. The process is started at time s 0 With ko entities in D. At subsequent times Sj’ j = l,2,...,r; kj entities enter 5 from a source external to the population. Each entity present in 5 at time 50 and the entities entering 5 due to external input at any time sj > s serve as ancestors of the entities generated 0 due to reproduction. The initial locations of these ancestors are denoted by the superscripted variable go. The movement of each entity in E ‘is characterized by a Markov transition probability x(A,t‘§,t0) denoting the probability that an entity at g at time to will be in the set A c:% at time t 2 t0. Moreover, the probabilities of producing any new entities due to reproduction at any time are also assumed to be Markovian in nature, i.e., the probability that an entity will produce any offspring at any time in future depends only on its present state and not the past history. It is also assumed that the input of entities to 5, their movement within. %, as well as their reproduction proceed independently of other entities in the population. It can be seen that the process described in this section includes the process discussed in Chapter IV as a special case. In the model considered in Chapter IV, E is 133 the nonnegative real line, the movement of an entity in 6 is char- acterized by a continuous Markov process satisfying (4.1) and the Markov probabilities of producing any offspring, or dying (i.e., of moving to §*) are defined in Section 4.1. For a complete description of the process the rate of input of entities and the points in 5 where they first appear must also be known. It is assumed that the arrivals of entities into the popula- tion due to external input occur according to k simultaneous independent nonhomogeneous Poisson processes, where k is finite. The Poisson inputs are characterized by the following axioms (see e.g., Parzen, 1962; Chapter 4): (i) The process has independent increments, i.e., the numbers of entities arriving during disjoint intervals of thne are independent of each other. (ii) For any time interval (no matter how small) there is a positive probability that some entities will arrive in 5, but it is not certain that an arrival will occur. (iii) The probability that one entity appears in 5 during the time interval (s,s+6s) from an external source through the i-th Poisson stream, i = l,...,k, is vi(s)6s + o(és), and (iv) The probability of more than one arrival in 5 during the time interval (s,s+és) from any of the input streams is o(és). It can be seen that in view of axiom (iii), the number of entities kj appearing in 5 at time sj must be equal to one for a finite number of Poisson inputs. When axiom (iv) is dropped so that the 134 probability of more than one arrival during (s,s+§s) is C(63), one obtains a generalized Poisson process. The material in this chapter can possibly be generalized to include an input in the form of a generalized Poisson process, but the analysis will be much more complicated. The location of an entity at the instant of its arrival in 5 at time s (which is obviously conditional on the fact that an arrival has taken place) from the i-th Poisson input is described by a probability density wi(§,s). The transition probability X(A,ti§0,s) can thus be looked upon as the probability that an entity is in the set A Cr? conditional on its first appearance in 5 at the point go at time 3. Remark 6.1: It can be seen that when the number of Poisson streams k is a constant, the net effect of the k Poisson pro- cesses with rates vi(s) and conditional probability densities wi(§,s), i = l,...,k, characterizing the location g where an entity entering 5 via the i—th stream appears first, is equivalent to a single Poisson process with rate k v(s) = v,(s) (6.1) i=1 1 and a conditional probability density k E Wi(§ )S)Vi(s) w(§,s) = i=—1————————— (6.2) k )3 vi(s) i=1 for the initial location of an entity at the time of its arrival. Similarly, when the process has a random number K of independent 135 Poisson inputs (fixed for each realization of the process) with all rates vj(s) and the probability densities wi(§,s) of the initial location of an entity being equal to V(s) and W(§,s) respectively, if K = k is given, the input streams can be replaced by a single stream with rate kV(s) and a probability density W(§,s) for the initial location of an entity at the time of its arrival. To summarize, the process described in this section is a Markov population process with a Poisson input. The main difference between the description here and the process considered by Radcliffe (1972) is that here an entity is assumed to appear at any point in 5 at the time of its arrival. Moreover, the locations of the entities produced by reproduction have been left unspecified in this description, and therefore can be arbitrary. As mentioned earlier in this chapter, Radcliffe assumed that all the entities are located at a single point in 5 at the time of their arrival, and that the entities reproduce according to a branching process. 6.2 The Probability Generating Functional of the Process For a quantitative description of the process, let 0 O O . O O . O 0 . s s ) N(A,ti§10;§20,--~,§k O: Ellyu'5gk11’ glr:"‘>gkrr’ 0’ ’ r o ’V‘; n . denote the number of entities in the set A C2D at time t, given . O that there was one ancestor arriving at each of the pOints gij ' ' ' t 0 denote at time Sj’ i = l,...,kj; J = 0,1,...,r. The p01n s gio the initial locations of the entities present in the population at 0 the initial location of the i-th entity arriving into D at time 0 . the start of the process at time s and gij, J = l,...,r, denotes S. from an external source. Using the conditional counting meaSure J 136 0 . k r’ r S ,sr), the PGF for the O 0 _ 0 N(A,t\§10,...,§k0,...;§1r,...,§ 0,.... 0 process can be defined as O 0 _ _ 0 O . _ C(9,ti§10,--',§k00,---:§1r,--~agkrr:SO,---:Sr) — 0 0 0 O 1 ... ;...; ,..., ; ,..., . E{exp L13 og e(§)[N(d§.ti§10. ,gk 0 510 gkrr s0 51):” 0 (6.3) In view of independence of the entities, 0 0 0 o _ _ N(A,tiglo,...,§k 0,...,glr,...,gk r, $0,...,sr) _ O r k. r J 0 2 >3 N(A,t\gi.; s.), (6-4) j=0 i=1 3 J from which it follows that O O _ . O 0 . = G(e,tlg10,...,gk 0,...,§1r,...,gk r, SO""’Sr) O r r k' J 0 n 1'1 G(9,ti§i.; s.), (6-5) j=0 i=1 J J where G(e,t\§2j,sj) is the PGF of a pOpulation starting with one ancestor at ggj at time 5j and allowing no external input. Equa- tions (4.7) and (4.6) are special cases of (6.® and (6.5) respectively when r = O and k0 = k. The co—ordinates g and g0 in (6.3) thru (6.5) correspond to z and y respectively in (4.5) thru (4.7). Note that when the input consists of a Poisson stream, by Axiom (iv) in the previous section k1 = 1 for i = l,...,r in (6.3) thru (6.5). For an input stream in the form of a generalized Poisson process, the numbers ki’ i = l,...,r can be greater than one. The conditional PGF (6.5) has to be summed over all possible initial locations g0 of the externally introduced entities at all instants 137 of arrival 3, 50 s s S t, and all input streams in order to obtain a complete description of the process. For the case discussed in Remark 6.1 the external input can be characterized by a single Poisson stream of entities. Thus, for the case where k is a constant, the PGF of the process can be written as 0 . 0 0 C(e)t‘§20,“°,§koo; 30) = E{exp Jfilog e(§)[N(dg:tig10)'"agkoo; 80]} o 0 , = Eiifi---ifiiexp [13100 e<§>£N]r which is the joint probability density of r independent random variables with densities v(si)/V(so,t), i = l,...,r. (Q.E.D.) In View of Remark 6.2, averaging (6.6) over sj conditional on r, the following equation is obtained for the case where the number of Poisson streams k is a constant: G(e,ti§20,...,§:0 0; so) k0 0v(S,)dS_ = H OG(e ”51030)“? H IEJ‘ OBI G(e, tigg ,Sj )W(§:.) ,Sj )dg? v(SO’t) i=1 ko _ X(e.:s ,t) - i10,I=1G(et|§io,so)0 V_(s_0.2)’_ (6.9) -¥__¥ 139 where t X(e,s0.t) = j j~cw<§°,s>vd§0ds. (6.10) 30 D Noting that r is a Poisson-distributed random variable with a X(e,so.t) parameter V(sO,t), averaging r ivzggjgy‘} over r yields E{[X(e.30,t)jr} 0 [W30’t)lrexi’[“’(so’t)][x(e’so’t)]r ——-——-—- = E v(s0,t) r00 r1 [v(so,t)jr = exp{X(e,SO,t) " v(sogt)}3 so that k 0 O G(e,ti§10,...,§:oo,so) = i21G(e,t\§iO,SO)exp{X(e,sO,t)-V(SO,t)}. (6.11) When w(§0,s) = 5(g0—g8), i.e., when all the entities arriving . 0 0 . . . in D first appear at g0 and when there are no entities in the population at time 50 (i.e., k0 = 0), the expression for the PGF reduces to C G(e,t|so) = exp{j v(s>[c(e,t\gg; s0) - 1]}, (6.11a) S 0 Which is essentially the result obtained by Bartlett (1966) for a population of individuals characterized by their age, as well as that of Radcliffe (1972). When the number of Poisson streams K is a random variable (fixed for each realization of the process) as specified in Remark 6.1, given K = k, the analysis is the same as that for a fixed k with v(s) and w(§,s) replaced by kv(s) and W(§,s) respectively. Thus, ‘140 k 0 0 O 0 X(93803t) r G(e,t\§10,---.§k 0;sO) = HG(e,ti§iO;so)E{[-\7—(—S—T] }, (6.12) 0 i=1 0’ where t 0 o 0 i(9,s0,t) = j f~G(e.t\§ ;s)0(§ ,s)kV(s)d§ ds 80 D t o o o = kf [~G(9,t\§ ;s)fi(§ ,s)§(s)d§ ds 3 D O E k§(9:30at) (6'13) and 0 t V(sO,t) = I kv(s)ds S 0 E k $280,t) . (6.14) Averaging over r yields 0 O C(e,ti§10.....§k 0; so) 0 k 0 ... ... = Ii C(e,t|gio;so)E{eXP[X(e,SO.t) - V(So»t)]} i=1 k0 = nC(e,t\gio;sO>E{exp[k(”§(e,sO.t> -’i7(so,t)>]} i=1 k 0 ,0 5% k 0 n g(e,t‘gio;sO)E{[exp(§(9,sO,t) - V(so,t))j }. (6.15) i=1 If K is characterized by a probability generating function Y(u) = E[(u)k], O s u s 1, then by averaging over k, (6.15) reduces to O 0 . C(e)t\§103'°'3§k 0330) o k 0 ~ m = n g(9,t\gio;sO)Y{expf§(e,so,t) - V(sO,t)]} . (6.16) i=1 141 6.3 Moment Distributions of the Population In order to obtain the factorial moment distributions of the pOpulation it is necessary to apply (2.13) to (6.11). Thus, when the number of Poisson inputs k is fixed, the expected number of F3 entities in the set A C1D is given by = a_- O O . M1(A, tlF10,...,§kOO;SO) 8g G(T‘+Qe\§10,'°‘:gkOOSSO)‘n=lag—.O k 0 i=1 t 0 O 0 + {f VI~M1(A,t\§ ,s)w(§ ,s)d§ ds}, (6.17) s D 0 when 9 = I(A\§). Thus, the expected number of entities in A<:<5 at time t is simply the sum of the expected number due to the entities present in the population at the initial time s and the ~ expected numbers arising from all the entities entering D during the interval (so,t) due to all the individual Poisson streams. The second factorial moment of the number of entities in the set N N A1 X A2,: D X D is 0 0 a__ 0 . (AIXAZat\§10’°°"gkOO’SO >= 2G(“+Qe‘§o ’°"’§k00’80)\n=1,C~0 M 2 () 662 k 0 _ o 0 o o _ iEIM(2)(A1XA2,t\§iO;sO )+ [O v(s)jD( M 2(A 1XA2,t\§ ;s)w(§ .s)d§ ds k + 20M (A t|§0 's f v(s)“ M (A t\§O-s)w(g0 s)d§0ds i=1 1 1, i0, O)S jg 1 2, 3 , k 0 O t +12 M 1.(A2,t:\tg 005302r v(s)j‘fi M 10(A1,t\§ ;s)w(§O ,s)d§0 ds =1 S0 142 kO kO 0 0 + 2 2 M1w(§ ,s)d§ ds i=1 k k 30 o o + z 2 M 1(A1.t\§ 00;s0)M1(A2,t\§.0;sO) i=1 j= —l J jti 144 + E n t {f v(s)I~Ml(Aj,t‘§0;3)W(§0,s)d§0ds} . (6.18a) J 1 s0 D The quantitative description of the population in terms of the PGF assumes a particularly simple form in the case where the entities in the population do not reproduce and when the number of entities in the population at the initial time 50 is zero. 6.4 External Input in Populations of Nonreproducing Entities As mentioned earlier, a large number of processes of chemical engineering interest involve populations of nonreproducing entities. It will be seen later that some problems related to biological populations can also be looked upon as the dynamics of nonreproduc- ing entities. . . . 0 ~ . Given a nonreproduc1ng entity at g E D at time 5, there will be only one entity in B at all times t 2 s in the absence of an external input, because if the entity is not in 5, it must be I * c a n u w in g , and thus the entity is always restricted to stay w1th1n D. 0 . . In view of this, the PGF G(e,t\§ ,s) for the populat1on start1ng . , 0 . . With one nonreproducing entity at g at time s and allow1ng no external input can be written as c = Ete<§),t\§°,s] L¥e(§)x(d€,t\§0,s), (6.20) D where 5 denotes the location of the entity (the "ancestor") at thne t (cf. Equation (2.9)). Substitution of (6.20) in (6.11) yields 145 0 C(e,t‘§10,...,§: O;sO) k0 O = E:(L§e(§)x(d§ tléo .8 O>>(€Xp{f V(S) f g X(dg, t\§% ,sO) 1: s 0 §€D § 6D w(§0,s)-1>d§°ds}> . (6.21) When there are no entities in the population initially, i.e., when k0 = O, and the number of Poisson inputs k is a constant, in view of Remark 6.1, C(e, tlso > = exp{jvw<§°,s>d§s = X(d§.t\§O.S)- D Case 2: When the initial location of each entity at the time of its first appearance in D due to the external input from the i-th Poisson stream is located on the boundary F of D with a proba- bility density wi(§,s), i = l,...,k; where E denotes a point on 148 P, the number of entities in any set A C5D at any time t 2 30 has a Poisson distribution with mean t . M1(A,c\so) = g: v(s);x(A,t\§,s)w(§,s)d§ds, (6.27) o where v(s) and w(s) are related to vi(s) and wi(s) by (6.1) and (6.2) reSpectively. In the case where the number of Poisson streams is a random variable with the process parameters V(s) and W(§,s) defined in Remark 6.1, substitution of (6.20) in (6.16) yields the expression t c = w(expif j j V(s)x(d§.t\§0.8)fi(§O.S)[e(§)-1]ds]} (6.28) 80 g6? 5.065 for the PGF of a population of nonreproducing entities with no entities in the population initially. When 9(5) is given by (6.22a) for the set A':(§, the PGF reduces to the probability gen— erating function ?(g,t|so) of the number of entities in the set A, given by t Y(Q,t\so) = ({exp(( j vx6<§°.s>d§°ds3 . (6.29) sogefi 6.5 Applications Most problems related to the population balance models of chemical engineering interest can be cast in the theoretical frame- work of the last section. For example, when a homogeneous nuclea-‘ tion of crystals or other suSpended particles (as in a flocculation process) occurs in a process vessel, a Poisson rate of production with a certain probability density for the location of the particle in the vessel at the time of formation may be a reasonable assumption. 149 When uniform conditions prevail throughout the contents of the vessel, a uniform probability distribution for the initial location of a particle at the instant of its formation can be assumed. A similar description can also be used to characterize atmospheric phenomena such as formation of smog. In many chemical processes such as crystallization, nucleate boiling of fluids, distillation and liquid- 1iquid extraction the particles (i.e., crystal nuclei, bubbles of a vapor or droplets of a liquid) are formed at the surfaces enclosing the fluid. In these cases a Poisson rate of formation with a certain probability density for the initial location of the particles may be a reasonable description. In case of processes carried out in perfectly mixed vessels, the size or weight of a particle is often the only quantity of interest -- in such cases the individual state Space for the particles in the population consists of the non-negative real line, and the formation of particles can be characterized as a Poisson input occurring at a single point in the state space, such as nuclei of zero size appearing in the vessel due to a homogeneous nucleation process. In many cases it may be necessary to consider the input of entities as consisting of a number of independent Poisson processes occurring in parallel. For example, in nucleate boiling of a fluid or in the case of a crystallizer where nucleation is taking place at the surfaces of the vessel, the local temperature at any point on the surface as well as the nature of the surface itself determines the rate of formation of the bubbles or crystal nuclei. Similarly, in processes such as distillation in a plate column the rate of bubble formation at the bubble caps or at the holes in a sieve plate will be dependent on the height of the liquid column 150 above the bubble cap or the respective set of holes in the sieve plate. Since the height of the liquid over a plate changes from the inlet to the outlet weir on the plate, the rates of bubble formation will be different at different points across the plate. As discussed earlier in this chapter, the stochastic descrip- tion of a population with an external input is quite involved in case of reproducing pOpulations and therefore application of the theory of Section 6.3 to a crystallization process will be quite difficult when homogeneous nucleation, nucelation at the surfaces of the vessel, as well as secondary nucleation are occurring simultaneously at significant rates. Most particulate processes of chemical engineering interest involve particles of nonreproducing entities, for which the simpler results of the previous section may apply. Thus, in general, it can be said that if (1) At the start of a particulate process there are no particles present in the process vessel. (ii) The input of the particles into the process vessel can be characterized by a fixed number k of Poisson streams, with a probability density assigned to the particles in each stream for their states (such as the location inside the vessel as well as other characteristics such as size, shape, weight, etc.) at the time of their first appearance in the vessel, (iii) The future states of a particle and also the probability of its removal from the population depend only on its present state and not the past history, 151 (iv) The changes in the state of a particle as well as its removal from the population proceed independently of other particles, then the number of particles in any given set of states such as a size class, a volume element inside the vessel, etc. at any time is a Poisson-distributed random number. In case of biological populations the external input commonly manifests itself in the form of a migration of individuals. As dis- cussed before, the analysis of the stochastic fluctuations in a population with reproductive as well as migratory inputs is quite complicated. However, the study of migration of a single generation of individuals is also of importance in many problems of practical interest and in these cases the simpler analysis of the previous section can be applied. The example of the cereal leaf beetle (CLB), Oulema melanopus (L.) may be cited in this context. The CLB is an economic pest of the small grain cr0ps in the North American con- tinent. It has only one generation every year. In early spring the hibernating adults emerge from their overwintering sites in wooded areas surrounding grain fields and move to the succulent grain crOps. They feed on the plant leaves for some time before starting their mating and egg-laying activity. Most of the damage to the plants is done by the larvae emerging from the eggs. Among other things, a detailed study of the movement patterns of the adults is necessary to arrive at efficient strategies for the control of the pest. Ruesink and Haynes (1972) studied the distribution patterns of the CLB adult pOpulations in grain fields under different population densities. Their observations indicate that the adults are scattered 152 in a field according to a Poisson distribution only at medium densities. The mean and variance of the number of adults caught in a number of samples in the same field are significantly different from each other at very high and very low densities indicating that the Spatial distribution is not Poisson at these densities. If the rate of emergence of the CLB adults is assumed to occur according to a Poisson process with the points of emergence (i.e., the overwintering sites) distributed along the boundary of the field with a certain probability density and if the movement of a beetle within a field as well as its death are assumed to occur independently of others according to a Markov transition probability, then the number of beetles in any section of the field at any time should have a Poisson distribution. This description may be a good starting point in the quantitative study of the migration of the CLB adults. It is likely that the assumption of a Poisson emergence rate is not justifiable at low densities (cf. Section 6.6) and there may be a significant interaction among the individuals at high densities. The analysis of external input to populations of nonreproduc— ing entities can also be applied to problems other than Spatial migration of biological populations. For example, in the study of the CLB pOpulation dynamics, the egg-laying activity of the over- wintering adults can be considered as an "external" input to the new generation of the insects. In this case the individual state Space E consists of the positive octant of a three dimensional Euclidean Space with co—ordinates 21,22 and 23; with 21 representing the maturity of an individual and 22,23 denoting its geographical loca- tion in the field. The egg-laying activity of an adult located at 153 the point (22,23) in the geographical domain constitutes an "external" input of individuals of new generation at the boundary point g = (0,22,23). Since the eggs and pupae are not translocated and Since the migration of larvae is small relative to the dimensions of the field, the (Markovian) transition probability X(A,t‘§,s) characterizing the movement of an individual in E can be looked upon as the transition probability of movement along the maturity co-ordinate alone. Thus, x(A,t|§,s) is independent of the geo- graphical co-ordinates of an individual. In View of this, the quantitative description of the population in a fixed geographical area can be written in terms of 21 alone. It should be noted that when the diffusion equation can be used to characterize the growth of an individual, the transition probability X(A,t\zl,s) is given by the solution of the backward diffusion equation BX(A9t\zl:S) _ as = a; x(A,t\z1,s) (6.30) 1 with the initial condition X(A,t\zl.t) = I amwz monEom mo .oz awn monEmm pmuwz LuHB mumm coHUHmoaH>o mHumom mme Hmmuoo mo mHthmc< HmoHumHumum .H.© pomH an .wam%Hmcm msu :4 pop . . omecoo mums am My CG>HM aficw w . CHMSU wwwm UHmH SUHLB monEmw mnu %Hco umSu uamoxm .H.© mHan pom pom: mums mfimm m£u mu . mHH anmu mwfiH # mm 3 Ho.¢ wq©.o N ¢m©.oH Hmm.w NN m H©.q mmH.m N qu.o on.q 0H m Hm.N mmm.m H Noo.H moo.N mH N Ho.q moq.N N woo.N Nwo.m NN o mu Hm.N wwm.MH H Hmo.m dom.N HN m 1.. H®.¢ HNq.m N mmo.n mo.o NN q H©.¢ mHo.m N qHN.m m.m NN m Ho.q Nqa.N N mHo.m NwH.m NN N Ho.q qwm.¢ N Hm.w m.q NN H Ho>mH «A . mocmonwcme NOH Awmusaeoov X Eowmmpm mU\NA ozv kmv\mwwm # AmoHan eoumV NX N mo mwouwmm mocmHum> cum: mmHmem mo .02 %mm moneom wchmHuwwm LuHB mama coHunoaH>o mHumom mmmH Hmonwo mo mHmszc< HmoHumHumum .N.o pouH can 160 Remark 6.5: It should be noted that the x2 test for good- ness of fit of a collection of random numbers to a given distribution is only an approximate test. Thus the results of the test for CLB oviposition data can only be interpreted to mean that there is no reason to reject the hypothesis that the number of eggs laid by a collection of ovipositing females each day does not differ significantly from the Poisson distribution. In the context of the particulate processes in chemical engineering systems, the number of "active sites" responsible for the generation of the particles (such as the number of surface imperfections on a heating surface in the case of nucleate boiling) is usually much larger than the number of "sources" (i.e., the CLB females) considered in the analysis above. Hence, a similar reasoning of a superposition of a large number of independent sequences of events converging to a Poisson process can be used to justify a Poisson input at least in an approximate sense in many particulate processes as well. CHAPTER VII AN EXPERIMENT BASED ON THE THEORY 7.1 Motivation It was shown in Chapter IV that a diffusion equation can be used to characterize the moment distributions of a population if the movement of each entity in its spatial domain 5 can be represented by a Markov transition probability satisfying (4.1), the birth and death processes are also Markovian -- defined as in Section 4.1 and the movement of an entity in D as well as its reproduction and death take place independently of others. It is therefore necessary to test the validity of these conditions before the use of the diffusion equa- tion can be justified in the modeling of a population. The inter- dependence between the members of the population may not necessarily be in the form of a physical interaction. For example, the insects feeding on the leaves of a grain plant may be located on different leaves of the plant with no physical interference between the feeding’ activity of each other. However, the feeding by each insect affects the condition of the plant in terms of the growth rate and succulence of the plant tissues, which in turn affects the new leaf surface area as well as the quality of food available to all the insects feeding on the same plant. The result obtained in Theorem 6.1 in the last chapter can be used to test whether the conditions of Markov transition probabilities and independence of entities in the pOpulation are 161 162 satisfied. To demonstrate this with the help of a simple problem of practical interest, a controlled experiment was performed on the cereal leaf beetle (CLB) populations. The life-cycle of the CLB can be divided into seven distinct stages -- eggs, four stages of development (commonly called "instars" by the entomologists)of larvae, pupae and adults. Nest of the damage to the crop is done by the larvae, mainly because they feed on the plant leaves in the early stages of development of the plant, thus severely affecting the plant growth at high larval densities. Thus, for the within-generation dynamics.of the CLB p0pulation and its impact on the plants the range of maturities of the CLB from the beginning of the egg stage to the end of the fourth larval instar is of interest. The survivals in the pupae and the adults are related to the supply of eggs for the succeeding season. When the egg-input occurs according to a Poisson process, the maturation and death of any individual are characterized by Markov transition probabilities and the change from one life-stage into another occurs at fixed levels of maturity, then by Theorem 6.1 the number of individuals in any life-Stage at any time must be a Poisson-distributed random number. Thus, if a Poisson input of eggs is generated and if the p0pulation of any particular life-stage at any time in a large number of replicate experiments do not have a Poisson distribution, then it can be said that the conditions mentioned above are not satis- fied for the population of that life-stage. 7.2 Methods and Materials The CLB eggs are about one millimeter in length. The larvae grow from approximately the size of an egg to a length of about six 163 millimeters. In an experiment aimed at testing the independence of individuals and Markov transition probabilities it is necessary to have a large number of replicate experiments and to count the popula- tions of different life-stages of the CLB at different points in time without disturbing the plants or the insects. Since the experi- ment is based on the measurement of the sample realizations of a stochastic process, the accuracy in counting the number of individuals is also very important. For the ease and accuracy of counting and to ensure a better control over the dynamic processes in the p0pula- tion, the experiment was conducted in a controlled atmosphere room in a laboratory. Two sets of experiments were performed: Experiment A and Experiment B. In Experiment A an abundant supply of food was provided to the insects so as to have little interaction between individuals. In Experiment B the food supply was severely restricted so that growth and death of an individual may be significantly in- fluenced by other individuals in the population. Thirty replicates of each experiment were run to have a sufficient number of sample realizations of the process. Although a large number of beetles laying eggs simultaneously lead to approximately a Poisson input of eggs (cf. Section 6.6), the experiment was performed with an artificially generated Poisson input by manually transferring the newly laid eggs onto the experimental plants to ensure a better control of the input. A sequence of 30 random numbers having a Poisson distribution with mean = 10 was cnmputed using a table of random numbers. Each number in the sequence represented the total number of eggs to be used in a replicate experi- ment. Adult beetles were allowed to lay eggs on plants enclosed in 164 a cage in the laboratory for a period of eight hours, thus generating approximately a cohort of eggs. The required number of eggs was then carefully transferred on the experimental plants. The eggs have a layer of an adhesive on the outside, which made the transfer of eggs relatively easy. The same sequence of Poisson-distributed random numbers was used to generate the egg input in Experiments A and B. Pots with barley plants having a uniform height of about 3.5 inches were chosen for the experiment. The number of stems per pot was chosen to be 40 for Experiment A and 25 for Experiment B. Plaster of Paris suspension was poured in each pot to cover the soil. The larvae have to go underground toikuuipupae. At the appropriate time a layer of vermiculite (a mineral normally used as a soil conditioner) was spread over the plaster of Paris surface. The plaster of Paris surface acted as a barrier restricting the insects to form pupae in the vermiculite layer. The required number of eggs were transferred to the plants and the pots were kept well separated from each other in large metal trays. Water was added to the trays every day to irrigate the plants through the holes at the bottoms of the pots. The atmospheric conditions in the room were set at 16 hours of daylight provided by fluorescent light and 8 hours of night, with the temperature being 80°F and 76oF during the day and night respectively. The relative hunidity in the room was kept at 40%. Under these conditions the eggs began to hatch on the seventh day. Several hours before the actual emergence of the larvae from the eggs, the eggs had turned dark brown and the heads of the larvae could be clearly seen through the shells. After the eggs were ready to hatch the number of individuals in each life-stage were counted every day. 165 Each larval instar of the CLB is characterized by its head capsule diameter, which changes only during transition from one instar to another, when the insect sheds its head capsule and a part of its skin, and forms a larger head capsule. The head capsule diameter changes from one instar to the next by only a few microns. It is therefore difficult to distinguish between an insect which would molt shortly and one which has just molted without a microsc0pe. More- over, the molting process lasts several hours. Since the counting was to be done without disturbing the insects or plants, the larvae had to be classified according to the instars only by a visual observa- tion. However, the presence of shedded head capsules and skins on the leaves often provided a clue regarding the number of individuals undergoing transition. The larvae can crawl from one plant to another. When the larvae grew large enough to crawl from one pot to another, the plants in each pot were enclosed in lantern globes covered with pieces of cloth. Plants in Experiment A were sprayed periodically with water by means of a spray bottle to simulate rain. However, this did not seem to affect the larval mortality. In spite of the limited food availability in EXperiment B, the larval mortalities were not noticeably different from that in Experiment A. In view of this the number of stems per pot in Experiment B were reduced progressively to four per pot. At this point in time however, all the available green foliage in almost all the replicate experiments was readily consumed by the larvae. After this stage only a limited amount of food (one stem per day or less) was supplied to the larvae in Experiment B. Figure 7.1 illustrates the difference in the condition of the plants in Experiments 166 m paw < muaoaHuooxm cmwzuwn aoHqucoo ucmHm oau «o newHummBoo .H.N ouswam 167 A and B on the 14th day for pots, each of Which had 14 eggs initially. The laboratory was also being used for other experiments, for which it was necessary to change the environmental conditions in the room. All the pots were therefore moved to growth chambers (i.e., small controlled—atmosphere cabinets) where the same temperature and day- light cycles were used. At the time of this transfer most of the insects had already formed pupae. All the pupae were left undisturbed until the adult emergence was complete, after which the pupal cases were removed from the vermiculite in the pots and Opened to make sure that any insects left inside the cocoons were dead. The important steps in the experimental procedure are Summarized in Table 7.1. The daily emergence of the adults did not Show the expected behavior —- in some replicates the emergence took place much earlier than the others, while in some replicates there was no emergence at all. The humidity in the growth chambers was not controlled. More- over, the presence of a stagnant mass of air enclosed by the lantern globe and the covering cloth over each pot created conditions Sig- nificantly different from those outside the lantern globes. The moisture evaporating from the soil was condensing on the inner sur- faces of many lantern globes, thus creating almost a complete reten- tion of moisture inside the globes. The presence of the condensate indicated that the air inside the globes was saturated with moisture. The temperature of vermiculite in the pots was also found to be sub- stantially different from that of the ambient air inside the growth chambers, and varied from pot to pot. The vermiculite temperature was found to be as high as 98°F in some pots. In general, the pots with drier soil tended to have higher temperatures. 168 Table 7.1. Important Steps in the Experiments with CLB P0pulations Day Procedural Step 1 Experiment started with 40 barley plants per pot 7 Eggs began to hatch -- daily counting of populations started 10 Put lantern globes on all pots '< 11 Sprayed water on the plants to simulate rain J.) G 12 Sprayed water on the plants to simulate rain 2 13 Sprayed water to simulate rain, covered lantern globes 'H with cloth H o 14 Sprayed water to simulate rain 2‘ 15 Added vermiculite to all pots for the insects to pupate “3 18 Moved all pots to growth chambers 33 Experiment terminated 1 Experiment started with 25 barley plants per pot 7 Eggs began to hatch -- daily counting of populations started 9 Put lantern globes on all pots 10 Reduced number of plants to 16 per pot a: ll Covered lantern globes with cloth U ‘12 Reduced number of plants to 9 per pot a 0 13 Added vermiculite to all pots for the insects to E pupate, reduced number of plants to 4 per pot -H H 14 Added one stem of barley as food 0 Q- 15 Added one stem of barley as food x a: 16 Moved all pots to growth chambers 31 Experiment terminated 169 The effect of moisture on growth of the CLB is not well understood. The vermiculite in many pots exhibited fungus growth. The fungus growth was more significant in pots with a higher moisture content. The fungus was probably responsible for a part of the pupal mortalities, and thus the nonuniform humidity in different replicate experiments may have resulted in nonuniform mortalities in the pupal populations in the different replicates. However, since Experiments A and B were conducted under almost identical environmental conditions, it is reasonable to assume that if the difference in the food supply to the larvae did not have any effect on growth and mortalities of the insects in the two experiments, the effect of nonuniformities in humidity in the different replicates of Experiments A and B should also be identical. The growth rate of an individual is directly related to the temperature (cf. Appendix A), which is believed to be a printipal reason for the deviations from the expected emergence behavior of the adults. However, the result of Theorem 6.1 is still applicable if the total number of adults that ever emerged from the pupae are chosen as the quantity of interest by considering the range of maturities corresponding to the entire pOpulation coupled with an additional absorbing state §** to characterize all the individuals that died during the adult stage (cf. Section 6.4) and by taking the time of observation to be such that no live individuals are present in any stage except the adults. In view of this only the total number of adults which emerged from each pot were analyzed instead of the daily number of emerging adults. 170 7.3 Results and Discussion The data on the daily population counts of each life-stage in Experiments A and B is listed in Tables 7.2 and 7.3 reSpectively. This data was analyzed using an IBM 1800 computer. The x? test was used to see whether the observed numbers in each life-stage at any given time had a Poisson distribution. The results of Experiments A and B are summarized in Tables 7.4 and 7.5 respectively. The com- puter program is listed in Appendix D. It can be seen from Tables 7.4 and 7.5 that the populations of all life—stages from the eggs to the pupae show a good fit to the Poisson distribution with few exceptions. The poor fit in these exceptional cases could be attributed to errors in distinguishing between successive instars. Better results were obtained in some instances when the dark-colored eggs which were ready to hatch were counted as larvae. Although the degree of starvation of the insects in quite a few replicates of Experiment B was extremely severe, the mortality as well as the X2 statistics of all stages up to the fourth instar did not exhibit any significant difference from that in Experiment A. All the fourth instar larvae which went underground were counted as pupae. Since the number of pupal cases recovered at the end of Experiment B was less than the number fourth instar larvae in many replicates, it is likely that some of the starved larvae could not form complete cocoons, or died in the early part of the pupal stage and were completely decomposed by microbial activity. The total pupal mortality in Experiment B was much higher than Experiment A. It can be seen from the results that while the final population of the emerging adults in Experiment A was not a good approximation to a Cereal Leaf Beetle P0pulation Data from Experiment A. Table 7.2. ----—--------—Populations in the Thirty Replicate Experiments--------- Life Stage ay D 7 12 9 8 5 4 7 8 5 6 10 10 14 8 9 l4 9 9 10 8 9 13 7 5 12 6 9 8 16 11 Eggs Unhatched 7 2 6 2 6 5 2 l 3 2 7 7 4 4 4 5 8 5 3 10 8 5 3 6 7 4 9 5 12 5 eggs Hatching 7 l 4 4 4 8 l 4 6 4 4 3 2 2 4 l 0 2 0 0 3 0 l 9 4 3 7 7 7 3 6 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 eggs Instar I l 7 8 Unhatched 3 2 1 0 2 2 1 l 4 0 l 0 0 2 0 l 2 0 0 2 0 1 3 eggs Hatching 8 171 1 0 0 l 0 0 l 0 7 9 6 3 9 5 4 10 8 l l 0 0 0 0 0 l l 0 1 0 0 0 0 O 0 6 9 8 7 3 3 7 5 3 4 6 7 12 4 8 ll 0 eggs 2 O O 2 3 5 7 13 11 Instar I 8 Unhatched 9 2 1 4 0 l 0 0 2 0 l 2 0 0 2 0 1 3 2 2 3 l 2 l 1 0 0 1 eggs Hatching 9 0 0 0 0 l 7 ll 4 8 10 8 7 9 6 3 9 5 3 8 O O l 0 0 O 0 0 0 l 2 0 0 0 0 0 l 0 0 eggs 2 0 0 l O l 3 5 7 14 ll Instar I 7 3 3 7 5 3 4 6 7 6 8 9 10 Unhatched 1 0 0 0 0 0 0 0 0 eggs 1 4 0 1 0 0 l 0 l 2 0 0 2 2 l l 2 2 0 l 0 Hatching 10 Od'c>d‘C>C>GIOJOJOJ OHQONNONT‘OON F‘CDQ'C>U)C>C>C>UWC>€) c>o~c>0\c>c>c>0\cwm v4<3l\C>d>C>¢)C>C>d'V)C>F~ 04C>SC>Q'W>C>S C>C>K>C>d>C>C>FiF~C>m> F*C>F)C>d'c>c>rid)w4e3 C>F‘O\C>53C>C>hl¢>0Jd> C>C>h~C>h-c>c>n4\o<3r\ nlC>€>C>€>C>C>d'°lC>¢> C>C>Q'C>Q'C>C>F'“3C><' C>C>F$C>FVC>C>OJF4U3C>C>P4‘3<3~¢ C>61K161U3C>C>°Q‘3!Hln C>FiolC>FSC>C>F4F46464 C>C>N3C>F3C>C>Pith>Ffi C>OJK)F4¢>C>C>~iUWcV~d C>C>¢>C>€>C>C>P4UWh~c>c3w4\o.4\o v4AIP~C>O\C>C>F)F-P4O\ ~a«)—JoJOJC>—4UN\D'H.A H4 H4 F4 Ov4\OCDP~C>C>d'd'C>h~ C>AIF)F*¢'C>C>OIFV¢3C>C>C>fll<3<fi eggs Instar I Instar II Instar I II Instar III Instar I Instar 111 Instar II Instar III Instar Instar II 10 10 ll 11 12 12 12 13 13 ll ‘-._ 172 I oomomOOOF-iOOChOOO‘OOChOON I I I OHMOGOMHOOQ'OOGOOQ'OQM I I : OOHOHOHOOOHOOHOOHOHH I : oomooooinmowxoooooooooowo I : COMOMOHNOOMOOMOOMOMN I : ooxooxoommooxooooooooor—I m 4; OMOOOOde'TOOOOOOOOOOOOOOOI‘ a) g ONdOOOOOOOOOOOOOOOOO 1.. g OOHONONOOMNOHQOONONO N [:1 OhflOOOOOOx‘I‘OOOOOOOOON H H H H H o . 4; OMdONONOONLfiOONOONONx‘T O :4 ONNHMOQOONNOHMOOQ‘OQM Q. g OSOHQHOOOOOOv-‘ONNOI-iooomxo f; OMQ'ONONOOHOOONOONONU‘) “ . E OMOHQOMOONMOHGOOUWOLnI—I a) ONNOQ‘OQ’OOHMOOdOOd'Ox'rm .C.‘ U ONOONONOOHHOONOONONI—i c'. "-1 OMHHmoq-OOI-Imooqoox‘roqq (D i g OQOHmHmOOd'NOv-‘mOOKOOO—i l ~.-I 4c; oHHHHHHoHov—IOONOONONI—I ...; a ONHOMOMOOHNOOMOOMOMM O ‘34 HQHNQNMHHHqHHqHHqu-N I I : ONd'r-imr-‘UWOHOMOOQOOOOOM I : OMQ'OOOKOOOHUWOOOOOOOGI—I I : OMNOmOd’r-JOOIAOOMOOUWOLOM I : Hx‘meOOONv—iOHBOOOOOOOOOOOO I I omMNONOOOthONOONONo I H H H H H I : OQ‘MONOQHOONOONOONONM I A : OMHOQ'OQ'OOHMOONTOOGOQ'Q' 'U I g : OHNomomHoomoomoomom—I C: ‘H H H H H H H H HH>H>H> I—I> I—I:> I—I.’> > g HHHHHHI—l HI—l I—II—I I—II—I H 0 HHHHt—IHH Hi4 Hi—I HM H U) V (D mmmmmmmmmmommommomou moo uuuuuuumuumuumuumumH ° ‘Hm mmmmmmmmmmmmmmmmamms “1 "49 cccccaasccgasgacacsp 5 H0) HHHHHHHQ-IHH HI—I HH Hn-Id: 9’ >3 H as dq-d'mmxoxoxonnnoooowo‘cmmoom % Q HHHHHHHHHHHHHHHHHNN ... i L 173 I S 2 ON or) H00 0 ooxoorxonOI‘IHoONO I I . . I In N NH N can 0 omoommoommoomo I I I I‘~ H HH 0 ON 0 ONOOHOOMOHHOOHO I v-| H H I—I I—II—I H I I 2 H I—I:> 0’ OJODbObOHGJODODbDI-IOJODDDbOHHHHHI—IHHHHHH *4 .EOOCIDD 4:00am 4:00:10!) ‘9 DID-HO)How-HGHUm-HOHHHHHHHHHHHH O a) u .2: mu 5 mu .t: mmmmmmmmmmmmm 000 com o um o um ,o uuuuuuuuuuuum .3. “:3 23s ‘6’ 2f “ 2i ‘6’ 2222‘28222228E? . m . [\I .40: but) :1: HD a: H: m I-IHHHI—II—II—II—IHI—II—II—Im 0’ >~ H m I—II\ I\ [\oo 00 com o oxoxOOHHNNmmx'fdx‘f g Q HHHHHHHHHHH [-4 H \l .L‘ I 0l\00\‘rMOI\0l\0I\0 I I I OMOOHNOMOMOMO I I I 0H000NmooI-I00r40 I H H H I I ONOOHMONONONO I I : c>ch>c>oIFIC>NIC>6IC>FIC> I : ONOOdMNmr-I0r-I0I-i m 4&1 Hnor—INOMLnI—IBOQO m E ONTOONNOdeTON'fO :4 0c: 0000000000000 é H0000HN0H0H0H r-4 H H H m g C)F-C>C)P~CDCVtnr4\0r4\0(3 o :3 C>WIC>C>P4<'C>U\C>U\C>V\C> D. ID 000005NO0NO0H00 In .H z» Hx'i'v-IOMC'WNN‘I'I-‘Lnoxoo H E; CJMIBICDQ'FIF4¢>C>F~C>F-C> m C>¢>C>C>cwaWConP~C>C>PI¢>C>F~C>F~C>P~C> I : c>EzC>c>mIrIq-¢>0Ia>04a>OI I : c:KIC>C>PI<'C>UIC>UIC>U)C> I : 0¢HOMN00000¢0 I : CIUIHIc>HIVIHIrIc>d-C>d'c> I I ONNoonmo. DI 9. >' g bib! hlhl PI PI PI PI 0 I-ILI LII-I I-I I-I I-I I-I \/ m m «30 m m m m m m w m m m won 44 U u-u u m u m u m u m u ° “IQ to m a.m m o.m n.m m.m m.m co 'HIJ G :15 G G 5 c D c 5 c 5 fi I; I403 hI+4nIhIhInIhIthInIhInqh4 m >» PI m InInIn\o\o\0I\I\uwH HH¢ qH mN.© MNN.H m mum.¢ mmo.N >H uwude qH mN.o mmm.w _m Ho¢.o MMN.m HHH naumcH ¢H mN.o mmm.m m 0mm.© 00H.o wMPMMH HH< mH mN.o nwm.N m nom.o ooN.m HHH uwumcH MH HN.N NoN.H H mm¢.o oo¢.o HH umumaH MH mN.o Nmm.m m mmm.o mmN.m mw>umH HH¢ NH mN.w me.N m ONm.¢ noo.d HHH umumcH NH Ho.q ooN.o N ONH.N mmm.H HH uwumcH NH mN.o owm.¢ m mwN.N mmq.o mm>umH HH< HH mN.o mwm.H m mwm.o Nom.m HH MmumaH HH mN.o mqm.¢ m 5mm.n Nom.c mm>umH HH< 0H Ho.¢ NmN.N N N¢N.H NOH.H H HmumcH +vmwwmu wcHnoumm 0H Ho.q oom.o N moN.H CON.H wwwm HH< OH mN.o Nmo.m m mam.o New.m HH HaywaH oH HN.N o¢N.H H Hom.o mmm.o wwwm wwsoumsaa OH m~.o qmw.m m qu.~ oom.o H “mumcfi + mmwm wafisoumm m HN.N wmq.N N mom.H mmq.H wwwm HH< m mN.o oo~.m m HHN.N Nom.o H umumaH m HN.N 500.0 H NmH.H OOH.H mwmw @msoumnab m m~.o . mmm.m m «qH.w NOH.N H “mumcfi + wwwm magnouwm m Ho.q HHN.o N Hmm.H mmo.H wwwm HH< w mN.o de.o m me.m 505.0 H umuwcH m :N 80.0 H owmd 03.0 $3 wfifiuwm m Ho.q NNN.H N HmH.H mmN.H .mwmm wosoumsa: w mN.o mom.mN m oom.o oo0.m H “mumafi + wwwm magnoumm N mN.o fifim.m m mqm.m oom.m mmwm HH< n mN.o qu.NH m omm.o noq.m wwwm magnoumm m mN.o mmN.m m “ma.“ mmo.m wwwm wmaoumna: N Hm>wH woconmHame NOH AcmusmEoov X Eovmmum wHummm msu AmmHan eoumv Nx Ho mmmuwmo mocmHum> cam: Ho mwmum mmHH mam < unmaflumaxm Scum muwm coHumHsaom mHummm mmmH Hamumu mo mHmkac< HmoHumHumum .¢.N mHan 176 H©.d MHm.HH N Nnm.¢ oow.N wuHDu< mm mN.m qnw.H m «do.m Noo.m mmasm 0N mN.© wa.H m md¢.o mm©.m mumsm mH 3.0 «$4 m End 35m $93 + 3.52 H2 3 mN.o moo.H m qo¢.o seq.m mmasm mH 2.0 gm: m .33. mfld 93% + 3:53 H2 2 3.0 H23 m H2. m oomé 92:5 S Ho.¢ «mm.HH N Hun.m mmq.H >H umumGH NH 3.0 «$4 m 30¢ m2.m mafia + $52 H2 3 HN.N wmm.H H BMH.H mmm.o mmasm 0H mN.@ Nwm.H m 000.0 mmm.q >H umumaH 0H mN.o qu.q m mom.o Nom.m >H umumuH mH HmpwH moconHchHw NOH Acmusaaoov x Eovwmum mHummm mSu N we mmmpwmo mocmHum> cmmz mo mmmum mmHH hmm AmmHame EoumV Nx AwmscHucoov .¢.n mHama 177 mN.o 0NH.N m Noo.o 000.0 mm>pwfi + wmasm mH HN.N NON.H H mm¢.o 000.0 amasm mH 3.0 mmNN 0 3.0.0 00H.0 >H .kumaH mH mN.o mNH.q m www.0 mmm.0 mw>hmH HH< «H mN.o me.m m MMH.w mmN.m >H umuwdH 0H HN.N mmN.o H H¢¢.N oow.o HHH umumcH dH mN.c oHn.N m OON.N mmo.© waPMmH HH< mH Ho.q mNN.¢ N mqo.N 00m.H >H kuwdH MH mN.0 NoN.N m 0Nm.n mmn.¢ HHH “wuwdH mH mN.0 NNH.m m mnw.0 ~0N.0 wm>pmH HH< NH mN.o omo.m m OHH.© 000.0 HHH uwumaH NH mN.o ooN.m m mNN.n 000.0 mm>mmH HH< HH 0N0 30.0 m 03.0 NON.N HHH 3005 HH mN.0 om¢.q m OON.N OON.¢ HH umuwcH. HH mN.0 Hum.H m on.N mmo.m wm>HmH HH< OH mn.m H¢o.H a 0N0.m Nom.0 HH umumaH oH wN.N NHH.H q mmo.w N©H.N wm>pwH HH< m Ho.¢ omq.NwH N www.mH ooo.m H HmumcH_+ wwww wdHaouwm m mN.o N¢O.NN m wmm.¢ OON.¢ HH umuwaH m H0.¢ onq.m0 N me.MH nmm.N H HmumcH m wN.N 0mm.m q 0HH.w mmq.n H HmumcH + wwwm waHsoumm w HN.N 0NN.¢ H 0H0.H 000.0 wwwo HH¢ w mN.0 mow.o m 0Nm.n HON.N H “wumcH m 3:5. 0:6 H Em; 02.0 00.. 3032.5 0 mN.0 HNN.N m ~00.0 N00.0 H HmumaH + mwmm wcHnoumm H wn.n ow¢.NH q 00N.N 000.0 wwww HH< N mN.0 nmn.n m Nm¢.0 N00.m H uwumaH N 3; 20.: 0 80.0 08.... .003 0.023% H 3.0 HmHH N 086 NEH .03 303200 H HmpmH wocmonHame NOH Eovmmum mHummm 0:» Amean Ecumv Nx Awmusgeoov Nx mo mmmuwmm mocmHuw> :00: Ho mmmum wHHH mmm m ucwEHuwaxm EOHH muwm coHumHzaom mHuwmm mmmH Hmwumo mo mmeHmc< HonumHumum .m.H mHan 178 El, Ho.q OHw.H¢ N mwo.m Noq.H muHsu¢ Hm mN.o oNo.N m oom.o Nom.w manna + mm>umH ON mN.o 00¢.m m OHH.0 000.0 mwmsm oN m.N.0 Oh0oN m 000.0 N0m.0 wanna + 2??ng mH mN.0 m00.¢ m mmm.m mmm.o mwasm 0H mN.o ¢NH.N m Noo.0 000.0 omada_+ wm>HMH mH mN.o o0N.N m mwN.m mmN.o mmasm wH mN.o qNH.N m N00.0 000.0 among + ww>umH NH mN.o Hmo.m m mmm.q ooN.m mmasm NH HN.N Nmm.o H MOd.H 000.0 >H HmumaH NH mN.0 ¢NH.N m Now.0 000.0 mwasa + mm>umH 0H mN.0 wwm.o m me.0 Nom.m ww>umH HH< 0H mN.0 on.0 m me.m mmo.m mmmnm 0H mN.0 0H0.o m Noo.0 mmm.m >H umumaH 0H m Hm>mH m mocon H0 H0 0 so mmu mHu mm msu .m. . NoH AvmusafiouV NX mo ammuwwm mommHum> cme we wwwum mHHH Nam AmmHan souuv Nx AvmscHucouV .m.N mHLwH 179 Poisson distribution (possibly due to the nonuniformity of temperature, humidity, etc. resulting in different death rates in the replicates), the adult population in Experiment B was even worse. Since the higher pupal mortality in Experiment B could be directly related to the amount of feeding done by the insects during the larval stages, i.e., the past feeding history, it can be said that the probability of death in the pupal stage was not Markovian, and therefore the Poisson char- acter of the adult population was further distorted in.Experiment B. Remark 7.1: In view of Remark 6.5, the goodness of fit of, the population data to the Poisson distribution can be inferred only in an approximate sense. The results of the x2 test can be taken to mean that the hypothesis regarding the independence of individuals and Markov nature of the dynamic processes in laboratory p0pulations of the CLB may be true for any life stage up to and including the fourth instar even under severely stressed conditions. If growth of an individual can be characterized by a continuous Markov process satisfying (4.1) and the probability of death can be defined as in Section 4.1, then the diffusion equation can be used to describe the changes in the population in the range of maturities from the beginning of the egg stage to the end of the fourth instar. This is the case when the egg input is known and only the impact of the beetle on the crOp for a given season has to be studied. Barr, Kharkar and Lee (1972) have used the diffusion equation to characterize the first moment density of field populations of the CLB. Remark 7.1 will hold for their model if the behavior of the insects in field populations is similar to that observed in the experiments. In practice, the CLB population densities in the grain 180 fields will never be allowed to reach a stage of a complete de- foliation of the plants. Chemical controls will normally be applied at a very much smaller level of damage to the plants than that observed in Experiment B. It is therefore possible that in field populations the conditions of independence of individuals and Markov transition probabilities will continue to hold even for the pupae and the adults at all the densities normally encountered, and the use of the diffusion equation may be justified for describing the entire life-cycle of the insect. However, it is also likely that under field conditions the solar radiation, wind and humidity will play a significant effect on the mortalities of larvae under stressed conditions. When there is a significant defoliation of the plants due to insect feeding, the larvae have to move larger distances on a plant or in a group of neighboring plants to find suitable food, and in this process they are more likely to be blown or shaken off the plant due to the action of wind. This increases the chances of a higher metabolic stress as well as mechanical damage to the tissues, re- sulting in a greater chance of dying at the time of molting due to inability of a larva to shed its head capsuleand skin. Moreover, the larvae on defoliated plants are directly exposed to the sun for longer periods of time and the effect of dehydration due to solar radiation is also likely to be more pronounced, particularly when the humidity of the ambient air is low. A careful investigation of the effect of these factors is needed before a definite conclusion about the justification of the diffusion equation to describe the CLB field populations can be drawn. CHAPTER VIII CONCLUSIONS 8.1 Summary of the Contributions of this Dissertation This work has cast the population balance models in a very general stochastic framework. Applications of the models considered in this dissertation include the residence time distribution analysis (RTD) of fluids in process vessels under steady flow conditions, dynamics of particulate and polymerization processes in chemical engineering systems as well as dynamics of populations in biological and social systems. A probabilistic treatment of the population balance models in chemical engineering has been restricted only to some simple cases thus far (e.g., Kane, 1971; Shinnar and Katz,r 1969; Schmalzer and Hoelscher, 1971). The "continuous stirred tank" and "dispersion" approaches in the RTD analysis have been shown to be applicable for a large class of general populations. Although analogous models had been prOposed and used for some biological populations, the generality of these approaches had not been fully appreciated so far. Since the methods of parameter estimation are well developed for the RTD analysis, the possibility of use of analogous models to other populations makes it possible to use analogous methods for the evaluation of parameters as well. Diffusion equations have been derived for the probability generating functional as well as the first two factorial moment 181 182 distributions and densities for a population of reproducing entities. These equations represent a more general situation than similar equa- tions previously reported in literature. The backward equations are used only very rarely in engineering applications -- the value of these equations in some cases of practical interest has been demon- strated. It is a common engineering practice to use Monte Carlo methods to characterize the stochastic fluctuations in the model of a process. These methods use models where the parameters are assumed to be random numbers with known distributions (the normal distribution is often used for this purpose). The values of the parameters are generated with the help of a random number generator and a number of solutions are obtained by using different values of the parameters thus generated. The mean and variance of the output of the model are then computed from the replicate solutions. COulman, Riece and Tummala (1971) have used this method in a model for a Species of freshwater shrimp to characterize the mean and variance of the population as functions of time. The equations for the moment dis- tributions derived in the present work represent an alternative to the Monte-Carlo techniques. The simple cases of practical interest for which the analytical and numerical solutions have been obtained in the present work represent new results and demonstrate the value of the stochastic approach to the populationbalance models. A very general analysis for populations of reproducing entities has been presented to include the effect of external input on the probability generating functional and the first two factorial moments of the p0pulation. In this analysis the external input is assumed to be in the form of a nonhomogeneous Poisson process, the 183 entities in the population are considered to be independent of one another and the movement of each entity in its state space is char- acterized by a Markov transition probability. An explicit partial differential equation has been obtained for the first moment density of a population of reproducing entities with an external input, with the movement of each entity in its one-dimensional state Space being characterized by a diffusion process. An important result (Theorem 6.1) has been obtained for the particular case of populae tions of nonreproducing entities, and its application to several cases of practical interest has been discussed. The models derived in this dissertation are based upon two key assumptions: mutual independence of entities in the population and Markov transition probabilities characterizing the death, re- production and movement of each entity in its state Space. Possi- bility of using the theory developed in this thesis in validating these assumptions for the particular case of nonreproducing popula- tions has been demonstrated by means of a simple experiment. 8.2 Areas for Future Research The analytical results in Chapter V were based upon the hypothesis that the dominant eigenvalue of the diffusion operator (with parameters defined by (5.1)) is unique and real. Although the existence of a real eigenvalue has been proved, it remains to be shown that this eigenvalue is the 23ly_dominant eigenvalue. The processes for which the present work does not apply include the size reduction processes in chemical engineering systems and degradation of polymer molecules, where every breakage of an 184 entity may result inoapoint process, i.e., a random distribution of sizes of the "pieces" may be generated at every breakage. Possi- bilities of extending the present analysis to describe these pro- cesses need to be explored. Another situation not covered in this dissertation is the case of a population with an external input in the form of a gen- eralized Poisson process. Incorporation of such an input in the analysis will greatly add to the generality of the present work. In industrial processes involving particulare matter the number of particles in the system is usually very large and con- sequently the stochastic fluctuations in the dynamic characteristics of the process system are almost always of a minor importance. However, the data needed for the design of these systems are usually collected in a small experimental setup where the number of particles is rather small, and the scatter in the data is often rather large. Use of stochastic models in the analysis of such systems instead of considering the scatter to be "noise" in the data may lead to a better understanding of these systems. Similarly, in systems such as fluidized beds often it is not possible to use a purely deter- ministic approach to characterize certain phenomena like the formation and breakup of large bubbles and the relationship of these phenomena to the characteristics of the fluidized particles. The stochastic framework provided in this thesis, or other formulations based upon the theory of stochastic population processes offer alternative ways to describe a broad range of such problems of chemical engineering interest. In order to fully appreciate the value of a stochastic approach such applications need to be explored in detail. 185 As mentioned before, the stochastic framework for the popula- tion balance models makes it possible to formulate the models based entirely upon phenomenological descriptions of the processes involved. In the analysis of crystallization processes the growth and nuclea- tion rates are commonly expressed by empirical kinetic expressions. Hypotheses treating these phenomena as stochastic processes would be a natural next step in gaining a better quantitative understanding of the crystallization processes. It would be beneficial to carry out similar investigations regarding the dynamics of other particulate processes also. The possibility of an application of the models developed in this thesis to social systems needs to be explored. For example, orthodox or enterprising nature of certain ethnic groups can possibly be quantified by using the concepts developed in Chapter III: Thus, if the post-world war II boom in the U.S. economy is considered as a known change (such as a step or a ramp change) in the available opportunities, a comparison of the economic (or educational) status of different ethnic groups can be made to quantify how well each group has used the opportunities. A group significantly surpassing (or lagging behind) the rest of the population in reaping the benefits from the opportunities can be seen to make the system a "nonideal population-flow system" (see Section 3.4), since all the individuals in the population do not progress at the same rate. The Opportunistic groups in the population can be characterized by the "bypassing" phenomenon whereas the groups lagging behind the rest of the popula- tion can be described by the phenomenon of ”dead space" in the population. The quantitative models for the social systems based 186 upon such considerations can then be used to derive alternatives to the current means of social control aimed at achieving a more equit- able progress in the society. 0f the possible boundary conditions for diffusion equations, the absorbing and reflecting boundaries have been used extensively in various applications. The residence time distribution analysis of flow in short packed beds offers a unique possibility of applica- tion of a third type of boundary condition describing the possibility that a fluid element at a boundary can jump to an interior point, as pointed out in Section 4.2. The application of such a boundary condition needs to be studied in greater detail. The formulations presented in this work are based upon the assumption that the entities in the p0pulation are independent of one another. However, many phenomena in the dynamics of biological populations result from interactions among individuals in the form of competition for a certain vital resource. Similarly, chemical engineering systems such as liquid-liquid extraction and gas absorption involve liquid droplets and gas bubbles which coalesce and break-up again, thus violating the assumption of independence of the entities. Application of the theory of stochastic population processes to pOpulations of interacting entities represents an un- charted area in the Study of stochastic processes. Development of even approximate techniques for the calculation of the first two moment distributions will represent an important milestone in the mathematical modeling of populations. REFERENCES REFERENCES Abkin, M.H. (1972). Policy making for economic development: a system Simulation model of the agricultural economy of southern Nigeria. Ph.D. Thesis, Michigan State University, East Lansing, Michigan. Adke, S.R. (1964a). A stochastic population diffusing on a finite interval. J. Indian Stat. Assoc. gz32-40. Adke, S.R. (1964b). The generalized birth, death and diffusion pro- cess. J. Math. Analysis and Applications 2:336-340. Adke, S.R. and Moyal, J.E. (1963). A birth, death and diffusion pro- cess. J. Math. Analysis and Applications ;:209-224. Barr, R.O., Kharkar, A.N. and Lee, K;Y. (1972). POpulation balance models in ecological systems. Presented at the 65th Annual Meeting of the American Institute of Chemical Engineers, New York City. Paper 82c. Bartlett, M.S. (1966). "An Introduction to Stochastic Processes". (Cambridge University Press, London). Bartlett, M.S. (1969). Distributions associated with cell populations. Biometrica 22:391-400. Beckman, W.C. (1940). Increased growth rate of rock bass, Amblqplites rupestris (R.), following reduction in the density of the popula- tion. Trans. American Fisheries Soc. 19:143-148. Billups, C., Wilkins, B. and Pike, R.W. (1971). Application of population balance models in the prediction of Shrimp size dis- tributions. Coastal Studies Bulletin No. 6, Special Sea Grant Issue, Center for Wetland Resources, Louisiana State University, Baton Rouge, Louisiana. Bischoff, K.B. (1966). Mixing and contacting in chemical reactors. Industrial and Engineering Chemistry §§: No. 11, pp. 18-32. Brockwell, P.J. (1966). Stochastic problems in tranSport theory. Report ANL-7l31, Argonne National Laboratory, Argonne, Illinois. Brockwell, P.J. (1972a). Department of Statistics and Probability, Michigan State University, East Lansing, Michigan. Personal communication. 187 188 Brockwell, P.J. (1972b). Deviations from monotonicity of a Wiener process with drift. Submitted to J. Applied Probability. Buffham, B.AH, Gibilaro, L.G. and Rathor, MQN. (1970). A proba- bilistic time delay description of flow in packed beds. AIChE Journal 16:218-223. Carslaw, H.S. and Jaeger, J.C. (1948). "Operational Methods in Applied Mathematics". (Oxford University Press, London). Connin, R.V. (1971). Department of Entomology, Michigan State University, East Lansing. Personal communication. Cooper, G.P. and Latta, W.C. (1954). Further studies on the fish p0pulation and exploitation by angling in Sugarloaf lake, Washtenaw county, Michigan. Papers of the Michigan Academy of Science, Arts and Letters 22:209-223. ' COOper, G.P., Latta, W.C. and Schafer, R.N. (1956). P0pulations of game fish and their exploitation by angling in several Michigan lakes. Presented at the Annual Meeting of American Fisheries Society, Las Vegas, Nevada. Coulman, G.A., Riece, S.R. and Tummala, R.L. (1971). Population modeling: a systems approach. Science 172:518-521. Cox, D.R. and Miller, H.D. (1965). "The Theory of Stochastic Pro- cesses". (Wiley, New York). Curl, R.L. (1967). Statistical treatment of particle processes. Industrial and Engineering Chemistry 22: No. 11, pp. 53-54. Daley, D.J. and Vere-Jones, D. (1972). A summary of the theory of point processes. In "Stochastic Point Processes" by Lewis, P.A.W. (ed.), (Wiley, New York), p. 299. Davis, A.W. (1965). On the theory of birth, death and diffusion processes. J. Applied Probabiligygz293—322. Davis, A.W. (1967a). Branching diffusion processes with no absorbing boundaries - I. J. Math. Analysis and Applications ;§:276-296. Davis, A.W. (1967b). Branching diffusion processes with no absorbing boundaries - II. J. Math. Analysis and Applications Lg:1-25. Feller, W.-(1954). Diffusion processes in one dimension. Trans. American Math. Soc. ;;:l-31. Feller, W. (1966). "An Introduction to Probability Theory and Its Applications -- Vol. II". (Wiley, New York). Funderburk, J.0. (1969). Design models for continuous emulsion poly- merization. Ph.D. Thesis, Iowa State University, Ames, Iowa. 189 Goodman, L.A. (1967). The probabilities of extinction for birth and death processes that are age-dependent or phase-dependent. Biometrika 24:579-596. Hall, D.J., Cooper, W.E. and Werner, E.E. (1970). An experimental approach to the production dynamics and structure of freshwater animal communities. Limnology and Oceanography g:839-928. Helgesen, R.G. (1969). The within-generation population dynamics of the cereal leaf beetle, Oulema melanopus (L.). Ph.D. Thesis, Michigan State University, East Lansing, Michigan. Hille, E. and Phillips, R.S. (1957). "Functional Analysis and Semi- Groups". American Math. Soc. Colloquium.PubliCations, Vol. 31. Himmelblau, D.M4 (1970). "Process Analysis by Statistical Methods". (Wiley, New York). Himmelblau, D.MJ and Bischoff, K.B. Simulation". (1968). (Wiley, New York). "Process Analysis and Hulburt, H.M. and Katz, s. (1964). same problems in particle technology -- a statistical mechanical formulation. Ebesieel Engineering Science L2:555-574. Indritz, J. (1963). "Methods in Analysis". (Macmillan, New York). Jagers, P. (1971). Notes on random measures and point processes. Unpublished notes, Department of Statistics, Stanford University, Stanford, California. Kane, 8.0. (1971). crystallizer. Secondary nucleation of ice in a stirred batch Sc.D. Thesis, Massachusetts Institute of Technology, Cambridge, Massachusetts. ' Katz, S. and Shinnar R. (1969). Particulate methods in probability systems. Industrial and Engg. Chemistry_§L: No. 4, pp. 60-73. Kendall, D.G. (1948). On the role of variable generation time in the develOpment of a stochastic birth process. Biometrika ggz3l6-330. Kendall, D.G. (1949). Stochastic processes and population growth. J. Royal Stat. Soc. Series B. ;;:230-264. Kharkar, AWN. (1971). Theoretical development of a mechanistic model for the dynamics of populations. 'MSU Systems Science/ Ecology Report, Michigan State University, East Lansing, Michigan. Kimura,‘M4 (1964). Diffusion models in population genetics. J, Agplied Probability ;:177-232. Leftkovich., LP. and Currie, J.E. (1963). The effect of food shortage upon larvae of Lasioderma serricone (F.). Bu11.Entomological Research 24:535-546. 190 Levenspiel, O. and Bischoff, K.B. (1963). Patterns of flow in chemical process vessles. Advances in Chemical Engineering gz95-198. Levenspiel, O. and Smith, W.K. (1957). Notes on the diffusion-type model for the longitudenal mixing of fluids in flow. Chemical Epgineeripg_Science 2:227-238. Moyal, J.E. (1962). The general theory of stochastic pOpulation pro- cesses. Acta Mathematica Lgézl-Bl. Moyal, J.E. (1964). Multiplicative pOpulation processes. J. Applied Probability ;:267-283. Parzen, E. (1962). "Stochastic Processes". (Holden-Day, San Francisco). Pawula, R.F. (1967). Generalizations and extensions of the Fokker- Planck-Kolmogorov equations. IEEE Trans. Information Theory Radcliffe, J. (1972). An immigration super-critical branching dif- fusion process. J. Applied Probability 2:13-23. Ramakrishna, D. (1972). Population balance modeling of mass transfer in lean liquid-liquid dispersion. 7lst National Meeting of the American Institute of Chemical Engineers, Dallas, Texas. Pre- print lle. Ramakrishnan, A. (1950). Stochastic processes relating to particles distributed in a continuous infinity of states. Proc. Cambridge Philosophical Society 22:595-602. Ramakrishnan, A. and Srinivasan, S.R. (1958). On age distribution in population growth. Bull. Math. Biophysics ggz289-303. Randolph, A.D. and Durando, A.K. (1971). Population balance models for predicting grinding mill performance. American Institute of Mechanical Engineers' Centennial Annual Meeting, New York City. Preprint 71-B-78. Randolph, A.D. and Larson, M.A. (1971). ”Theory of Particulate Pro- cesses -- Analysis and Techniques of Continuous Crystallization". (Academic Press, New York). Read, KJLHQ. and Ashford, J.R. (1968). A system of models for the life cycle of a biological organism. Biometrika §§:211-221. Ruesink, W.C. and Haynes, D.L. (1972). Sweepnet sampling for the cereal leaf beetle, Oulema melanOpus (L.). To appear. Schmalzer, D.K. and Hoe1SCher, H.E. (1971). A stochastic model for packed bed mixing and mass transfer. AIChE Journal 17:104-110. WIIIIIIIIIII 191 Sevast'yanov, B.A. (1958). Branching stochastic processes for particles diffusing in a bounded domain with absorbing boundaries. Theory of Probability and Applications 3:111—126. Sevast'yanov, B.A. (1961). Extinction conditions for branching stochastic processes with diffusion. Theory of Probability and Applications 2:253-263. Sinclair, C.G. and MCNaughton, K.J. (1965). Residence time prdbability density for complex flow systems. Chemical Engineering Science Sinco, J.W. and Streifer, W. (1969). Applying models incorporating age-size structure of a pOpulation to daphnia. Ecology 2g:608-615. Sinco, J.W. and Streifer, W. (1971). A model for populations re- producing by fission. Ecology 2;:330-335. Srinivasan, S.K. and Mehata, KHM. (1972). A Markov chain model of a packed bed. AIChE Journal ;§:650-652. Stuart, R.N. and Merkle, T.C. (1965). Calculation of treatment schedules for cancer chromatography, part II. E.0. Lawrence Radiation Laboratory, Report UCRL-14505. University of California, Berkeley, California. Takahashi, M. (1968). Theoretical basis for cell cycle analysis -- II. J. Theoretical Biology_l§:195-209. Tsuchiya, H.M., Frederickson, A.G. and Aris, R. (1966). Dynamics of microbial cell populations. Advances in Chemical Engineering 9:125-206. Wang, J.Y. (1960). A critique of the heat unit approach to plant response Studies. Ecology 2;:785-790. Weiss, G.H. (1968). Equations for the age structure of growing populations. Bull. Math. BiOphysics £92427-435. Wellso, 8.0. (1972). Department of Entomology, Michigan State University, East Lansing. Unpublished data. White, E.T. and Wright, P.G. (1971). 'Magnitude of size dispersion effects in crystallization. Chemical Engineering Prggress Symposium Series 2;: No. 110, pp. 81-87. Wright, P.G. and White, E.T. (1969). Size distribution studies in sugar crystallization. Proc. Queensland Society of Sugarcane Technologists, 36th Conference, Maryborough, Australia. pp. 299-310. Yun, Y.M. (1967). Effects of some physiological and biological factors on the reproduction, development, survival and behavior of the cereal leaf beetle, Oulema melanopus (L.), under laboratory con- ditions. Ph.D. TheSis, Michigan State University, East Lansing, Michigan. APPENDICES APPENDIX A SIMPLIFICATION OF THE DIFFUSION MODEL FOR SOME BIOLOGICAL POPULATIONS The physiological processes such as growth and reproduction in poikilothermic ("cold-blooded") organisms exhibit a common char- acteristic dependence on temperature. The rate of these processes is negligible below a certain threshold temperature TL, Which de- pends upon the species and the process, and increases more or less linearly with temperature up to an upper threShold Tu. The rate decreases with temperature above the upper threshold. However, for a Species adopted to a particular climate, the threshold temperatures are often such that the total length of time spent by the organism at temperatures above the upper threshold during a year is rather small. For example, the upper and lower threshold temperatures for the rates of growth and egg-laying for the cereal leaf beetle, a major pest of the small grain crops in midwestern United States, are about 48°F and 90°F reSpectively (Yun, 1967). The number of days in a year with temperatures in the nineties is indeed small in this part of the United States. The commonality in the temperature dependence of poikilo- thermic organisms is used to define a physiological time in terms of ”degree days”. The use of degree days often leads to clearer bio- logical insights into a system (see e.g., Wang, 1960), and is widely used by biologists to study population data. The profound effect of I 192 193 temperature on the physiological processes in poikilothermic organisms implies that except for constant temperature studies, the para- meters in the mathematical models of such populations must be functions of temperature. Barr, Kharkar and Lee (1972) have shown that the temperature dependence of the model parameters can often be eliminated by using degree days instead of chronological time in the quantitative descriptions. They have discussed this Simplifica- tion for the forward diffusion equation for the first moment density of a population where each individual is characterized by only one descriptor, namely, maturity. The same argument can be used to show that the diffusion equation for the probability generating functional (PGF) can also be greatly simplified by the use of degree days, and hence the simplification applies to all the moment distributions of the population. Consider the backward diffusion equation for the PGF with y denoting the initial location of an ancestor at time s: accessing) ._ azuusi 32G(e,tly,§l + r( AG(e,gly.S) .. . — ' Y’s) as 2 2 BY BY " M(Y,S)G(9>C\Y,S) + ”(Yas) ' z Ai(y,S)G(est\YaS) i=1 + z xi(y,s>c(e,t|o,s>ic (4.12) i=1 with the initial condition G(e.t|y.t) = 9(y) (4.16) and boundary conditions AGCSLtIYfiQ 0, (4,14) ay ‘y=0 = 194 and aG(e.tLy¢0‘ _ 3y y=5 for a reflecting barrier at y C(e,t\§,s) = l (4.13) for an absorbing barrier at y 2 o I E-, r, M and xi, i = 1,2,... be piecewise linear functions of 5. Moreover, let the parameters temperature with the same lower threshold temperature, i.e., for temperatures T(S) below the upper threshold let 2 2 113$)- = f[T(s>]§-—. - (A.l) r(y.S) = f[T(s)]f(y). (A-Z) u(y.S) = f[T(S)]E(y), (A-3) and xi(y,s) = f[T(s)]ii(y), i = 1,2,. . (A.4) where T(S) - TA for T(S) 2 TL f[T(s)] = (A.5) 0 for T(S) 3 TL and define the physiological time variable ? as s E = l f[T(w)]dw , (A’6) for a fixed t. Since d? = f[T(s)]ds, it follows that (4.12) thru (4.16) can be written as 195 ~ ~2 2 ~ 0 _ AG(Q:_t LXL’T) = LL11 LG ( eitl‘lz'r) + "f (y) ACKGit ly LT) aFF 2 BY2 ay - a + a - z ii(y>c(e,t\y,$> i=1 + 2 iic(e,t\o,$)ic = e, (A.8> aGiefifla,f) = A 9 by \y=0 0’ ( ° ) aG£9,tLa.?)‘ _N = o (A.10) ay y—a and G(9,t\§,¥) = 1 (A.11) respectively. Note that the temperature dependence of the parameters has been eliminated in (A.7) thru (a.ll) because f[T(s)] has been cancelled from each term in the equations. In view of the fact that xi(y,s), i = 1,2,... essentially describe a physiological process, (A.4) is a realistic approxima- tion for reproduction in poikilothermic organisms. The actual temperature dependence of growth or maturation would very likely be such that the random increment by in the maturity of an arbitrary individual during a time interval 63 ,is a piecewise linear func- tion of temperature. The incremental growth 6y is related to the diffusion and drift coefficients by (4.1), from which it follows that r would be characterized by a piecewise linear function as 2 in (A.2) and %_ would be a quadratic function of temperature. Analysis of the data of Yun (1967) and Helgesen (1969) for maturation of the cereal leaf beetle lends some Support to this conjecture. 196 In spite of the quadratic dependence of the diffusion coefficient on temperature, it may be possible to obtain a fair piecewise linear 2 approximation for g—' as shown in (A.1). The term u(y,s)ds is the probability of an individual of maturity y being removed from the population by death or emigration during the time interval (s,s+ds). In fact, it can be viewed as a sum of three different terms representing loss of individuals to the population from (i) emigration, (ii) death due to old age, disease, etc. and natural enemies such as parasites and predators, and (iii) artificial removal by harvesting, pesticide application, etc. The following discussion considers (i), (ii) and (iii) in sequence as if each occurred alone: When emigration of individuals occurs solely due to locomo- tion, (A.3) can be justified in view of the fact that locomotive activity of an individual is essentially a physiological process and hence dependent upon temperature. For natural death and death due to natural enemies, (3.4) does not at first seem to be an accurate re- presentation. However, in some species these mortalities are often considered to occur at those particular discrete maturities which define the transitions between the various life stages of an organism. For such a situation (4.12) may be used to describe the behavior between each of these transitions with u(y,s) being 0. Thus (A.3) is trivially Satisfied with fi(y) = 0. When this form of analysis is employed, these mortalities are used to determine boundary condi- tions at each of the transitions. The assumption of mortalities 197 occurring only at maturity transitions seems to be particularly good for natural death and parasitism, as compared to predators. Never- theless, since a prey may be less mobile at transitions between life Stages, it could be more vulnerable to predators. Artificial re- moval by harvesting and pesticides is essentially independent of temperature and thus for this case (A.3) is invalid, and the func- tion u(y,s) in (A.7) must be replaced by fLEQTLE-g—J— . (A.12) If n(y,s) = 0 whenever f[T(s)] = 0, then the analytical form of (A.3) and the resulting (A.7) are well posed. This supposition is true when harvesting, pesticide application, etc. do not occur at temperatures below the lower threshold TL. When f[T(s)] > 0, (A.5) describes a one-to-one relationship between s and E, and thus in (A.12) 3 may be replaced by ?. The resulting form of (A.7) is simpler than (4.12) in the sense that only one parameter is dependent upon temperature. When (i), (ii) and (iii) occur simul- taneously, all the remarks above must hold for the simplification to be possible. When (A.1) thru (A.4) are valid even after a suitable choice of a maturity variable y to remove the dependence of parameters on y, (A.7) reduces to an equation with constant coefficients, and the analysis of Sections 5.2 and 5.3 can be used to solve the moment equations derived from the diffusion equation for the PGF.’ In practice, 2 the parameters gr, r, n, and *1 have to be taken as being in- dependent of maturity over the respective maturity intervals. For example, the data of Yun (1967) and Helgesen (1969) on the growth 198 characteristics and egg laying behavior were obtained in the form of the mean and standard deviation of the develOpmental times for each of the life stages, and mean number of eggs per adult per day respectively. In such cases elimination of temperature dependence of the parameters again leads to a diffusion equation with constant or piecewise constant parameters. APPENDIX B NUMERICAL COMPUTATION OF THE "CRITICAL LENGTH" As discussed in Section 5.4, the critical location a of an absorbing barrier for a pOpulation of reproducing entities char- acterized by a diffusion process is given by solving (fi/i)[le exp(va)-v2P exp(le)] + va exp[v1(P-Pb)] - le exp[v2(P-Pb)] + VIP - vZP = 0, (5.64) with _ -1 +«/1 + apt/i3 V1 2 and -1 -\/1 + Apt/P \) = a 2 2 A 2" 2 for P = -§%- and substituting for the diffusion coefficient %- O and the drift coefficient r. The computer program CRITL was used to obtain the family of curves in Figure 5.1 by solving (5.64) for P with a number of values of other parameters. In the program the variables XLM, PB, XMP and X1 correSpond to the quantities I/fi, Pb, fi/P and the assumed value of P (for an iterative solu- tion) reSpectively from (5.64). Since P must be larger than PB for a physically meaningful solution (cf. Section 5.4), the iterative procedure in CRITIC is started with X1 = PB. The function on the left hand side of (5.64) 199 200 is generated by the function subprogram FNC. The root P of (5.64) is computed by using the subroutine CVGl available on magnetic tape on the CDC 6500 computer at Michigan State University. The subroutine works as follows: X1 and X2 are two trial values of a root of (5.64). If both X1 and X2 lie on the same side of the root, the values Y 1 and Y2 of the function FNC evaluated at X1 and X2 will have the same sign. When this is the case, the subroutine sets IND to -l and returns to the main program, to take a new value of X2. This pro- cedure is repeated until X1 and X2 lie on opposite sides of a root. If this procedure leads to too large a value of X2 to be physically meaningful, the iteration is stopped and after printing a message to that effect, the program proceeds to the next set of values of para- meters. Once a value of X2 is found such that the root lies between X1 and X2, the subroutine computes the point where the Straight line connecting the points (X1,Y1) and (X2,Y2) crosses the X-axis. The value of X corresponding to this point is taken as the new value of X2 and the iteration is continued until Y1 or Y2 is sufficiently close to zero. This iterative procedure is illustrated in Figure B.l. Since the subroutine CVGl gives only one root of the function depending on the initial guesses of X1 and X2 and the iteration scheme for changing the values of X2, the reliability of this computer pro- gram hinges on the uniqueness of the solution to the problem of critical length. The shape of the function was studied for a large number of values of the parameters using an IBM 1800 computer coupled with a plotter and it was found that only one root i of (5.64) occurred for P > Pb. Although the uniqueness of the solution to the problem of critical length was not proved theoretically, these results were taken as an indication of the uniqueness. Y = FNC -—-——%F> 201 Figure B.l Iteration scheme for obtaining a root of Equation (5.64) 202 H¢u*.xN.*a n\Ih.~>.H>.H0H HzHao mszHHzoo MN 0» 00 ~0H 0H 00 ~0H HzHao o~.o~.m~..0~IHeo.nH HuozH 0.H+~xu~x o~.-.H~.0zH.nH ~>\Hx*m.~uH¢0 HozH.~x.H00¢.Hx.H0>0 HH<0 .Hx.m>.~>.H>.02nuH00¢ HuozH .H+qu~x ~>uHx .¥.ozxum> .o.0au~> ...zaqu> 0.NH.Hux Hm 00 «no N». Hizhwzwg hH dhmm a IhIa0uIa Ha.~.>.x.uzu onNuzzu APPENDIX C COMPUTATION OF EXTINCTION PROBABILITY FOR A POPULATION 0F REPRODUCING ENTITIES As discussed in Section 5.5, the extinction probability for a p0pulation of reproducing entities characterized by a diffusion process with parameters as in Table 5.2 (given one ancestor initially at y) was computed by solving the backward diffusion equation for the PGF with the initial condition G(9,0\y) = 0. In view of the nonlinear nature of the diffusion equation, the solution was obtained numerically using finite differences. The difference equations for the simulations with conditions described in Table 5.2 are G(9.w+oriiox)-G(ei¢116y) =,gi_{C(QrTl£i+l)GY)-ZG(05T116YI+GLOLTl(i'1)6Y)} 6T 2 6y + r{G(e’TL(i+l)g§)‘GgfllTLiQY)} + ull - G(e.T|16y)] - i1(ioy)G(e,¢\16y)[1-c(e,T\O)j for 6y s i6y s a-oy, (0.1) G®4fi6T10)-G(9JLQI = 02{G(QLT!OY) ‘G(Q:TTLCD4} + r{G(QLTl61)-G(QLTLQ)} 2 5T 6y 6y + uEI-G(e,w\0)] - A1(O)G(0,T\0)[1‘G(0,T\0)] (c.2) and C(e.I+6715)-Gieiwlfib = 02{G(61TL5-6yg-G(9,Ti§)} 6T 6y + r{G(elTl5‘5Y)'G(91TL51} + u[l-G(0,T\5)] - x (a)c(e,w\5)[1-G(e,¢\0>] 6y 1 (c 3) 204 205 when there is a reflecting barrier at y = a, or G = 1 ((2.4) for all T > 0, when there is an absorbing barrier at a. As mentioned above, the initial condition for the difference equations is G(e.0ly) = 0 - (C-S) Equations (C.l) thru (0.5) were solved using program PROB. The notation for PROB is as follows: The arrays G(I) and XG(I) contain values of G(9,T‘i5y) and G(e,T+6T‘i6y) respectively. Similarly, the quantities BH(I), s, R, DH, DT and DX denote x1(i5y), g5, r, u, 67, and 6y respectively in (0.1) thru (C.5). The number M represents the total number of mesh points in the discretizing scheme with Mel equal divisionsof the interval [0,3]. The program PROB given on the following pages characterizes only one set of parameters from Table 5.2. 206 awhm< owthan wm< mw34<> th e0 Hsza .oH.HOH.HuH..H.Im..0o thaa m0 thma 00 thaa .HH.HuH..HOHx..No th10 00 Nsza 00 Nsza 10.00H thaa x0*.0H+.HIH.hxu.H.Hx m HH.NuH m on .ouHHONx monN¢w>m*****u .H0+H mth hq monoma wxb no onhqhaazou***a*u .om*.~010mmu.H.Im 0H .on0*.HIH.uN z.HuH 0H 00 HHarm >x x0xp0u>ho .xoaxovxhouxho omuz mNuxz HIzqu HoHux .ouxh \Ho..N..o..0.H.mflfixhoox0.10.a.m flowomm.* mH whda IhHH00H<2 ..~.mn.xH.HH.x-.Haaa< HIH mH<0a0eeneno .Hqucox mzzHHzoo mm .H00+H0nHH0e.H.0*HH010-..H.0H I.H.nH0n:0+.HHcolHH+H000n>H0nm+l.HIH00+.H00*.NIHH+H.0.nmexH0uHH.0x Hz.NuH mm 00 .H.0+.H .H.0e.H.me100oH0nHH0+..H00IHN00.*>H0*0+W*XH0*HAH00:.N00.*.~u.H_0x .0 u > Ha 20HHH0200 >e<0zaom ozHeomHnma HIH mmHaaoaaoozH HH.0xee*ee0 HH00I.HnHH0 xz.Huo mo 00 .oH.x.HuH..H.0..zH.H0H HzHaa N.Huxo 0N 00 e.an0x HN 00 '2‘"--- ' 207 APPENDIX D ANALYSIS OF EXPERIMENTAL DATA The experimental data on the CLB p0pulations was tested for goodness of fit with a Poisson distribution using the x? test. The program CHSQR was written to do this analysis on an IBM 1800 computer equipped with a keyboard, typewriter, card reader and printer. The logical unit numbers 1, 2, 3 and 6 in the READ and WRITE statements refer to the typewriter, card reader, line printer and keyboard respectively. The function IDTSW takes Signals from the data switches on the main console to go from one part of the program to another. The program CHSQR.calculates the x? Statistics for N replicate populations in Experiments A and B to test the goodness of fit of the data with the Poisson distribution having the same mean as the data. The program was thus used to compute * 2 7- (ni ' n1) x_=§i * (D.1) n . 1. where ni the observed frequency of the number i (or a combina— tion of numbers i1, i2,..., etc.) and n. = the frequency of occurrence of the number i (or the same combination of numbers as in ni) computed using the distribution against which the sequence of numbers is to be tested. 208 209 This value of X? was then compared with that from a standard table of x? distribution for a 10% level of significance and k-2 degrees of freedom, where k is the number of terms in the summation in (D.1). If the computed value of x? was greater than that from the tables, the hypothesis that the sequence of numbers (representing the populations in different replicates) has a Poisson distribution was rejected. For a greater reliability, the numbers in the sequence to be tested are usually combined so that each ni represents one or a combination of more than one numbers in the sequence such that each ni is always greater than or equal to five. This convention was followed for the computations with few exceptions, in which one of the ni's was equal to four. In program CHSQR IX(I,J) contains a two-dimensional array of pOpulations of different life stages on a given day, with I de- noting the serial number of the replicate experiment and J repre- senting the life stages of the CLB. The array IX(I,J) is to be read in as data from cards. The working array X(I) consists of populations of a particular life Stage or a combination of the numbers of individuals in different life stages in the replicate experiments. The program uses Subrouting MOMNT to compute the first four central moments (i.e., mean, variance, skewness and kurtosis, cf. Himmelblau, 1970) of the numbers in X(I). Next, the frequencies of occurrence of various numbers in the sequence X(I) are computed -- IKQ(J) represents the number of times KQ(J) = 0,1,2,..., etc. occurs in the sequence X(I). This as well as the calculation of the (Poisson) probabilities of occurrence POIS(I) of KQ(I) is done using subrouting POISN. To evaluate x2 as in 210 (D.1), different terms from the sequence KQ(I) are to be combined to make each ni greater than or equal to five. This is done by selecting the appropriate numbers from the sequence “0(1) and feeding these as data for evaluation of each ni using the key- board. The individual terms in the summation in (D.l) are computed and printed as CHI(I). The variable CHISQ denotes the final value of x2. Mr 211 onxxz Hz<.z.oxH.ox.mHoa.zz.u.x02wwwflaw“«w u 024 mmHQZm: mam «anna on ozHozoammmaou Az<.z.x.hzzo: HH<0 u >40H .0 szszmnxm mam mama .< szszmaxwaaaaau mom w4mm awrbo >z< «o .meUHHzm oHu zOszOm¢u wIH Hmzpo< O» mmHOszOmam w>Om< 2H mzmwp mo .OZ wIN mH Huaaaau mmqnom HIO mH¢¢xdda Ihai hmwhlmmdsow HID hdwawd ..hHmma¢<4 mHOZH o..buzw30mam thOo< Oh..xm¢..+..h<110m om ..NHovao.. ooz mmq OmZHmzou wtmmh ..hdtmom Om ..>..m.Om.. uz.z<.z.NOH.m00NHmz .mIz.\:Omu:N¢< .H.>*.H.>+x0mux0m 0H onuH 0H O0 .ouzam ..mHmONNOx. Hmeox H<¢NZmu Ihxzom mHOazouuaauau .NIz.\zOmuzwxm .H.N*.H.>+znmn:0m o z.HuH o 00 mmwzzwxm whzaxoo «asauu .HIz.\zOmum<> .H.>+z:mutnm 0 z.HuH 0 00 .ouzOm wuz whaazouuanaao .H.N*.H.~u.H.> N z..H.x onmzwzHo Hz<.z.x.bzzox mthnommnm 216 Ith onthHthHO zommH00 4 «Cu .H.Ox no H+0HusH meOa*.0H.mHOau.sH.mHOa 0xxx<*.0HHmHOau.pH.mHOa ¢m HH.Hu0 an O0 .0000quH meOauH000m~Oa mm 0m.mm.0m..wH.axIH.uH OOH.HuH 0m O0 Nm H+wHupH .xz4I00xmu.aH.mHOa 0N Nmoomon..H.ax.uH HusH .moaxtdlvaxwumeOn zw*****O HHVOXH..HVOX.HHOmHOa.HHHx onmzwzHO Hz<.z.OxH.Ox.mHanzz.zoxOZmHOa wthDOmmDm ozm szhwx .moOHuOHexH.h