Physiologically Structured Population Dynamics: A Modeling Perspective B. W. Kooi F. D. L. Kelpin QUERY SHEET Q1 Q2 Q3 Q4 Q5 Q6 Q7 Au: Au: Au: Au: Au: Au: Au: symbol missing? OK as meant? from where? give location pls. update available? add publisher please Please specify pages 3b2 Version Number File path Date and Time : : : 7.00a/W (Sep 25 2000) P:Santype/Journals/Taylor&Francis/Gctb/8(2-3)/8(2-3)-6627/ctb8(2-3)-6627.3d 5/3/03 and 15:20 Comments on Theoretical Biology, 8: 1–44, 2003 Copyright # 2003 Taylor & Francis 0894-8550/03 $12.00 + .00 DOI: Physiologically Structured Population Dynamics: A Modeling Perspective B. W. Kooi F. D. L. Kelpin Department of Theoretical Biology, Institute of Ecological Science, Faculty of Earth and Life Sciences, Vrije Universiteit, Amsterdam, The Netherlands 5 In classical population dynamic models all the individuals in a population are assumed to be identical or only population averages are considered. The state of the population is therefore the number of individuals or the total amount of biomass in the population. The mathematical model is a first-order ordinary differential equation (ODE) that specifies the time derivative of this population state. Subsequently, the dynamics of an ecosystem where populations interact with each other and with their abiotic environment is usually described by a system of ODEs. The classical models are sometimes called unstructured population models, as opposed to structured population models, which take into account the differences between individuals and characterize an individual by some state. The development of this state is modeled from birth till death, often using a first-order ODE that specifies the time derivative of the individual state. The model is complemented with models for the birth and death of individuals. The structured population model is derived straight forwardly from the individual model using balance laws. The state of a population is no longer a single number; it is the distribution of the individuals over their possible states. If the number of individuals is large, this distribution is continuous instead of discrete. If, in addition, the distribution is a density with respect to the individual state, the population model can be written as a partial differential equation (PDE). The first structured Address correspondence to B. W. Kooi, Department of Theoretical Biology, Institute of Ecological Science, Faculty of Earth and Life Sciences, Vrije Universiteit, De Boelelnaan 1087, 1081 HV Amsterdam, The Netherlands. E-mail: kooi@bio.vu.nl We thank Stijntje Kelpin-Vlaskamp for improving the English. 1 10 15 20 25 2 B. W. Kooi and F. D. L. Kelpin population models use age for the state of the individuals. For many organisms, however, a difference in age alone does not explain the differences in individual behavior and other state variables such as size or energy reserves are more suitable. This has led to physiologically structured population models (PSPs). The life cycle of the individuals may consist of a number of life stages, such as egg, juvenile, and adult. The distribution of the number of individuals with respect to the state variable may be irregular. The introduction of these two extra components, namely, several life stages and nonsmooth distributions, gave mathematical difficulties in the PDE formulation. Recently, an alternative cumulative formulation in terms of renewal integral equations has been proposed that deals with these difficulties. In this review we describe the various model formulations from a modeling perspective. The dynamics of a population of worms that propagate by division into two unequal daughters serves as an example. We give an overview of numerical methods for structured population models. 30 35 40 45 Keywords: age-structured, discrete-time model, continuous-time model, physiologically structured population model, size-structured A population is the collection of individuals belonging to one species living in a spatial region which is either closed or has known immigration and emigration fluxes of individuals across its borders. A description of the dynamics of a population is therefore essentially a description of the dynamics of all of its individuals. There is always individual variability, which is usually incorporated in the model as a random variable. For instance, individuals may be affected differently by environmental stochasticity or mutations, and by random event behavioral aspects, such as motion in space in search of mates or food. As a consequence, also intra- and interspecific encounters, which are important in the predator– prey interplay, are stochastic. If the number of individuals is limited, the individual-based modeling (IBM) (DeAngelis and Gross 1992) approach can be used, where each individual is followed in time. Birth, growth, death but also kinds of behavior, such as defence and mating, can be incorporated in a straight forward manner. Upon birth, the newborn individual is assigned parameter values that are random perturbations of those of its parent(s). These parameters describe the development of the individual. Due to the stochastic effects, each run will give different outcomes. Therefore, the results of one simulation run are of limited interest but several simulations can be done. The final results then are statistics of these runs, for instance, mean values of population statistics, or the extinction probability. 50 55 60 65 Physiologically Structured Population Dynamics 3 The prospects of applying the IBM approach in ecology is in the possibility of applying numerical simulations. In general, no mathematical tools are available to obtain analytical solutions nor to do mathematical analysis. This type of population models can, however, sometimes be analyzed mathematically with the theory of branching processes (Jagers 1975; Arino and Kimmel 1989; Jagers 1991; Arino and Kimmel 1993). This theory was applied to the analysis of models for cell proliferation by fission. Jagers (1991) describes how the classical Galton–Watson branching process, where the population evolves from generation to generation by individuals getting an independent and identically distributed (IID) number of children, yields a dichotomy between extinction or unlimited growth. This dichotomy results from the fact that the population grows, on average, geometrically with a net reproductive rate equal to the mean number of offspring. If the expected number of offspring is one or less than one, the population dies out. If it is larger than one, the population grows beyond all bounds. This means that the population cannot stabilize. If we want stabilization to occur, a description of the reproduction of individuals does not suffice and the interaction with their environment must be included in the model. More elaborate models with different types of individuals can be modeled using multipoint processes. The dynamics of these populations can be studies by individual-based simulations. However, such a simulation is not necessary since the properties of the population can be analyzed by straightforward generalizations of the analysis for the single-type process (see Caswell and John 1992) where a semelparous organism (they die after first reproduction) with fixed generation length is considered. The Galton–Watson branching process models only the branching of the family trees in the population. No information is specified about the times at which the branches are formed. However, the age of the individuals and the ages at child bearing, can be included. Then a mother’s age as well as her type (i-state) determine the starting platform from which the newborn’s life is to evolve. This form of dependence between two next generations leads to a Markov model for the population growth. A superposition of the resulting point processes, essentially Markov renewal theory, provides information such as the average population growth rate (see Jagers, 1991). The life cycle of some species can be divided into a number of stages, thus splitting the population into classes. If the number of individuals in each class can be written as a linear combination of the class sizes some time step ago, the vector containing all class sizes can be multiplied by a so-called population matrix to yield the next population state. The elements in this matrix are transition probabilities from one class into the other, so the matrix is nonnegative. In the special case where the classes are age classes the matrix is called the Leslie matrix and the elements are fecundities on the first row and survivabilities on the subdiagonal. If the matrix is constant, the model is linear and the Perron–Frobenius theorem provides conditions for the 70 75 80 85 90 95 100 105 110 4 B. W. Kooi and F. D. L. Kelpin convergence to a population state that is the eigenvector belonging to the dominant eigenvalue. Once it has reached this state, the population grows geometrically. Each time step, the number of individuals in each class is multiplied with the dominant eigenvalue. This yields a trichotomy of the following three possibilities. If the dominant eigenvalue is less than one, the population goes extinct. If it is larger than one, the population grows without bounds. If it equals one, the population stabilizes. The entries in the matrix can depend on the population size or explicitly on time. Then the matrix equations are nonlinear with possibly positive bounded asymptotic behavior. The mathematical analysis of these matrix models is worked out in Caswell (1989; 1997), Tuljapurkar (1997), and Cushing (1997). Stage-structured population models emphasize the stage structure of the populations. The development of the population is described by a number of life-cycle stages, such as eggs, larvae, pupae, and adults. The rate of change of the stage size must equal the difference between total recruitment and loss rates for the stage. Each term in these equations is related to stagerate processes such as recruitment, death, and maturation, for each stage. The parameters in these equations are stage-dependent mean value rates, for instance, death and birth rates. These balance equations form the dynamic equations for each stage where the number of individuals in the different stages are the population state variables, called p-state variables. The resulting population dynamic model is generally described by delay differential equations (DDEs). This type of model is worked out in Nisbet (1997). Thieme (2003) discusses a stage-structured population model where the definition of a density-dependent transition rate from one stage into the other enables a formulation in terms of ODEs, which are easier to handle than DDEs. In the more general physiologically structured population models, each individual is described by its so-called i-state, which contains a number of physiological traits such as age, energy, or size. If the number of individuals is sufficiently large, a deterministic model can be used (see Metz and de Roos 1992). The individual model is described by individual rates, such as birth, growth, and death rates that describe the change of the momentary i-states of the individual. The deterministic ODE model for the individual together with balance laws directly leads to a PDE model formulation at the population level. In special cases the PDE model can be reduced to an equivalent DDE or an ODE for a p-state variable such as the total number or biomass of all individuals that make up the population. The later case presents unstructured models that average over the i-state distribution in the population without loss of information. It is a special case of a technique referred to as the linear chain trick. Then the model is reduced to a number of coupled ODEs being a closed, finite-dimensional system for moments of the biomass weighted with powers of size (see Metz and Diekmann 1989). 115 120 125 130 135 140 145 150 155 Physiologically Structured Population Dynamics 5 The interplay with the environment is of paramount importance for the techniques that can be used to analyze the resulting model. If the environment is assumed to be constant, the model is linear and mathematical techniques are available to prove existence, uniqueness, and positivity of a solution. As we mentioned earlier, linear models show a trichotomy between extinction, stabilization, and unlimited growth of the population. Stabilization only occurs for specific parameter values. If the population has a time-dependent food source or suffers from a nonconstant predation pressure, the food or the predator has to be modeled explicitly as well. The resulting nonlinear model is autonomous. If the environment changes explicitly with time (seasonal or daily fluctuations) the resulting model is nonautonomous. Generally, in these last two cases the model is nonlinear and consequently more difficult to analyze. Often a middle course is adopted in which the vital rates, such as birth rate, growth rate, and death rate, are density dependent; that is, they depend on the population size (number of individuals or total mass), thus implicitly modeling the effect of the environment on the individual life history. These models are often called semilinear. In nonlinear and semilinear models, biologically realistic positive bounded solutions can occur for large parameter ranges. Age has often been chosen to describe the state of an individual; in addition, the birth and death rate depend solely on age. For continuous-time problems without any interplay with the environment the resulting model is a linear first-order PDE for the density function. Webb (1985) analyzes nonlinear versions of this PDE model. Properties such as the existence and uniqueness of the solution are studied and, with respect to long-term dynamic behavior, existence of a so-called equilibrium or stable age distribution is important. For many biological systems, age is not appropriate to characterize the state of an individual. Other i-state variables such as size, energy, or mass have been proposed leading to a class of models called physiologically structured population models (PSPs), which are formulated and analyzed in Metz and Diekmann (1986) and de Roos (1997). If the environment is constant it is, however, always possible to use age as the individual state variable. Traditionally, the PDE was used. Thieme (1988) works out an example that shows that PDE model formulations can be ill-posed. Diekmann et al. (1998; 2001) propose a cumulative approach where a description of the life cycle of each individual is the starting point for the derivation of the population model formulation. Here ingredients from the theory of multitype branching processes are used. Mathematically this leads to renewal integral equations (RIEs). In the 2001 paper, the nonlinear case is discussed, where interplay with a time-dependent environment is allowed. There is considerable literature on the dynamics of a population consisting of individuals that propagate by binary fissions. Division can be seen as a birth and death at the same time. The mother dies but two or more daughters 160 165 170 175 180 185 190 195 200 6 B. W. Kooi and F. D. L. Kelpin are born at the same time. Division can occur at different i-states. This situation corresponds to a stochastic life length of the cell. Bell and Anderson (1967), Sinko and Streifer (1971), Metz and Diekmann (1986) assume that the number of individuals is large and that there is a population density and analyze the deterministic PDE formulation. Diekmann et al. (1993) assume that the population size is large, but allow for measures and analyze a cell cycle model using the theory based on renewal equations describing the cumulative effect of divisions. Models for the cell cycle of the majority types of cells including bacteria, unicellular eukaryotic organisms, and mammalian cells have been developed and described in Lasota and Mackey (1984) and in Tyson and Hannsgen (1985; 1986). They propose and analyze a model where cells divide into two equal daughter cells. The random division time of the daughters is correlated to that of the mother. In Metz and Diekmann (1986) and Diekmann et al. (1993), the division time of the daughters is also random, but with the same distribution for all cells. In both formulations there is no stable age-distribution when the cells grow exponentially. Lasota and Mackey (1984) discuss cell cycle models in terms of Markov operators. These operators connect the density function of the age or mass distribution at birth over generations. The age-dependent cell generation time case is worked out in Tyrcha (2001). It is assumed that the mass dynamics obeys a first-order autonomic ODE. The generation time is the sum of a random variable plus a constant time after the initiation of the DNA synthesis that triggers the cell to transverse through subsequent phases of the cell cycle that have approximately a constant length. This approach gives alternative proof of the fact that there is no stable mass distribution of the cells that grow exponentially. Many sciences study the link between what is happening at the singleparticle scale and what can be observed among large populations of particles. We mention population balance models (PBMs) which have been used in chemical engineering since the 1960s (see Fredrickson, 1991; Ramkrishna 2000; Fredrickson and Mantzaris 2002; and references therein). In this review we restrict our attention to ecology, in particular to one population that lives in a spatially homogeneous environment. Only physiological traits, such as age, size, and energy reserves, serve as individual state variables. We discuss various model formulations, starting with age-structured models and proceeding to PSP models. The interplay between the choice of the mathematical model formulation and the features of the biological system is emphasized. Generally, no closed-form solutions can be derived for the transient or short-term dynamics and numerical techniques have to be used to approximate the solution. We describe the most popular numerical methods. We elaborate one example in great detail in which the various model formulations as well as their mathematical and numerical analysis are elucidated. We consider the dynamics of a population of worms which 205 210 215 220 225 230 235 240 245 Physiologically Structured Population Dynamics 7 propagate by division into two unequal daughters. We restrict ourselves to the case where there is an abundant food supply. In this type of continuous-time problems it is sometimes advantageous to look for an equivalent discrete-time model formulation. For long-term population dynamics the ratio of the interdivision times of the two sisters is crucial. Convergence to a stable age or size distribution depends on whether this ratio is a rational or irrational number. 250 255 STRUCTURED POPULATION MODELS All individuals in a population share the same species-specific behavior. This behavior generally depends not only on environmental conditions but also on the internal state of the individual. Unstructured population models assume that all the individuals are identical. Structured population models, however, allow for differences between individuals. Each individual has its own internal i-state x. This i-state may vary in time. It necessarily contains all intrinsic information needed to define individual behavior for given environmental conditions. The model for the individual defines the individual’s behavior and the development of its i-state as a function of its environment E and its internal i-state x. In the previous paragraph, the precise meaning of the word ‘‘behavior’’ was deliberately not specified. Whole ranges of models can be thought up to describe one particular population, depending on what is considered to be its relevant behavior. These models differ in their level of detail, but also in the mechanisms that they focus on. There is not one single correct model for a population. The choice of model, and as a result of that also the contents of the i-state, depend on the model’s intended use. First, we discuss age-structured models, which have one single i-state variable, the age a of an individual. They are a subgroup of physiologically structured population models (PSP models), in which the i-state variables represent physiological properties of the individual. 260 265 270 275 Age-Structured Population Models Consider the following age-structured population model. Individuals are born with age 0 and may live up to a maximum age o. The survival function F ðaÞ specifies the probability that an individual survives to age a. So F is a nonincreasing function with F ð0Þ ¼ 1 and F ðoÞ ¼ 0. The maternity function bðaÞ is defined such that bðaÞ da is the number of newborns mothered by an individual from age a to age a þ da. We look at several classical continuous-time and discrete-time formulations of this model, where we partly follow Kot (2001). 280 285 8 B. W. Kooi and F. D. L. Kelpin Continuous-Time Renewal Integral Equation Lotka’s model (Sharpe and Lotka 1911) is a continuous-time model that tracks the population’s birth rate. Let bðtÞ dt be the number of births during time interval t to t þ dt. Then Lotka’s model reads bðtÞ ¼ Z 290 t bðt
aÞF ðaÞbðaÞ da þ GðtÞ ð1Þ 0 where GðtÞ, the contribution of all individuals already present at t ¼ 0, is defined by GðtÞ ¼ Z o
t nð0; aÞ 0 F ða þ tÞ bða þ tÞ da F ðaÞ ð2Þ Here nð0; aÞ equals the age density distribution of the initial population; that is, nð0; aÞ da equals the number of individuals of age a to a þ da at time 0. Notice that GðtÞ ¼ 0 for t > o, in which case the resulting integral equation, Eq. (1), is homogeneous. 295 Discrete-Time Renewal Equation The time t and also age a are now integers. Mathematically the renewal equation yields a difference equation of the form bt ¼ t X bt
a F a ba þ Gt 300 ð3Þ a¼1 where Gt ¼ o
t X a¼1 n0;a F aþt b F a aþt ð4Þ The functions bt and Gt have the same meaning as in the continuous-time case except now at fixed census time intervals. An example of a renewal equation is the Fibonacci (1202) series bt ¼ bt
1 þ bt
2 305 ð5Þ where the number of newborn rabbits bt is the sum of the number of births one (bt
1 ) and two time intervals ago (bt
2 ). In Eq. (3) we have F a ¼ 1 for all a (no mortality) and b1 ¼ 1, b2 ¼ 1, and ba ¼ 0 for a > 2 (one newborn at two preceding time points). 310 Physiologically Structured Population Dynamics 9 Discrete-Time Leslie Matrix Here the state is described by the age distribution in discrete points in time and age. The classical Leslie matrix reads 1 0F n0 ðtÞ 0 P0 B n1 ðtÞ C B B C B B .. C B B C¼B 0 B . C B . B .. A @ . @ . . no
1 ðtÞ 0 0   .. P1 . .. .. . .  0  0 .. . F1 0 0 Po
2 1 Fo
1 0 n0 ðt
1Þ 1 0 C CB n1 ðt
1Þ C C .. CB .. C B . C C B . CB C .. C@ A . 0 A no
1 ðt
1Þ 0 315 ð6Þ We give the Fibonacci series as a example. Then o ¼ 2 and F0 ¼ F1 ¼ P0 ¼ 1, leading to  n0 ðtÞ n1 ðtÞ  ¼  1 1 1 0  n0 ðt
1Þ n1 ðt
1Þ  ð7Þ which is equivalent to the two equations 320 n0 ðtÞ ¼ n0 ðt
1Þ þ n1 ðt
1Þ ð8Þ n1 ðtÞ ¼ n0 ðt
1Þ ð9Þ Substitution of n1 ðtÞ of the second equation in the first equation gives the familiar Eq. (5), now for n0 : n0 ðtÞ ¼ n0 ðt
1Þ þ n0 ðt
2Þ. Indeed, the term n0 ðtÞ is associated with the birth bt rate of the previous paragraph. This is the general procedure to derive from a higher order difference equation a set of first-order difference equations, similar to the procedure in ordinary differential equations. Let NðtÞ denote the total number of individuals at time t. Then we have NðtÞ ¼ n0 ðtÞ þ n1 ðtÞ ¼ n0 ðt þ 1Þ; that is, the number of individuals at time t equals the number of individuals that will be born at t þ 1. Indeed, in this biological interpretation both classes reproduce at the same time and give birth to one child. This simple example shows that the number of individuals in the population is just the summation of the number of individuals born and that survive (all survive in this example) whereby the summation is over all possible instants individuals were born (two in this example). Continuous-Time Partial Differential Equation Let nðt; aÞ denote the population age-distribution function and hence nðt; aÞ da the number of individuals of age a to a þ da at time t. The McKendrick–von Foerster equation (McKendrick, 1926; von Foerster 325 330 335 10 B. W. Kooi and F. D. L. Kelpin 1959) reads 340 @nðt; aÞ @nðt; aÞ þ ¼
dðt; aÞnðt; aÞ @t @a ð10Þ where dðt; aÞ is the per capita mortality rate. The boundary condition at age zero is Z o nðt; 0Þ ¼ nðt; aÞbðt; aÞ da ð11Þ 0 where bðt; aÞ is the per capita birth rate. Observe that we should multiply the left-hand side with da=dt but since it is dimensionless and equals one, it is generally omitted. The initial condition completes the mathematical formulation and reads ð12Þ nð0; aÞ ¼ n0 ðaÞ 345 The resulting model is a first-order linear PDE and mathematical theory exists dealing with the most important issues such as existence and asymptotic age distribution. To that end the problem is formulated in the framework of the theory of semigroups of linear operators (see Webb 1985). In order to derive a link with Lotka’s model [Eq. (1)]; we assume that the birth rate and death rate depend only on age, that is, b ¼ bðaÞ and d ¼ dðaÞ. The survivorship F ðaÞ is linked to the death rate by dF =da ¼
dðaÞ where F ð0Þ ¼ 1. With bðtÞ ¼ nðt; 0Þ we have for t > o, nðt; aÞ ¼ bðt
aÞF ðaÞ. Substitution into Eq. (11) yields the Lotka equation, Eq. (1). Similarly, the equivalence can be shown for t  o. The case of asynchronous exponential growth, where the number of individuals grows or dies out exponentially but their distribution over the i-state space stabilizes, is often observed, for instance, in models where all individuals at all time contribute to reproduction and all individuals interact through one food resource. This can be seen as follows. Suppose there is exponential growth, hence bðtÞ ¼ expfrtg. Substitution into Eq. 1 where t > o gives the characteristic equation, called the Euler–Lotka equation Z o F ðaÞbðaÞ expf
rag da ¼ 1 ð13Þ 350 355 360 365 0 It has exactly one real root r ¼ m, the population growth rate. In Figure 1 the solution nðt; aÞ of Eq. (10) is shown where the initial age distribution, seen on the left above the line t ¼ 0 (the age axis), is uniform. This unnatural initial condition leads to a discontinuity of the density which propagates along the line t ¼ a. For t > o this discontinuity disappears and eventually the distribution converges to the stable age distribution, seen on the Rright. The rates bðaÞ and dðaÞ are such that the population size o NðtÞ ¼ 0 nðt; aÞ da blows-up as is clear from Figure 2 where the size NðtÞ is plotted as function of the time t. 370 375 Physiologically Structured Population Dynamics 11 FIGURE 1 Solution of the linear McKendrick–von Foerster equation, Eqs. (10), (11),  0 0a<a (12), with a ¼ 0:5, o ¼ 1, bðaÞ ¼ , dðaÞ ¼ 1, and n0 ðaÞ ¼ 1. 5 aa1 FIGURE 2 Total number of individuals NðtÞ for the solution of the linear McKendrick–von Foerster equation plotted in Figure 1. 12 B. W. Kooi and F. D. L. Kelpin McKendrick (1926) turned his attention to a single population and considered the linear case, in which the birth and mortality terms were both linear functions of the state variable: the instantaneous distribution of the population. Later on, Gurtin and MacCamy (1974) worked on the nonlinear case, in which the birth rate bðNÞ and the mortality rate dðNÞ depend on the dynamic states of the populations involved. Due to this density dependence there is a stable age distribution and the total number of individual stabilizes. This is illustrated in Figures 3 and 4 where the solutions for the age distribution nðt; aÞ and total number of individuals NðtÞ is shown where death is density dependent similar to the Lotka–Volterra model. 380 385 Reduction to Unstructured Model Formulation In special cases a set of ODEs or DDEs can be constructed which is equivalent to the PDE formulation. Here we give two examples. Reduction to an Ordinary Differential Equation. Suppose that both the birth rate and mortality rate do not depend on age but on the total number of individuals in the population, FIGURE 3 Solution of the nonlinear McKendrick–von Foerster equation, Eqs. (10),  0 0a<a (11), (12), with a ¼ 0:5, o ¼ 1, bðaÞ ¼ , dðt; aÞ ¼ N ðtÞ
1, and 5 aa1 n0 ðaÞ ¼ 1. 390 Physiologically Structured Population Dynamics 13 FIGURE 4 Total number of individuals N ðtÞ for the solution of the nonlinear McKendrick-von Foerster equation plotted in Figure 3. NðtÞ ¼ Z o nðt; aÞ da ð14Þ 0 then the PDE Eq. (10) becomes @nðt; aÞ @nðt; aÞ þ ¼
dðNÞnðt; aÞ @t @a ð15Þ with nðt; 0Þ ¼ Z o nðt; aÞbðNÞ da ¼ bðNÞN and nð0; aÞ ¼ n0 ðaÞ ð16Þ 0 Integration by parts of the second term of Eq. (15), with nðt; oÞ ¼ 0 and the boundary conditions of Eq. (16) gives dN ¼ ½bðNÞ
dðNÞN dt and Nð0Þ ¼ N0 ð17Þ Hence, if the two vital rates that are defined at the individual level, namely, the birth and death rates, do not depend on age, the age structure is superfluous and we can simply use an ODE formulation. It is, however, not possible to recover the age distribution at all times. This is the price we have to pay for dealing with a much simpler model formulation that gives only the time dependence of the total number of individuals. 405 14 B. W. Kooi and F. D. L. Kelpin Reduction to a Delay Differential Equation. Now we assume that there are two life stages, juvenile and adult. The individual mature at age a. The birth rate and mortality rate functions are only positive in the age interval ða; oÞ. Thus, we assume that the juveniles do not die. We perform a reduction procedure similar to the one used in the previous section. There are two agedistribution functions, and their dynamics are described each by a PDE: @nJ ðt; aÞ @nJ ðt; aÞ þ ¼ 0; @t @a @nA ðt; aÞ @nA ðt; aÞ þ ¼
dðNA ÞnA ðt; aÞ @t @a 410 ð18Þ with boundary conditions nJ ðt; 0Þ ¼ bðNA Þ Z o nðt; aÞ da nJ ðt; aÞ ¼ nA ðt; aÞ nA ðt; oÞ ¼ 0 ð19Þ a Ra R o Let NJ ¼ 0 nJ ðt; aÞ da denote the number of juveniles and NA ¼ a nA ðt; aÞ da the number of adults. Integration of Eq. (18) using the boundary conditions of Eq. (19) gives dNA ¼ bðNA ÞNA ðt
aÞ
dðNA ÞNA ðtÞ dt ð20Þ with initial function NA ðtÞ ¼ NA0 for
a  t  0. In a similar way we obtain a DDE for the number of juveniles NJ ðtÞ dNJ ¼ bðNA ÞNA ðtÞ
bðNA ÞNA ðt
aÞ dt 420 ð21Þ As this involves NA ðtÞ, one first has to calculate the number of adults NA ðtÞ and afterward the number of juveniles NJ ðtÞ. This last equation shows that the change in the number of juveniles equals the difference in population birth rate at t (the newborn) minus the population birth rate at t
a (the individuals that mature at t and were born a ago). 425 Physiologically Structured Population Model In Kooijman (2000) it is argued that age is in many cases not an appropriate state variable to describe the growth and development of a population. As an example, consumption of food depends on sizes and the waterflea Daphnia magna typically has a length of 2.5mm at the onset of reproduction. Size differences are also essential in models of cannibalistic populations (Diekmann et al. 1986; van den Bosch et al. 1988; Kirkilionis et al. 2001; Claessen et al. 2000; Claessen and Dieckmann 2002) where large adults eat small juveniles. For bacteria there is a similar situation. The DNA duplication is triggered upon exceeding a fixed cell size and DNA duplication lasts 430 435 Physiologically Structured Population Dynamics 15 a fixed time period (see Donachie, 1968). Hence the division of an individual does not occur at a specific age, but size is also involved. Furthermore, many species go through different life stages. In Kooijman (2000) three stages in the life history of individuals are recognized: egg, juvenile, and adult. Eggs neither eat nor reproduce, juveniles eat but do not reproduce, and adults both eat and reproduce. In the PSP model formulation, age is therefore replaced or augmented with other state variables such as size=mass or energy content. The type of information stored in the i-state varies from model to model, but in many cases x is an n-dimensional vector. Each number in this vector then represents a physiological trait of the individual. Recently, models have been proposed where the i-state is a twodimensional vector. Droop (1973) introduced the concept of cell nutrient quota (internal nutrient density) besides biomass in a model for phytoplankton growth. In this model, growth is decoupled from nutrient uptake. Phytoplankton cells store nutrients and the growth rate of a cell depends on the amount of stored nutrient. Hallam et al. (1990) describe a model in which an organism consists of storage: lipids, the major source of energy, and a structural component—nonlipid dry material consisting mostly of proteins and carbohydrates. In Persson et al. (1998) the components are called the reversible and irreversible mass. In reversible mass they include energy reserves such as fat, muscle tissue, and gonads, and in irreversible mass they include compounds like bones and organs, which cannot be starved away by the consumer. The Dynamic Energy Budget (DEB) model (Kooijman 2000), deals with storage materials (reserves) and structural biomass. Here the storage materials are carbohydrates and lipids but also a part of ribosomal RNA and of the proteins (see Hanegraaf et al. 2000). The storage materials have a dual function as a reserve pool for energy and elementary compounds. These storage materials are continuously used and replenished, while structural materials are permanent and require maintenance. The reserves adapt to the external food availability, that is, to the environment. If this food availability is constant or varies slowly compared to the adaptation process, the energy reserves will be in or near their equilibrium state. Individuals will not be spread out all over the i-state space, but concentrated on or near the line through these equilibrium values for energy reserves. Given a PSP model, the set of all possible i-states is called the i-state space O. If the individual goes through different life stages, O may be divided into subsets Oi (e.g., OE ; OJ ; OA for the egg, juvenile, and adult stage). The boundaries of these subsets mark points where the individuals move from one life stage to the next, such as where they hatch or mature. As we saw in the age-structured model formulation, the p-state of the population can be a number or mass density that depends on the i-state variables. This leads to a PDE formulation for the population with boundary and initial conditions similar to the age-structured case described in the previous section. 440 445 450 455 460 465 470 475 480 16 B. W. Kooi and F. D. L. Kelpin In Diekmann et al. (1998) the mathematical difficulties that may arise in this formulation are explained in biological terms as follows. In order to keep models parameter scarce, one wants to allow for discontinuities with respect to i-state of rates (think of waterflies that start to reproduce upon reaching a critical size; see Thieme 1988). Now consider a situation in which the i-state of some individual moves in O for an extended period of time along a line of discontinuity of the rate of offspring production, that is, the line that separates OJ and OA . Then the model is not acceptable as a model and one should not expect that existence and uniqueness of solutions at the p-level holds. Whether this phenomenon actually occurs or not in a specific model, is hidden in the rates. It is the combined global effect of the rates that makes the difference between the model being ill- or well-posed. In this section we first introduce the model ingredients. Then we show how in one big step the behavior of all the individuals in a population is lifted to the population level. We start with the classical formulation of PSP models (Metz and Diekmann 1986), which uses densities and leads to a partial differential equation with boundary and initial conditions. Then we continue with the more recently developed cumulative formulation leading to a renewal integral equation. 485 490 495 500 i-State Formulation If the individuals go through life stages, the population is split into subpopulations that contain the individuals in each life stage. We now describe the dynamics for one single population or subpopulation. Later on, boundary conditions will connect the subpopulations at the borders of their respective subspaces Oi . In the PDE formulation, an i-state space O (or Oi ) always is a subset of Rn , for some n. The i-states in different life stages may have different interpretations and their dimension may vary from one stage to another. Typical instantaneous vital rates that describe the temporal change of the i-state variable are the growth, birth, and death rates, which will occur as per capita rates at the population level. Growth Rate. Since we assume that the environment is homogeneous in space, the spatial coordinates of an individual are irrelevant for its behavior. The coordinates of its i-state x, however, follow a trajectory through O. The growth rate gðE; xÞ 2 Rn defines a vector field on O that determines the speed and direction of this movement through i-state space. This implies that individual growth is deterministic. Note that entries of g do not refer to the physiological ‘‘growth’’ rate if the corresponding i-state is not size. If gi represents growth it is sometimes allowed to have gi ðE; xÞ < 0 >, that is, a shrinking individual. This can lead to hard mathematical problems and therefore it is generally assumed that gi ðE; xÞ > 0; that is, the direction of the movement is never reversed and the movement never stops. 505 510 515 520 Physiologically Structured Population Dynamics 17 Death Rate. Death removes an individual from the population. Death is modeled by specifying the death rate dðE; xÞ  0. The environmental condition that influences the death rate generally is predation pressure. Emigration or harvest of individuals from the region under consideration can also be modeled using the death rate. 525 Birth Rate. Different reproduction modes exist. Some of these can be modeled using an individual birth rate bðE; xb ; xÞ  0. This is the rate at which an individual with state x produces offspring with state at birth xb . Often the support of this function is some OA  O, the subset of adults. 530 Influence on the Environment=Feeding=Predation. The environment EðtÞ can represent food consumed by the population. When food is abundant the growth rate does not depend on the food (gðxÞ). On the other hand, when food is limiting, growth depends on the food availability (gðEðtÞ; xÞ). In the latter case the amount of food available depends on the population size since this determines how much food is consumed. This feedback interplay is typical for predator–prey, or similar trophic, interactions in ecological communities. In a similar way the population can be consumed by another population in which case the death rate stands for natural death rate plus predation rate (dðEðtÞ; xÞ). If food is assumed to be an abiotic nutrient which is supplied to the region where the population lives in, the dynamics of this nutrient is generally described by an ODE. In case of interplay with other populations they can be described by one or more ODEs or even PDE. The growth process is new with respect to age-structured population models. That is, when gðEðtÞ; xÞ ¼ 1, where x is one-dimensional, we get essentially an age-structured population model. As mentioned earlier, the environment sometimes implicitly modeled as density-dependent birth and death rates. Murphy (1983) studies a model with a density-dependent growth rate. 535 540 545 550 p-State Formulation The state of the population (the p-state) is now represented by a continuous number density function nðxÞ. The integral of this density function over an area of i-states equals the total number of individuals in that area. To be a R bit more precise, for any measurable set o  O, o nðxÞ dx equals the number of individuals in the population with i-states in o. To express the fact that this p-state changes in time, we will denote it as nðt; xÞ, where x is a vector if the i-state space is multidimensional. Field Equation. The p-state can be shown (see for instance Metz and Diekmann, 1986, 14), to satisfy the PDE 555 560 18 B. W. Kooi and F. D. L. Kelpin or shorter @nðt; xÞ X @ þ ðgðEðtÞ; xÞnðt; xÞÞ ¼
dðEðtÞ; xÞnðt; xÞ @t @xi i ð22Þ 565 @nðt; xÞ þ H  ðgðEðtÞ; xÞnðt; xÞÞ ¼
dðEðtÞ; xÞnðt; xÞ @t ð23Þ This basically is a transport equation with velocity g, but with a nonzero right-hand side, which represents the death of individuals. Boundary Conditions. Individuals may be born on @O, the border of O. Only part of this border, @ þ O, is of interest. This is because the growth vector gðE; xÞ should point inward. Otherwise newborn individuals will leave O, which they are not supposed to do. In mathematical terms: Birth only occurs in i-states xb 2 @ þ O, where the inner product of g with the inward-pointing normal n is positive. Then the boundary condition reads nðxb Þ  gðEðtÞ; xb Þnðt; xb Þ ¼ Z bðEðtÞ; xb ; xÞ nðt; xÞ dx xb 2 @ þ O 570 ð24Þ O This is sometimes also called the renewal equation. Initial Conditions. Initial conditions complete the mathematical formulation: nð0; xÞ ¼ n0 ðxÞ x 2 O0 ð25Þ where O0  O is the support of the density n. In Murphy (1983) the PDE and its boundary and initial conditions are transformed from numbers density into biomass density. The resulting PDE formulation is the same as Eqs. (23), (24), and (25), but the interpretations of some coefficients differ. We remark that in general this transformation is not needed. If the dynamics of the environment, for instance, nutrients or predators, are described by an ODE, the initial amount of these substances completes the initial conditions and that ODE is solved together with the PDE for the population. Just as in density-dependent age-structured models, also in PSP models the birth, growth, and death rates can be assumed to be density dependent in order to avoid modeling the food or predator explicitly. For populations that propagate by binary fission birth rate and death rate (division rate) are coupled. If division is into two equal daughters the resulting PDE reads: @nðt; a; vÞ @nðt; a; vÞ @gðt; a; vÞnðt; a; vÞ þ þ @t @a @v ¼
½bðt; a; vÞ þ dðt; a; vÞnðt; a; vÞ þ 2½2bðt; a; 2vÞnðt; a; 2vÞ ð26Þ 580 585 590 Physiologically Structured Population Dynamics 19 where bðt; a; vÞ is now the per capita division rate and dðt; a; vÞ the death rate, that is, here the disappearance of individuals by causes other than division. The terms of the left-hand side are the same as before, but those of the right-hand side are sink and source terms associated with division. If an individual with size v divides, the first term, two newborns appear at size v from an individual with size 2v. Characteristics. The formulation in Eqs. (23), (24), and (25), is based on the Eulerian approach where a description is given of the flow of individuals around a fixed point in the i-state space O 2 Rn . The Lagrangian approach follows a frame moving with the individuals along their trajectories in the istate space O. These trajectories correspond to the so-called characteristics of the PDE. Characteristics are curves in the space O  Rþ . The following equation defines the characteristic, parameterized by a (which can be age but is not necessarily so), that is, tðaÞ and xðaÞ. dt ¼1 da dx ¼ gðEðtÞ; xÞ 2 Rn da 600 605 ð27Þ The PDE Eq. (23) can be rewritten as @nðt; xÞ þ Hnðt; xÞ  gðEðtÞ; xÞ ¼
nðt; xÞ½H  gðEðtÞ; xÞ þ dðEðtÞ; xÞ ð28Þ @t where we used H  ½gðEðtÞ; xÞnðt; xÞ ¼ Hnðt; xÞ  gðEðtÞ; xÞ þ nðt; xÞ½H  gðEðtÞ; xÞ ð29Þ Along the characteristic of the PDE Eq. (23), that is, tðaÞ and xðaÞ satisfy Eq. (27), we obtain with the chain rule the result 615 dn @nðt; xÞ dt dx ¼ þ Hnðt; xÞ  ¼
nðt; xÞ½H  gðEðtÞ; xÞ þ dðEðtÞ; xÞ ð30Þ da @t da da The initial conditions tð0Þ; xð0Þ are now formulated. Initially we have tð0Þ ¼ 0 and xð0Þ ¼ x, substituted in Eq. (25). For t > 0 the initial values are related to the appearance rate of the newborn individuals, which is fixed by the boundary condition for the PDE [Eq. (16)] where tð0Þ ¼ t and xð0Þ ¼ xb 2 Oþ . Loosely speaking, the PDE problem with boundary and initial conditions is rewritten as an initial value problem for a set of infinite number of ODEs. In the age-structured case we have gðxÞ ¼ 1 and the solution of Eq. (27) gives characteristic curves that are straight lines in the ða; tÞ plane with slope 1. If the growth rate gðxÞ depends on the i-state x and not on time (for instance, in the case of constant food availability), the characteristic curves 620 625 20 B. W. Kooi and F. D. L. Kelpin tðaÞ and xðaÞ can be calculated by solving the ODE Eq. (27) without knowing the density distribution nðt; xÞ. This is in the linear case. However, in the nonlinear case where gðEðtÞ; xÞ depends on time or on the state of the population itself with density dependent i-state rates, the characteristic curves and the density nðtðaÞ; xðaÞÞ have to be calculated simultaneously. If newborns only start at the boundary xb 2 @ þ O it is possible to transform a size-structured population density nðt; xÞ to an age-structured density mðt; aÞ. This technique is known as the Murphy trick (Murphy 1983). Cumulative and RIE Formulation. A more general, so-called ‘‘cumulative’’ formulation exists in terms of renewal integral equations (Diekmann et al., 1998; 2001). Here the ingredients that serve to describe the processes at the i-level are not rates, such as the birth rate, but quantities at the generation level where the whole life history of an individual is taken into account, such as the lifetime expected number of offspring. The RIE formulation can be split into different steps. First, the environmental conditions are assumed to be known. With these conditions as an input, individual development, reproduction, and survival can be modeled in a manner that resembles multitype branching processes. From this description for one individual, an expression is derived for its entire ‘‘clan.’’ This clan contains an individual plus all its descendants. The definition of the dynamics of the population now comes down to summing up all clans of the individuals in an initial population. Finally, the population influences its environment. This output follows from the population dynamics. It should of course equal the input that was initially assumed. Only the fixed-point problem of closing the feedback loop between the environment and the population is nonlinear. This explains why the cumulative formulation is successfully used for the continuation of equilibria, where the environment is constant (see Kirkilionis et al. 2001 and Diekmann et al. 2003). i-State Development and Survival. Diekmann et al. (2001) give the cumulative formulation for the populations from the previous subsection, which we copy here for reference. The input from the environment is a function I that is defined on a time interval ½0; ‘ðIÞ. So ‘ðIÞ is the length of the input. Two operations on inputs are of importance here, the restriction operator rðsÞ, which restricts the input to the time interval ½0; s, and a shift yð
sÞ which is defined by ðyð
sÞIÞðtÞ ¼ Iðt þ sÞ; 0  t < ‘ðIÞ
s. This means the input is shifted s forward in time. For a given environmental input I, one can construct the trajectory of an individual with initial i-state xb . Diekmann et al. (2001) define XI ðx0 Þ; as the i-state of an individual at time ‘ðIÞ; given that it had i-state x0 at time zero, it experienced input I and it survived. For populations with instantaneous growth and death rates, the function t 7! XrðtÞI ðx0 Þ is the solution of the initial value ODE problem 630 635 640 645 650 655 660 665 670 Physiologically Structured Population Dynamics d xðtÞ ¼ gðIðtÞ; xðtÞÞ xð0Þ ¼ x0 dt 21 ð31aÞ The survival probability F I ðx0 Þ of this individual equals F I ðx0 Þ ¼ expf
Z ‘ðIÞ dðIðtÞ; XrðsÞI ðx0 ÞÞ dsg ð32Þ 0 The development measure uI ðxb ; oÞ combines this information. It specifies the probability that, after experiencing input I , an individual with initial state xb is still alive and has an i-state in o  O.  F I ðx0 Þ if xI ðx0 Þ 2 o 0 otherwise ¼ dxI ðx0 Þ ðoÞF I ðx0 Þ uI ðxb ; oÞ ¼ ð33Þ ð34Þ where d denotes a Dirac delta measure. In general, though, u need not be confined to one single point. Stochastic movement through i-state space can therefore be modeled using this more general cumulative formulation. This is not possible using the classical PDE, which requires an instantaneous growth rate g and therefore assumes that growth is deterministic. Reproduction Kernel. The lifetime reproduction measure LI ðxb ; oÞ equals the expected number of offspring with states of birth in o. The formulation is the same whether individuals are born in internal points of O or on the boundary Oþ LI ðxb ; oÞ ¼ Z Z o 680 685 690 ‘ðIÞ bðIðtÞ; xb ; XrðtÞI ðx0 ÞÞ F rðtÞI ðx0 Þ dxb dt ð35Þ 0 Then a next generation operator is defined, which is the expected number of offspring but now for each distribution of individuals. This leads to a generation expansion, that is, the iteration of the reproduction rule to specify the expected total offspring of the entire clan. In a similar way a next state operator is defined that is the distribution of the individuals for a specific distribution of individuals at time zero. For the definition and mathematical treatment of these two operators the reader is referred to (Diekmann et al. (2001). Combining the two ingredients, namely, the reproduction together with i-state development and survival, specifies the structured population model by adding the contributions of all individuals. 695 700 22 B. W. Kooi and F. D. L. Kelpin NUMERICAL TECHNIQUES FOR PHYSIOLOGICALLY STRUCTURED POPULATION MODELS 705 Extensive descriptions of numerical methods for the time integration of structured population models are given in the literature. Finite difference techniques are derived from the mathematical model by replacing derivatives by differential quotients. With other techniques biological interpretation is taken into account with developing the numerical procedure. In full discretization schemes, time and i-state space are discretized simultaneously. In classical finite difference schemes (Richtmeyer and Morton 1967), derivatives are replaced by differential quotients based on Taylor series expansion in grid points. For the nonlinear density-dependent models the classical Lax–Wendroff method is used in Sulsky (1993) for the agestructured model and in Sulsky (1994) for size-structured models. In Sulsky (1993; 1994) a fixed grid is used and the resulting scheme is an adapted version of the classical Lax–Wendroff method. For age-structured models the support for the density distributions is known for the initial conditions and subsequently for all times in the time interval of interest. For sizestructured populations, the support of the size distribution is not known beforehand and even several magnitudes in the size of the individuals with in the population may occur. Also, the support of the size density may increase with time. This leads to wasted work over much of the grid at early times in the calculation when this density function is zero in a large portion of the computational domain. In the box method (Fairweather and López-Marcos 1991; Angulo and López-Marcos 2002), the starting point is the balance law for a cell that is between x and x þ Dx, and t and t þ Dt: Z xþDx nðt þ Dt; xÞ
nðt; xÞ dx 710 715 720 725 x þ Z tþDt Z tþDt gðEðtÞ; x þ DxÞnðt; x þ DxÞ
gðEðtÞ; xÞnðt; xÞ dt t ¼
t Z xþDx dðEðtÞ; xÞnðt; xÞ dx dt ð36Þ x The difference scheme is obtained by discretization of the integral terms in this equation. Another implicit scheme is proposed and analyzed in Ackleh and Ito (1997). In semidiscretization schemes, initially only the i-state space is discretized. In Gurney and Nisbet (1998) a number of these schemes (upwind difference scheme and central difference scheme) are discussed. The resulting set of ODEs is then further discretized in time using an ODE solver. Similar to the full discretization scheme, these methods work with a discrete representation of the density function. 735 Physiologically Structured Population Dynamics 23 Perhaps the simplest scheme for the age-structured model formulation, which is a full discretization scheme as well as a semidiscretization scheme, is an explicit finite difference scheme on a regular grid. As the grid size in both time and age direction are taken equal the method is also equivalent to a time-integration along the characteristics with equidistant time steps. At each time step the number of newborns with age zero is obtained as a sum of adults weighted with the reproduction rate. This simple scheme was used with the calculation of the results presented in Figures 1–4. For age- and size-structured PDE models a similar approximation scheme is proposed in Ito et al. (1991), Milner and Rabbiolo (1992), Angulo and López-Marcos (2000, 2002), and Kostova (2002). This procedure comprises three steps. First a grid is constructed such that each point belongs to the same characteristic. This is done by numerical integration of the individual growth equation, Eq. (31). In principle any ODE solver can be used. Solving the ODEs of Eq. (30) with initial points on the characteristic curve gives the solution for the next time step. Again, any ODE solver can be used. The final step is to calculate the population birth rate occurring in the boundary condition, which is obtained by employing a numerical quadrature for the integral term in Eq. (16). If the environment is constant the characteristic curves are fixed in time; that is, all individuals born at points in @ þ O follow trajectories that originate in these points. If, on the other hand, the birth and death rates are density dependent, a similar approach can be used where the last two steps are more complicated. Kostova (2002) gives a short overview of existing numerical methods with proven global error estimates and describes an explicit method of third order. In Angulo and López-Marcos (2002) the efficiency of three numerical methods—the box method, the Lax–Wendroff, and a characteristics scheme—are assessed by considering five problems with diverse degrees of freedom. They conclude that for an autonomous problem no best method exists and that all specific features involved in the numerical integration have to be considered in order to choose a specific numerical scheme for a particular problem. For nonautonomous problems the box method was more efficient than the Lax–Wendroff scheme in four of the five test problems. These schemes are related to the Eulerian approach. The escalator boxcar train (EBT) method (de Roos 1988) uses the Lagrangian approach. This method can be used for a broad class of PSP models. The EBT method (de Roos 1988) uses a semidiscretization scheme. In order to discretize the i-state space, individuals with similar i-states are grouped into so-called cohorts. Each cohort is represented by a delta peak placed in the average i-state of the individuals in the cohort. The height of the peak equals the number of individuals in the cohort. The p-state is approximated by the sum of these delta peaks. Figure 5 gives an example of the discretization of a one-dimensional i-state space where the p-state is represented by a continuous density 740 745 750 755 760 765 770 775 780 24 B. W. Kooi and F. D. L. Kelpin FIGURE 5 An i-state discretization of a continuous density function nð0; xÞ into n cohorts. Individuals with i-states in the sameRrange Oj form a cohort. This cohort is represented by a delta peak with height lj ¼ Oj nð0; xÞ dx (the number of individuals R in the cohort), placed at Zj ¼ l1j Oj xnð0; xÞ dx (the average i-state of all individuals in the cohort). function. The method also works, however, for a higher dimensional i-state space and for p-states given by a measure on the i-state space. Each individual follows a trajectory through the i-state space, which is a solution of the ODE that describes the individual development. The righthand side of the ODE is a function of the individual’s i-state and the environment. Initially, cohorts consist of individuals with nearly identical i-states. If the solutions of the development ODE do not diverge, the individuals in a cohort will stay close in the time interval of interest. As a consequence, the trajectory starting in the initial average i-state of a cohort will at all times be a good approximation for the i-states of the individuals in that cohort. The death rate on this trajectory will also be a good approximation for the death rates of all individuals in the cohort, as death rates too depend only on i-state and environment. A treatment of classical type boundary conditions with continuous reproduction by all adult individuals can be found in de Roos (1988). The introduction of one or several special boundary cohorts preserves the order of the numerical integration technique. These boundary cohorts are filled with 785 790 795 800 Physiologically Structured Population Dynamics 25 newborn individuals. Once they are large enough they are transformed into ordinary cohorts. Due to the continuous introduction of new cohorts, the total amount of cohorts may grow so large that the numerical solution of the system is slowed down significantly. If this happens, we can reduce the number of cohorts by merging cohorts with nearly the same i-state or by removing cohorts with a sufficiently small number of individuals. Since all finite-difference methods are based on Taylor-series expansion, it is assumed that the frequency distribution as well as the birth and death rate distribution functions are sufficiently smooth. For a model with pulsed reproduction, that is, a model where individuals are born in a specific period of the year or when a certain threshold for the mothers is satisfied, we successfully used the simple forward Euler scheme in Kelpin et al. (2000). We showed that for a nontrivial test problem, an age-structured population model with pulsed reproduction, the calculated results obtained with a first-order forward Euler scheme agree well with those obtained with the EBT method. In Ackleh and Deng (2000) a (monotone) approximation, based on an upper and lower solutions technique, is developed for the nonautonomous size-structured model. 805 810 815 820 WORMS Naidids occur in a wide range of aquatic habitats such as rivers, lakes and estuaries, sewage filterbeds, slow sand filters of water works, and drinkingwater distribution systems. In sewage treatment plants the number of worms present is inversely related to sludge disposal; however, their presence is still unpredictable and uncontrollable (see Ratsak, 2001). Naidids are often important detrivores in organic polluted river systems and therefore may be an important element in the transfer of energy from the heterotrophic bacteria associated with organic matter to higher trophic levels such as fish. The length of Nais elinguis is generally less than 12 mm and the mean diameter is 0.15 mm. It is a hermaphrodite, under normal conditions it propagates by division into an anterior and a somewhat smaller posterior part. It was observed that both the anterior and the posterior worm divide when they reach the same particular volume. Individual growth of two sisters with constant, abundant food availability is shown in Figure 6. The naidids were collected from a sewage treatment plant and fed with activated sludge particles from the same plant. The anterior worm initially grows faster than her posterior sister but the difference diminishes and is negligible when they divide. Ratsak et al. (1993) give a model for the feeding, growth, and division of worms, based on DEB theory (Kooijman 2000). The model gives an explanation for the initial growth spurt of anterior worms after division. While food flows from mouth to anus its nutritional value diminishes as nutrients are 825 830 835 840 26 B. W. Kooi and F. D. L. Kelpin FIGURE 6 Growth curves of sister naidids. The upper and lower curves represent the fit of the DEB model to the measured data of anterior () and posterior naidid (u), respectively; see Ratsak et al. (1993). Observe that the worms grow in length only, so length is directly proportional to volume. The length at division is vd ¼ 8:41 mm. assimilated into the energy reserves of the worm. These energy reserves are used for maintenance and growth. Close to the mouth, the food still contains many nutrients and the energy reserves are richer. After division, the anterior worm therefore has richer energy reserves and initially grows faster than her posterior sister. We first give the mathematical model for the individual, in dimensionless form. Then we show how a series of time-scale assumptions greatly reduce the complexity of this model. 845 850 Individual Model Formulation Naturally, the i-state of a worm in this model contains its length or volume, but we use, according to the (DEB) model (Kooijman 2000), its volume v and its energy density e. Ratsak et al. (1993) parameterize the worm along its length using a parameter x 2 ½0; 1, which equals 1 at the mouth and 0 at the anus. A function eðt; xÞ is the energy reserve density at time t and at position x. A severe drawback is that the i-state becomes infinite dimensional, since it includes the function eð; xÞ: ½0; 1 ! R. The population dynamics can no longer be described in the familiar PDE formulation. The ingestion through the mouth of the worm equals Im f ðtÞv, where Im is the maximum ingestion rate and f ðtÞ is the scaled functional response. A more rigorous version of the model would also model the digestion of food inside the gut. Here we simply assume that the nutritional value of the food in the gut of the worm decreases linearly from 1 at the mouth (x ¼ 1) to 0 at the anus (x ¼ 0). Then the energy reserve density at location x is assumed to follow 855 860 865 Physiologically Structured Population Dynamics @eðt; xÞ ¼ nð2f ðtÞx
eðt; xÞÞ @t 27 ð37Þ The growth rate of a worm is computed by integrating the local growth rate along its length: dv ¼v dt Q1 Z 0 1 neðt; xÞ
mg dx eðt; xÞ þ g 870 ð38Þ When a worm’s volume reaches v ¼ vd , it divides. The volumes of the two daughters are a fixed fraction of that of the mother, not necessarily one half. The anterior daughter has volume va ¼ ð1
xd Þvd ; the posterior daughter has volume vp ¼ xd vd ; see also Figure 6. So just prior to the division the range of the anterior part is from 1 to xd , and of the posterior part from xd to 0. At division we redefine the variable x, such that it still has the value 1 at the anterior end and 0 at the posterior end of the daughters. Thus, for the anterior part position x of the daughter corresponds with position xð1
xd Þ þ xd of the mother. For the posterior part, position x of the daughter corresponds with x xd of the mother. Observe that the distribution of the nutritional value for both anterior and posterior worm changes abruptly at division. Let us assume that the energy density is at equilibrium at the moment of division eðtb ; xÞ ¼ 2f ðtb Þx. Then the initial condition for the anterior and posterior daughters, respectively, is ea ðtb ; xÞ ¼ 2f ðtb Þ½xd þ x ð1
xd Þ ð39Þ ep ðtb ; xÞ ¼ 2f ðtb Þxd x ð40Þ 875 880 885 and Hence, all individuals belonging to the same time have an unique state at birth, that is the same volume vi and energy density distribution ei ðtb ; xÞ, i 2 fa; pg. That is, the mother decides when a daughter is born but passes not her i-state information. The growth rate does not depend on the x explicitly, only on the integral over a function that depends on the energy density distribution over x. For all x we have Z t n½2f ðtÞx
ei ðt; xÞ dt ð41Þ ei ðt; xÞ ¼ ei ðtb ; xÞ þ 890 tb Substitution of this expression in Eq. (38) gives the instantaneous growth rate, which depends on the history of the environment f ðtÞ and the energy reserves since birth. In the next sections we describe various model formulations for the population dynamics of the worm where mortality is assumed to be zero. We 900 28 B. W. Kooi and F. D. L. Kelpin assume that food supply is abundant, that is, the environment is timeindependent and f ¼ 1. This allows for an age-dependent model formulation. We start with a continuous-time formulation and we continue with a discrete-time formulation. Continuous-time Age-Dependent Model Formulation 905 We assume that there are two-populations with densities ma ðt; aÞ defined on ½0; 1Þ  ½0; ta  and mp ðt; aÞ on ½0; 1Þ  ½0; tp  with interdivision times ta and tp . The PDEs read @ @ m ðt; aÞ ¼
m ðt; aÞ @t a @a a @ @ m ðt; aÞ ¼
m ðt; aÞ @t p @a p ð42Þ We arrive at a boundary condition for this coupled PDEs system, Eq. (42), at a ¼ aa ¼ ta and a ¼ ap ¼ tp : ma ðt; 0Þ ¼ ma ðt; ta Þ þ mp ðt; tp Þ ð43Þ mp ðt; 0Þ ¼ mp ðt; tp Þ þ ma ðt; ta Þ which describe division. The interdivision times ta and tp are implicitly given by ta 1 nea ðt; xÞ
mg dx dt ea ðt; xÞ þ g 0 0 Z tp Z 1 nep ðt; xÞ
mg dx dt lnfvd g ¼ lnfvp g þ ep ðt; xÞ þ g 0 0 lnfvd g ¼ lnfva g þ Z Z ð44aÞ ð44bÞ where ea ðt; xÞ and ep ðt; xÞ are given by Eq. (39) and Eq. (40) with constant food f ¼ 1. In this continuous-time model we introduce the number of individuals in the populations: Na ðtÞ ¼ Z ta ma ðt; aÞ da and Np ðtÞ ¼ 0 Z tp mp ðt; aÞ da ð45Þ 0 Integration by parts of Eq. (42) and the use of the boundary conditions Eq. (43) yields dNa ¼ mp ðt; tp Þ dt and dNp ¼ ma ðt; ta Þ dt ð46Þ 910 Physiologically Structured Population Dynamics 29 These results have a clear biological interpretation; when an anterior individual divides an extra new posterior worm is born, and vice versa, when an posterior individual divides an extra new anterior worm is born. Let T ðtÞ and Ai ðaÞ be defined such that mi ðt; aÞ ¼ T ðtÞAi ðaÞ. Then separation of variables yields that the two functions are determined by dT j ¼ cj T j dt T j ð0Þ ¼ 1 925 ð47Þ and 930 dAaj da ¼
cj Aaj and dApj da ¼
cj Apj ð48Þ with boundary conditions Aaj ð0Þ ¼ Apj ð0Þ ¼ Aaj ðta Þ þ Apj ðtp Þ ð49Þ This gives T j ¼ expfcj tg and Aaj ðaÞ ¼ K expf
cj ag and Apj ðaÞ ¼ K expf
cj ag ð50Þ where K is a normalization factor. Substitution of these results into Eq. (49) yields the Euler–Lotka equation expf
cj ta g þ expf
cj tp g ¼ 1 ð51Þ The roots cj give the following solutions of the PDE, Eq. (42): ma ðt; aÞ ¼ mj expfcj tg expf
cj ag and mp ðt; aÞ ¼ mj expfcj tg expf
cj ag, where mj is a constant fixed by the initial conditions. The principle of superposition of these solutions can, at least formally, be applied to describe the transient solution of the linear, PDEs Eq. (42). Then the Fourier coefficients (the mj terms) are fixed by the initial age distribution, which is assumed to be sufficiently smooth. In practice it turned out that for a reasonable approximation many terms are needed and this restricts applicability for practical purposes. This completes the analysis of the short-term dynamics. We assume that the age distributions for both populations are the stable distributions belonging to the real root c0 of the Euler–Lotka equation, Eq. (51), Aa0 ðaÞ and Ap0 ðaÞ. Then both populations grow exponentially, so T ðtÞ ¼ expfc0 tg or, if Nw ðtÞ ¼ Na ðtÞ þ Np ðtÞ denotes the total number of worms, dNa ¼ c0 N a dt and dNw ¼ c0 Nw dt ð52Þ 940 945 950 30 B. W. Kooi and F. D. L. Kelpin With Eq. (46) and Eq. (51) we get c0 Na ¼ T 0 ðtÞAp0 ðtp Þ 955 and c0 Nw ¼ T 0 ðtÞAp0 ð0Þ ð53Þ Then, with Eq. (50) and Eq. (51) we obtain the result  t =t Na ðtÞ Na ðtÞ a p þ ¼1 Nw ðtÞ Nw ðtÞ ð54Þ Under the condition that for j 6¼ 0, Re cj < c0 , Eq. (54) holds true asymptotically for t ! 1. We have Re cj ¼ c0 if and only if the ratio ta =tp is rational, that is, there are multiple poles denoted by j with Re cj ¼ c0 , j 6¼ 0. We elaborate this case later with the discrete-time formalism. The case where ta =tp is irrational is dealt with in the next section where a renewal equation for the cumulative birth rate is derived. Continuous-Time RIE Model Formulation 965 El Houssif (2001) derives a renewal equation for the cumulative birth rate. The Renewal Theorem is used to derive an expression for the asymptotic behavior of the cumulative birth rate and for the asymptotic composition of the population. He obtains Eq. (54) but now for irrational ratios ta =tp . As in the last part of the previous section, so is food assumed to be constant. In El Houssif (2001), mortality was taken into account. For the sake of simplicity we assume mortality zero. Let BðtÞ denote the cumulative number of individuals that passed through size va during the time interval ½0; tÞ. Then the renewal equation reads for t > tp , that is, when all individuals belonging to the initial population have divided at least once, BðtÞ ¼ Bðt
ta Þ þ Bðt
tp Þ þ GðtÞ ð55Þ The term GðtÞ is the contribution of the population at time t ¼ 0 to BðtÞ. It is the sum of three contributions given by GðtÞ ¼ Z va H½t
tp þ ta þ aw ðvÞnp ð0; vÞ ðdvÞ vp þ Z vd  va 960  Hðt
tp þ aw ðvÞÞ þ Hðt
2tp þ ta þ aw ðvÞÞ nw ð0; vÞ ðdvÞ ð56Þ where H is the (Heaviside) step function defined by HðxÞ ¼ 0 for x < 0 and HðxÞ ¼ 1 for x > 0. The population size and composition at time 0 is described by a Borel measure nw ð0; vÞ ¼ np ð0; vÞ for vp < v < va and 970 975 Physiologically Structured Population Dynamics 31 nw ð0; vÞ ¼ na ð0; vÞ þ np ð0; vÞ for va < v < vd . The function aw ðvÞ is the time it takes for an individual to grow from vp to v. Thus tp ¼ aw ðvd Þ and ta ¼ aw ðvd Þ
aw ðva Þ. The first term is the contribution of all posterior worms with size vp < v > va , the second term all those who are anterior worms after division, and the third term those who are posterior worms after division. The explicit expression for the solution of Eq. (55) reads BðtÞ ¼ GðtÞ þ k¼n   1 X X n n¼1 k¼0 k Gðt
kðtp
ta Þ
nta Þ Z 990 ð57Þ where GðtÞ ¼ 0 for t < 0. The renewal theorem (Feller 1971 and Jagers 1975, section 5.2) is used to deduce the asymptotic behavior for BðtÞ for t ! 1. Let GðsÞ be a real function then the Laplace–Stieltjes transform ^f ðaÞ is defined as ^ ðaÞ ¼ G 985 995 1 expf
asgG ðdsÞ ð58Þ 0 The characteristic equation Eq. (51) is now defined using this Laplace– Stieltjes transform technique, ^0 ðaÞ ¼ expf
ata g þ expf
atp g ¼ 1 m ð59Þ There exists a unique real root a1 that is positive. Two cases are distinguished. In the lattice case the ratio of the interdivision times, ta =tp , is rational, that is ta =tp ¼ l=k with l, k 2 N where for the greatest common divisor (gcd) of l and k we have gcdðl; kÞ ¼ 1. In the nonlattice case this ratio is irrational. For both cases El Houssif (2001) derives the asymptotic behavior of BðtÞ where t ! 1. He formulates expressions for Na ðtÞ and Np ðtÞ in terms of BðtÞ and thereafter he derives the asymptotic behavior for Na ðtÞ and Np ðtÞ given in Eq. (54) for the nonlattice case whether the death rate is positive or zero, but in the lattice case only when the death rate is zero, and it is conjectured that the relation holds also for the lattice case when the death rate is positive. In the next section we make a link with a discrete-time formulation where we derive the asymptotic behavior given in Eq. (54) in the lattice case where death rate is zero. Discrete-Time Model Formulation We assume that ta ¼ l=ktp for k ¼ 1; 2; . . . ; 1 and l ¼ 1; 2; . . . ; k
1, while gcdðk; lÞ ¼ 1, hence the lattice case defined in the previous section. 1005 1010 1015 32 B. W. Kooi and F. D. L. Kelpin The number of anterior daughters after the i-th division is denoted by Na ðita Þ, the number of posterior ones by Np ðita Þ, and the total number of individuals by Nw ðita Þ. 1020 Discrete-Time Renewal Equation The number of anterior daughters after the i-th division follows the generalized Fibonacci series Na ðita Þ ¼ Na ðði
lÞta Þ þ Na ðði
kÞta Þ for i  k ð60Þ with, for instance, one of Na ð0Þ, Na ðta Þ; . . . ; Na ððk
1Þta Þ equal to 1. This follows from the following two relationships: Np ðita Þ ¼ Np ðði
lÞta Þ þ Na ðði
lÞta Þ ð61Þ Na ðita Þ ¼ Na ðði
kÞta Þ þ Np ðði
kÞta Þ ð62Þ 1025 and This implies for the number of posterior daughters Np ðita Þ ¼ Na ðði þ k
lÞta Þ 1030 ð63Þ The Fibonacci series, Eq. (60), converges to a geometrical one with Na ðita Þ ¼ aNa ðði
1Þta Þ, with a given by the characteristic equation ak
ak
l
1 ¼ 0 or 1 1
¼1 ak al ð64Þ With c0 ¼ ta
1 ln a Eq. (64) is equivalent to Eq. (51) being the characteristic equation for the continuous-time case. In Kooi and Boer (1995) this approach is worked out in detail; here we continue with a population matrix formulation. 1035 Discrete-Time Leslie Matrix We continue with a description of an approach using life cycle graphs and Leslie matrices (see Caswell 1989). We define life stages here as follows. At time il ta the length interval ½1=kvd ; vd  is divided into k subintervals each representing a life stage. The number of individuals at time il ta in life stage j with volume j=kvd  v  ðj þ 1Þ=kvd is denoted as njp ðil ta Þ. The anterior worms live only in life stages k
l þ 1  j  k. The duration of one time step is now tla . The number of anterior and posterior worms in the same life stage are lumped together with nkw ðil ta Þ ¼ nkp ðil ta Þ þ nka ðil ta Þ as shown in Figure 7 where the life-cycle graph (see Caswell 1989) for a one-population and two-population formalism is given. 1040 1045 Physiologically Structured Population Dynamics 33 FIGURE 7 Life-cycle graphs for the one- and two-population formalism. As an example of the two-population formalism we consider the special case where l ¼ 1 and k ¼ 2. The interdivision time of the posterior worm tp is a multiple (two times) of the interdivision time of the anterior worm ta . This problem is directly related to the Fibonacci series by the following population matrix: 1 0 0 n1p ðita Þ C @ B 2 @ np ðita Þ A ¼ 1 0 n2a ðita Þ 0 1 0 1 1 10 1 1 np ðði
1Þta Þ C B 0 A@ n2p ðði
1Þta Þ A 1 n2a ðði
1Þta Þ 1050 ð65Þ Now Na ðita Þ ¼ n2a ðita Þ and Np ðita Þ ¼ n1p ðita Þ þ n2p ðita Þ and this leads directly to the Fibonacci series, Eq. (60). The asymptotic behavior can be obtained here straightforwardly by the study of the population matrix in Eq. (7) where n0 ðita Þ ¼ n2a ðita Þ þ n2p ðita Þ and n1 ðita Þ ¼ n1p ðita Þ. In the one-population formalism the population projection matrix for the general case reads 1055 1060 34 B. W. Kooi and F. D. L. Kelpin 1 0 i nlþ1 0 w ðl ta Þ C B1 B . C B B .. C B B nk ð i t Þ C B 0 B w l a C B . C¼B B 1 i B B np ðl ta Þ C B .. C B . .. C B B A @ .. @ . 0 nlp ðil ta Þ 0  0 .. . .. . .. . 1 0    0 .. .. .. . . . .. .. .. . . . .. . 1 0  0 1  10 lþ1 i
1 1 nw ð l ta Þ 1 C B .. 0C C B . C B .. C C C k i
1 . CB n ð t Þ C B w a l .. C C B 1 i
1 t Þ n ð C B a p .C l CB C CB . C .. 0 A@ A l i
1 0 np ð l t a Þ ð66Þ The nonzero elements in the population matrix are two elements of the first row that model the division and the subdiagonal, which models survival to the next class. The Perron–Frobenius theorem applies (see Caswell 1989, 58, 60), because the graph is irreducible (i.e., there is a path in the graph from every node to every node) and the graph is primitive (i.e., it is irreducible and gcdðl; kÞ ¼ 1). The characteristic equation for this matrix equation is the same as Eq. (64). The column eigenvector corresponding to the eigenvalue a reads ak
1 ; ak
2 ; . . . ; a; 1ðTÞ k
1 X ai1 ¼ i¼k
l k
1 X Np ¼ i¼0 ai1 k ak
l
1 1
a1 ¼ 1
a1 1
a1 1
ak1
ak
l 1 ¼ 1
a1 1
a1 1070 ð67Þ The existence of a real eigenvalue larger than 1 follows from the evaluation of the function ð1
al Þak
l . So, there is one real dominant eigenvalue a1 > 1, a root of Eq. (64) with largest real part. This multiplier a1 relates to the asymptotic population growth rate m as m ¼ ta
1 ln a1 . Kooi and Boer (1995) show that for the population projection matrix given in Eq. (66) there is a real eigenvalue with greater absolute value than any other eigenvalue. The simpler case l ¼ 1, which is the interdivision time of the posterior worm tp , is a multiple of the interdivision time of the anterior worm ta , dealt with in Ratsak et al. (1993). Frauenthal (1986) proves, using spectral decomposition, that if one starts with a strictly positive vector and repeatedly applies the linear transformation defined in Eq. (66), the resulting vector converges to the column eigenvector, Eq. (67), corresponding to the dominant real eigenvalue. The relationship Eq. (54) is obtained with Na ¼ 1065 and ð68Þ where a1 is again the positive real dominant eigenvalue. In combination with the characteristic Eq. (64) for the dominant eigenvalue a1 > 1, this implies 1075 1080 1085 35 Physiologically Structured Population Dynamics Na 1 ¼ k
l Np a1 or Na 1 1 ¼ ¼ k k
l Nw ð1 þ a1 Þ a1 ð69Þ Consequently, we have the relationship Eq. (54) as time goes to infinity. 1090 Comparison of the Discrete-Time and the Continuous-Time Model For worms living in abundant food environments that divide into two possibly unequal daughters, the anterior and posterior ones with interdivision times ta and tp , the asymptotic population growth rate m is given by Eq. (51) with m ¼ c0 . This is also the characteristic equation Eq. (64), with m ¼ ta
1 ln a. Hence this equation holds whether the ratio of the interdivision times ta and tp is rational or irrational. Asymptotically for the fraction of anterior worms, that is, Na divided by the total number of worms Nw , relationship Eq. (54) holds. The results obtained with the discrete-time model can be used to clarify results of the continuous-time model. Consider the case of division in a fixed partition, where ta =tp ¼ 7=19. In Figure 8 we show the trajectories of the worms when we start out with one cohort of worms. The scaled age of the posterior individual a=tp is shown as a function of that of the sister, the anterior individual, a=ta . That is, 7 revolutions along the center line of a torus representing the growth of the anterior individual are accompanied by FIGURE 8 The increase in the scaled age of the posterior worm, a=tp , as a function of the scaled age of the anterior worm, a=ta . Dotted line for rational ratio of interdivision times ta =tp ¼ 7=19 and solid line for irrational ratio ta =tp ¼ 1
ðln 2Þ=ðln 3Þ. 1095 1100 1105 36 B. W. Kooi and F. D. L. Kelpin 19 revolutions along the torus. Then the line is continued along the path travelled at time 19ta ¼ 7tp ago, showing periodicity. We conclude that the age distribution does not converge to a stable age distribution. The asymptotic behavior is cyclic. When the ratio of the interdivision times is rational, both the continuoustime and the discrete-time model apply. We lump the number of individuals into the life stages introduced earlier, so at time il ta in life stage j we have, R j=kvd np ðil ta ; vÞ dv. The relationship between the growth rate njp ðil ta Þ ¼ ðj
1Þ=kv d m ¼ c0 and the dominant real eigenvalue a1 is: expfmta g ¼ a1 . Convergence to the eigenvector of Eq. (67) means that the distributions of the total numbers in the life stage converge. So, within a life stage there is no convergence, and this is known from the solution of the continuous-time model, but the distribution of the total numbers of individuals within the life stages converges to an asymptotic distribution. If the ratio is irrational the dimension of the matrix becomes infinite. Therefore we study a series ln =kn with kn ! 1 for n ! 1 such that this series of rational numbers converges to this irrational number. Then the length of the life stages converges to zero, and therefore we have effectively convergence to a length distribution. However, the rate of convergence becomes small because the absolute magnitude of the eigenvalues of the second largest one is almost equal to the largest real eigenvalue. Hence, there is convergence, but the rate becomes infinitely small. In Figure 8 we show the trajectories also for the case ta =tp ¼ 1
ln 2=ðln 3Þ and thus irrational. We do not have periodicity now as we had before in the rational case (where ta =tp ¼ 7=19 in Figure 8). 1110 1115 1120 1125 1130 CONCLUDING REMARKS The following observation can be made. First unstructured population models were replaced by simple structured ones where age served as istate variable. Mathematical techniques were developed to prove existence, permanence, and stability properties of the ultimate age or size distribution. Then the model for the individual underlying the structured population model became more complex for the need to get more biological realism. Later other i-state variables then age were used and growth depending on the food availability was taken into account, as well as different life stages. Because food is consumed by the population but on the other hand growth depends on the food availability, generally the food dynamics has to be modeled as well. The individual model was described by individual rates in deterministic ODEs, and balance laws directly lead to a PDE model formulation at the population level. The mathematical techniques were adapted to cope with these more complex models and these techniques seemed initially to be successful. It was possible to take food dynamics into account and to formulate 1135 1140 1145 Physiologically Structured Population Dynamics Q2 37 the interaction via a feedback. However, rather subtle problems arose and very different mathematical techniques had to be developed to tackle these problems. The step for i-state to p-state is now postponed by dealing with cumulative quantities at the generation level. Stochastic movement through i-state space can therefore be modeled. This approach is partly based on techniques from branching processes. We have mentioned that the growth rate is generally assumed strictly positive, gðE; xÞ > 0. Then the direction of the movement along the characteristics neither reverses nor stops. Growth occurs when available internal energy reserves or external food is sufficient to meet maintenance costs which are needed to maintain the state of living. Consequently, the individual dies when gðE; xÞ ¼ 0; see Kooijman (2000), Muller and Nisbet (2000), and Hanegraaf and Kooi (2002). In Kooi and Kooijman (1999) we studied the dynamics of a rotifer population consuming algae in a chemostat. During a transient there was insufficient food to support life of the whole population for a relatively short period. When an unstructured model formulation is used where all individuals are identical, death of one member would directly lead to extinction of the whole population. In Kooi and Kooijman (1999) we showed that when food becomes scarce for a relatively short period, some unlucky individuals with insufficient energy reserves will die, leaving more food for survivors. Later when food recovered the population recovers as well and continuously persists. In Gyllenberg et al. (1997) and Persson et al. (1998), populations with both continuous, namely, death, and discrete, namely, reproduction, elements are investigated. In Gyllenberg et al. (1997) the dynamics of the population with pulsed reproduction are formulated by a difference equation for one state variable, namely, the population size. The influence of the environment is taken into account as a density-dependent within-season mortality. This formulation is appropriate for simple ecological systems of seasonally breeding populations with non-overlapping generations. In the worms case study, the age at which reproduction occurs determines the projection time interval in the Leslie matrix formulation. In Kooi et al. (2001) a similar model formulation was used to analyze the asymptotic dynamics of a waterflea Daphnia magna population with iteroparous propagation (they reproduce more than once and die immediately after last reproduction). In that model all individuals are born with a species specific biovolume and are clones of the parent. In the fixed period between consecutive reproductive events the adult density-dependent mortality is a continuous process. Natural death occurs when adults reach the species specific maximum age immediately after their last reproductive event. With the Leslie matrix model formulation it was possible to clarify longterm behavior results obtained with numerical integration of the continuoustime model. Cushing (1998) derives the continuous-time PDE model from the discretetime matrix model by the limiting process of infinitesimally refining the 1150 1155 1160 1165 1170 1175 1180 1185 1190 38 B. W. Kooi and F. D. L. Kelpin classes and the census time interval. This approach shows in the case of agestructured populations the relationship between the McKendrick–von Foerster formulation and the Leslie matrix formulation. The relationships between the various discrete- and continuous-time models for the dynamics of PSP models are elucidated in de Roos et al. (1992), where a numerical scheme for the solution of physiologically structured populations is derived from the Leslie matrix model. In that formulation there is one life stage and therefore the time step is equal to the fixed time period between two consecutive reproduction events. Most models presented and analyzed in the literature deal with one isolated structured population or one in interaction with unstructured resources or predators. The modeling of the trophic interactions is a very important issue when we want to deal with multispecies systems such as communities and biosystems, where in the latter case the interaction with a physical environment is explicitly modeled. In principle, each population can be described by a PDE together with boundary and initial conditions. For several agestructured species, discrete- and continuous-time models are proposed in Cushing (1998). With size-structured populations, modeling interactions among species, predator=prey as well as competition interactions are even more troublesome. Numerical difficulties with simulations of spatially inhomogeneous populations with physiologically distinguishable individuals were discussed in Gurney et al. (2001). In the worms case, the division into an anterior and a posterior worm disturbs the equilibrium state of the energy reserves, as the anterior daughter gets a bigger part of the reserves than the posterior one. It is assumed that the energy dynamics are fast enough to smooth out the disturbance caused by one division before either of the daughter worms divides. In other words, both daughter worms reach the equilibrium value for energy reserves before they in turn divide. This means that all worms have the same i-state when they divide. All anterior worms are born identical and therefore the population of anterior worms is concentrated on one line in i-state space. The same is true for the posterior worms, so there are only two possible states of birth. In community ecology, predator=prey interaction, host–parasitoid interaction, inducible defense, competition, mutualism, symbiosis, nutrient recycling, invasion, succession, persistence, stability, and so on play an important role in theory building. These theoretical concepts are used in the study of ecological issues such as the stability–diversity of ecosystems debate, the maximum length of food chains, and occurrence of succession. Unstructured ODE population models have been used frequently to address these issues whereby the short- and long-term dynamic behavior is studied. Spatially structured population models (Durrett and Levin 1994, Gurney et al. 2001), in which the i-state contains the spatial coordinates of an individual, are not discussed in this article but are important in many applications. Structured population models can also be combined with adaptive dynamics 1195 1200 1205 1210 1215 1220 1225 1230 1235 Physiologically Structured Population Dynamics 39 theory to study evolutionary processes. In Claessen and Dieckmann (2002) an age-structured model for predators is used and resource polymorphisms and adaptive speciation are studied. For general structured metapopulation models, fitness is defined in Gyllenberg and Metz (2001). The theory of structured populations can be applied to metapopulations straightforwardly. There is analogy between local populations or dispersers and individuals, and between metapopulation and population. Birth is now colonization of empty patches by dispersers and death is extinction of a patch population (Gyllenberg and Metz 2001). Thieme (2003) links age-structured models with endemic disease models. There are a limited number of papers dealing with the bifurcation analyses of (PSP) models. We mention the work of Cushing (1989) and Kirkilionis et al. (2001). A computer package BASE (Kirkilionis et al. 2001) is available for continuation of steady states of PSP models. For finite dimensional maps and ODE systems, computer packages are available. We mention AUTO (Doedel et al. 1997) and CONTENT (Kuznetsov and Levitin 1997; Kuznetsov 1998). These packages are in these days heavily used for the numerical bifurcation analysis of unstructured ecosystem (food chain and food web) models. Although these packages can be used for continuation of bifurcation points of PSP models in principle, this is far from easy in practise. The availability of a bifurcation package for PSP would stimulate potential users to apply these models for populations as part of food chain or food web ecosystems. It would be interesting to have the efficiency of the three numerical methods—the box method, the Lax–Wendroff, and a characteristics scheme—be compared with that of the EBT method and those based on upper and lower solutions techniques. Especially when the support of the density function is time dependent, the EBT method works fine, while methods using fixed grids in the i-space are less efficient in those situations because this grid has to span the largest support in the time interval of interest. In the case where the population density vanishes almost everywhere in the i-state space, integration along the characteristics and the EBT algorithm follows cohorts of individuals through the i-state space and will work fine. Another possibility is to reduce the dimension of the i-state space and to transform to a size or an age-structured population with normal population densities using the so-called Murphy trick (Murphy 1983). Sometimes mathematical techniques can be used to reduce the PSP model to equivalent unstructured ODEs or to stage-structured model DDEs (see Nisbet 1997). In Kooi and Kooijman (1999) and Bocharov and Hadeler (2000), conservation laws play an important role in the derivation of the reduced model. When this applies, assumptions often made implicitly during the construction of unstructured models are made explicit and the parameters of the reduced model can be expressed as lumped parameters of the structured model. 1240 1245 1250 1255 1260 1265 1270 1275 1280 40 B. W. Kooi and F. D. L. Kelpin The aggregation technique is an analytical tool which has been used for unstructured models and is at present applied to structured population models. Arino et al. (2000) study a model of an age-structured population with two time scales. The slow scale corresponds to the demographic process and the fast one describes the migration process between different spatial patches. Singular perturbation technique is used to show that the solutions of the system can be approximated by means of the solutions of a simpler (unstructured) problem. In some papers PDE formulations were already used in case studies addressing the issues mentioned above, but in general this is a tremendous task. In many studies the PSP model formulation has proven to be essential for realistic population modeling; we mention the study of cannibalism. Mathematical and=or biological insight is needed to know beforehand under which conditions the use of a structured instead of an unstructured population model is essential to have realism that is required for the purpose for which the mathematical model is used. 1285 1290 1295 REFERENCES Ackleh, A. S., and K. Deng. 2000. A monotone approximation for a nonautonomous sizestructured population model. Q. Appl. Math. 58(2–3):103–113. Ackleh, A. S., and K. Ito. 1997. An implicit finite difference scheme for the nonlinear sizestructured population model. Numer. Funct. Anal. Optimiz. 18(9–10):865–884. Angulo, O., and J. C. López-Marcos. 2000. Numerical integration of nonlinear size-structured population equations. Ecol. Model. 133:3–14. Angulo, O., and J. C. López-Marcos. 2002. Numerical integration of autonomous and nonautonomous non-linear size-structured population models. Math. Biosci. 177:30–71. Arino, O., and M. Kimmel. 1989. Asymptotic behavior of a nonlinear functional-integral equation of cell kinetics with unequal division. J. Math. Biol. 27:341–354. Arino, O., and M. Kimmel. 1993. Comparison of approaches to modelling of cell population dynamics. SIAM J. Appl. Math. 53(5):1480–1504. Arino, O., E. Sanchez, R. Bravo de la Parra, and P. Auger. 2000. A singular perturbation in an age-structured population model. SIAM J. Appl. Math. 60(2):408–436. Bell, G. I., and E. C. Anderson. 1967. Cell growth and division I. A mathematical model with applications to cell volume distributions in mammalian suspension cultures. Biophys. J. 7:431–444. Bocharov, G., and K. P. Hadeler. 2000. Structured population models, conservation laws, and delay equations. J. Diff. Eq. 168:212–237. Caswell, H. 1989. Matrix population models. Sunderland, MA: Sinauer Associates. Caswell, H. 1997. Matrix methods for population analysis. In Structured-population models in marine, terrestrial, and freshwater systems, ed. S. Tuljapurkar and H. Caswell, 19–58. New York: Chapman & Hall. Caswell, H., and A. M. John. 1992. From the individual to the population in demographic models. In Individual based approaches and models in ecology, ed. D. L. DeAngelis and L. J. Gross, 36–61. Springer-Verlag. Claessen, D., A. M. de Roos, and L. Persson. 2000. Dwarfs and giants: Cannibalism and competition in size-structured populations. Am. Naturalist 155(2):219–237. 1300 1305 1310 1315 1320 1325 Physiologically Structured Population Dynamics Q3 Q4 41 Claessen, D., and U. Dieckmann. 2002. Ontogenetic niche shifts and evolutionary branching in size-structured populations. Evol. Ecol. Res. 4:189–217. Cushing, J. M. 1997. Nonlinear matrix equations and population dynamics. In Structuredpopulation models in marine, terrestrial, and freshwater systems, ed. S. Tuljapurkar and H. Caswell, 205–243. New York: Chapman & Hall. Cushing, J. M. 1998. An introduction to structured population dynamics, vol. 71. Philadelphia: Society for Industrial and Applied Mathematics. DeAngelis, D. L., and L. J. Gross. 1992. Individuals-based models and approaches in ecology. New York: Chapman & Hall. de Roos, A. M. 1988. Numerical methods for structured population models: the escalator boxcar train. Num. Meth. Part. Diff. Eq. 4:173–195. de Roos, A. M. 1997. A gentle introduction to physiologically structured population models. In: Structured-population models in marine, terrestrial, and freshwater systems, ed. S. Tuljapurkar and H. Caswell, 119–204. New York: Chapman & Hall. de Roos, A. M., O. Diekmann, and J. A. J. Metz. 1992. Studying the dynamics of structured population models: A versatile technique and its application to Daphnia. Am. Naturalist 139(1):123–147. Diekmann, O., M. Gyllenberg, H. Huang, M. Kirkilionis, J. A. J. Metz, and H. R. Thieme. 2001. On the formulation and analysis of general deterministic structured population models. II. Nonlinear theory. J. Math. Biol. 43(2):157–189. Diekmann, O., M. Gyllenberg, and J. A. J. Metz. 2003. Steady state analysis of structured population models. Technical report, preprint. Diekmann, O., M. Gyllenberg, J. A. J. Metz, and H. R. Thieme. 1998. On the formulation and analysis of general deterministic structured population models. I. Linear theory. J. Math. Biol. 36(4):349–388. Diekmann, O., M. Gyllenberg, H. R. Thieme, and S. M. Verduyn Lunel. 1993. A cell-cycle model revisited. Technical report, CWI-Report MAM-R9305, Amsterdam, the Netherlands. Diekmann, O., R. M. Nisbet, W. S. C. Gurney, and F. van den Bosch. 1986. Simple mathematical models for cannibalism: a critique and a new approach. Math. Biosci. 78:21–46. Doedel, E. J., A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, and X. Wang. 1997. Auto 97: Continuation and bifurcation software for ordinary differential equations. Technical report, Concordia University, Montreal, Canada. Donachie, W. D. 1968. Relationship between cell size and time of initiation of DNA replication. Nature 219:1077–1079. Droop, M. R. 1973. Some thoughts in nutrient limitation in algae. J. Phycol. 9:264–272. Durrett, R., and S. Levin. 1994. The importance of being discrete (and spatial). Theor. Pop. Biol. 46:363–394. El Houssif, N. E. 2001. A deterministic size-structured population model for the worm Naidis elinguis. J. Math. Biol. 42:424–438. Fairweather, G. and J. C. López-Marcos. 1991. A box method for a nonlinear equation of population dynamics. IMA J. Num. Anal. 11:525–538. Feller, W. 1971. An introduction to probability theory and its application, 2nd ed. New York: Wiley. Frauenthal, J. C. 1986. Analysis of age-structured models. In Mathematical ecology, ed. T. A. Hallam, and S. A. Levin, 117–147. Springer-Verlag. Fredrickson, A. G. 1991. Segregated, structured, distributed models and their role in microbial ecology: A case study based on work done on the filter-feeding ciliate Tetrahymena pyriformis. Microb. Ecol. 22:139–159. Fredrickson, A. G., and N. V. Mantzaris. 2002. A new set of population balance equations for microbial and cell cultures. Chem. Eng. Sci. 57(12):2265–2278. 1330 1335 1340 1345 1350 1355 1360 1365 1370 1375 42 Q5 B. W. Kooi and F. D. L. Kelpin Gurney, W. S. C., and R. M. Nisbet. 1998. Ecological dynamics. New York: Oxford University Press. Gurney, W. S. C., D. C. Speirs, S. N. Wood, E. D. Clarke, and M. R. Heath. 2001. Simulating spatially and physiologically structured populations. J. An. Ecol. 70(6):881–894. Gurtin, M. E., and R. C. MacCamy. 1974. Nonlinear age-dependent population dynamics. Arch. Rat. Mech. Anal. 54:281–300. Gyllenberg, M., I. Hanski, and T. Lindström. 1997. Continuous versus discrete single species population models with adjustable reproductive strategies. Bull. Math. Biol. 59:679–705. Gyllenberg, M., and J. A. J. Metz. 2001. On fitness in structured metapopulations. J. Math. Biol. 43(6):545–560. Hallam, T. G., R. R. Lassiter, J. Li, and L. A. Suarez. 1990. Modelling individuals employing an integrated energy response: Application to Daphnia. Ecology 71:938–954. Hanegraaf, P. P. F., and B. W. Kooi. 2002. The dynamics of a tri-trophic food chain with twocomponent populations from a biochemical perspective. Ecol. Model. 152(1):47–64. Hanegraaf, P. P. F., B. W. Kooi, and S. A. L. M. Kooijman. 2000. The role of intracellular components in food chain dynamics. C. R. Acad. Sci. Ser. III Sci. Vie–Life 323:99–111. Ito, K., F. Kappel, and G. Peichl. 1991. A fully discretized approximation scheme for sizestructured population models. SIAM J. Numer. Anal. 28(4):923–954. Jagers, P. 1975. Branching processes with biological applications. Probability and mathematical statistical–applied probability and statistics. New York: Wiley. Jagers, P. 1991. The growth and stabilization of populations. Stat. Sci. 6:269–283. Kelpin, F. D. L., M. A. Kirkilionis, and B. W. Kooi. 2000. Numerical methods and parameter estimation of a structured population model with discrete events in the life history. J. Theor. Biol. 207(2):217–230. Kirkilionis, M. A., O. Diekmann, B. Lisser, M. Nool, A. M. de Roos, and B. P. Sommeijer. 2001. Numerical continuation of equilibria of physiologically structured population models I. Theory. Math. Models Methods Appl. Sci. 11(6):1101–1127. Kooi, B. W., and M. P. Boer. 1995. Discrete and continuous time population models, a comparison concerning proliferation by fission. J. Biol. Syst. 3(2):543–558. Kooi, B. W., T. G. Hallam, F. D. L. Kelpin, C. M. Krohn, and S. A. L. M. Kooijman. 2001. Iteroparous reproduction strategies and population dynamics. Bull. Math. Biol. 63(4):769– 794. Kooi, B. W., and S. A. L. M. Kooijman. 1999. Discrete event versus continuous approach to reproduction in structured population dynamics. Theor. Pop. Biol. 56(1):91–105. Kooijman, S. A. L. M. 2000. Dynamic energy and mass budgets in biological systems. Cambridge: Cambridge University Press. Kostova, T. 2002. An explicit third-order numerical method for size-structured population equations. Numer. Methods Partial Differ. Eq. in press. Kot, M. 2001. Elements of mathematical ecology. Cambridge: Cambridge University Press. Kuznetsov, Y. A. 1998. Elements of applied bifurcation theory, vol. 112 of Applied mathematical sciences, 2nd ed. New York: Springer-Verlag. Kuznetsov, Y. A., and V. V. Levitin. 1997. CONTENT: Integrated environment for the analysis of dynamical systems. Centrum voor Wiskunde en Informatica (CWI), Amsterdam, the Netherlands, 1.5 edition. Lasota, A., and M. C. Mackey. 1984. Globally asymptotic properties of proliferating cell populations. J. Math. Biol. 19:43–62. McKendrick, A. G. 1926. Application of mathematics to medical problems. Proc. Edinburgh Math. Soc. 44:98–130. 1380 1385 1390 1395 1400 1405 1410 1415 1420 1425 Physiologically Structured Population Dynamics Q6 43 Metz, J. A. J., and A. M. de Roos. 1992. The role of physiologically structured population models within a general individual-based modeling perspective. In Individuals-based models and approaches in ecology, ed. D. L. DeAngelis and L. J. Gross, 88–111. New York: Chapman & Hall. Metz, J. A. J., and O. Diekmann. 1986. The dynamics of physiologically structured populations, vol. 68 of Lecture notes in biomathematics. Berlin: Springer-Verlag. Metz, J. A. J., and O. Diekmann. 1989. Exact finite dimensional representations of models for physiologically structured populations. I: The abstract foundations of linear chain trickery. In Differential equations with applications in biology, physics and engineering, ed. J. A. Goldstein, F. Kappel, and W. Schappacher, 269–289. New York: Marcel Dekker. Milner, F. A., and G. Rabbiolo. 1992. Rapidly converging numerical algorithms for models of population dynamics. J. Math. Biol. 30(7):733–754. Muller, E. B., and R. M. Nisbet. 2000. Survival and production in variable resource environments. Bull. Math. Biol. 62:1163–1189. Murphy, L. F. 1983. A nonlinear growth mechanism in size structured population dynamics. J. Theor. Biol. 104:493–506. Nisbet, R. M. 1997. Delay-differential equations for structured populations. In Structuredpopulation models in marine, terrestrial, and freshwater systems, ed. S. Tuljapurkar and H. Caswell, 89–118. New York: Chapman & Hall. Persson, L., K. Leonardsson, A. M. de Roos, M. Gyllenberg, and B. Christensen. 1998. Ontogenetic scaling of foraging rates and the dynamics of a size-structured consumer-resource model. Theor. Pop. Biol. 54:270–293. Ramkrishna, D. 2000. Population balances—Theory and applications to particulate systems in engineering. San Diego: Academic Press. Ratsak, C. H. 2001. Effects of Nais elinguis on the performance of an activated sludge plant. Hydrobiologia 463:217–222. Ratsak, C. H., S. A. L. M. Kooijman, and B. W. Kooi. 1993. Modeling the growth of an oligocheate on activated sludge. Water Res. 27(5):739–747. Richtmeyer, R., and K. Morton. 1967. Difference methods for initial-value problems. New York. Sharpe, F. R., and A. J. Lotka. 1911. Numerical bifurcation analysis of a tri-trophic food web with omnivory. Philos. Mag. 21:435–438. Sinko, J. W., and W. Streifer. 1971. A model for populations reproducing by fission. Ecology 52:330–335. Sulsky, D. 1993. Numerical solution of structured population models I Age structure. J. Math. Biol. 31:817–839. Sulsky, D. 1994. Numerical solution of structured population models II Mass structure. J. Math. Biol. 32:491–514. Thieme, H. 2003. Mathematics in population biology. Princeton, NJ: Princeton University Press. Thieme, H. R. 1988. Well-posedness of physiologically structured population models for Daphnia magna. J. Math. Biol. 26:299–317. Tuljapurkar, S. 1997. Stochastic matrix models. In Structured-population models in marine, terrestrial, and freshwater systems, ed. S. Tuljapurkar and H. Caswell, 59–87. New York: Chapman & Hall. Tyrcha, J. 2001. Age-dependent cell cycle model. J. Theor. Biol. 213(1):89–101. Tyson, J. J., and K. B. Hannsgen. 1985. Global asymptotic stability of the size distribution in probabilistic models of the cell cycle. J. Math. Biol. 22:61–68. Tyson, J. J., and K. B. Hannsgen. 1986. Cell growth and division: A deterministic=probabilistic model of the cell cycle. J. Math. Biol. 23:231–246. 1430 1435 1440 1445 1450 1455 1460 1465 1470 44 Q7 B. W. Kooi and F. D. L. Kelpin van den Bosch, F., A. M. de Roos, and G. Gabriel. 1988. Cannibalism as a life boat mechanics. J. Math. Biol. 26:619–633. von Foerster, H. 1959. Some remarks on changing populations. In The kinetics of cellular proliferation, ed. F. Stohlman. New York: Grune and Stratton. Webb, G. F. 1985. Theory of nonlinear age-dependent population dynamics, vol. 89 of Monographs and textbooks in pure and applied mathematics. New York: Marcel Dekker. 1475 1480