A Hierarchical Kinetic Theory of Birth, Death and Fission in Age-Structured Interacting Populations.

We develop mathematical models describing the evolution of stochastic age-structured populations. After reviewing existing approaches, we formulate a complete kinetic framework for age-structured interacting populations undergoing birth, death and fission processes in spatially dependent environments. We define the full probability density for the population-size age chart and find results under specific conditions. Connections with more classical models are also explicitly derived. In particular, we show that factorial moments for non-interacting processes are described by a natural generalization of the McKendrick-von Foerster equation, which describes mean-field deterministic behavior. Our approach utilizes mixed-type, multidimensional probability distributions similar to those employed in the study of gas kinetics and with terms that satisfy BBGKY-like equation hierarchies.

Introduction

Ageing is an important controlling factor in populations of organisms ranging in size from single cells to multicellular animals. Age-dependent population dynamics, where birth and death rates depend on an organism’s age, are important in quantitative models of demography [29], biofilm formation [3], stem cell differentiation [40, 44], and lymphocyte proliferation and death [52]. For example, cellular replication is controlled by a cycle [36, 38, 50], while higher organisms give birth depending on their maturation time. For applications involving small numbers of individuals, a stochastic description of the age-structured population is also desirable. A practical mathematical framework that captures age structure, intrinsic stochasticity, and interactions in a population would be useful for modeling many applications.

Standard frameworks for analyzing age-structured populations include Leslie matrix models [6, 31, 32], which discretizes ages into discrete bins, and the continuous-age McKendrick-von Foerster equation, first studied by McKendrick [28, 34] and subsequently von Foerster [47], Gurtin and MacCamy [17, 18], and others [24, 49]. These approaches describe deterministic dynamics; stochastic fluctuations in population size are not incorporated. On the other hand, intrinsic stochasticity and fluctuations in total population are naturally studied via the Kolmogorov master equation [7, 27]. However, the structure of the master equation implicitly assumes exponentially distributed event (birth and death) times, precluding it from being used to describe age-dependent rates or age structure within the population. Evolution of the generating function associated with the probability distribution for the entire population has also been studied [4, 8, 39, 41]. While this approach, the Bellman-Harris equation, allows for age-dependent event rates, an assumption of independence precludes population-dependent event rates. More recent methods [19, 22, 23, 26] have utilized Martingale approaches, which have been used mainly to investigate the asymptotics of age structure, coalescents, and estimation of Malthusian growth rates.

What is currently lacking is a complete mathematical framework that can resolve the age structure of a population at all time points, incorporate stochastic fluctuations, and be straightforwardly adapted to treat nonlinear interactions such as those arising in populations constrained by a carrying capacity [45, 46]. In a recent publication [16], we took a first step in this direction by formulating a full kinetic equation description that captures the stochastic evolution of the entire age-structured population and interactions between individuals. Here, we generalize the kinetic equation approach introduced in [16] along two main directions. First, we quantify the corrections to the mean-field equations by showing that the factorial moments of the stochastic fluctuations follow an elegant generalization of the McKendrick-von Foerster equation. Second, we show how the methods in [16] can be extended to incorporate fission processes, where single individuals instantaneously split into two identical zero-age offspring. These methods are highlighted with cell division and spatial models.

In the next section, we give a detailed overview of the different techniques currently employed in age-structured population modeling. In Sect. 3, we use previous results [16] to show how the moments of age-structured population size obey a generalized McKendrick-von Foerster equation. In Sect. 4, we expand the kinetic theory for branching processes involving fission. In Sect. 5, we demonstrate how our theory of fission can be applied to a microscopic model of cell growth. In Sect. 6, we demonstrate how to incorporate spatial effects. Conclusions complete the paper.

Age-Structured Population Modelling

Here we review, compare, and contrast existing techniques of population modeling: the McKendrick-von Foerster equation, the master equation, the Bellman-Harris equation, Leslie matrices, Martingale methods, and our recently introduced kinetic approach [16].

McKendrick-von Foerster Equation

It is instructive to first outline the basic structure of the classical McKendrick-von Foerster deterministic model as it provides a background for a more complete stochastic picture. First, one defines such that is the expected number of individuals with age within the interval . The total number of organisms at time t is thus . Suppose each individual has a rate of giving birth that is a function of its age a. For example, may be a function peaked around the time of M phase in a cell cycle or around the most fecund period of an organism. Similarly, is an organism’s rate of dying, which typically increases with its age a.

The McKendrick-von Foerster equation is most straightforwardly derived by considering the total number of individuals with age in [0, a]: . The number of births per unit time from all individuals into the population of individuals with age in [0, a] is , whilst the number of deaths per unit time within this cohort is . Within a small time window , the change in N(a, t) isIn the limit, we findUpon taking of Eq. 2, we obtain the McKendrick-von Foerster equation:The associated boundary condition arises from setting in Eq. 2:Finally, an initial condition completely specifies the mathematical model.

Note that the term on the right-hand side of Eq. 3 depends only on death; the birth rate arises in the boundary condition (Eq. 4) since births give rise to age-zero individuals. These equations can be formally solved using the method of characteristics. The solution to Eqs. 3 and 4 that satisfies a given initial condition isTo explicitly identify the solution, we need to calculate the fecundity function B(t). By substituting Eq. 5 into the boundary condition of Eq. 4 and defining the propagator , we obtain the following Volterra integral equation:After Laplace-transforming with respect to time, we findSolving the above for and inverse Laplace-transforming, we find the explicit expressionwhich provides the complete solution when used in Eq. 5.

The McKendrick-von Foerster equation is a deterministic model describing only the expected age distribution of the population. If one integrates Eq. 3 across all ages and uses the boundary conditions, the rate equation for the total population is . Generally, n(t) will diverge or vanish in time depending on the details of and . In the special case , the population is constant.

What is missing are interactions that stabilize the total population. Eqs. 3 and 4 assume no higher-order interactions (such as competition for resources, a carrying capacity, or mating patterns involving pairs of individuals) within the populations. Within the McKendrick-von Foerster theory, interactions are typically incorporated via population-dependent birth and death rates, and , respectively [11, 17, 18]. The McKendrick-von Foerster equation must then be self-consistently solved. However, as shown in [16], this assumption is an uncontrolled approximation and inconsistent with a detailed microscopic stochastic model of birth and death.

Master Equation Approach

A popular way to describe stochastic birth-death processes is through a function defining the probability that a population contains n identical individuals at time t. The evolution of this process can then be described by the standard forward continuous-time master equation [7, 27]where and are the birth and death rates, per individual, respectively. Each of these rates can be population-size- and time-dependent. As such, Eq. 9 explicitly includes the effects of interactions. For example, a carrying capacity can be implemented into the birth rate through the following form:Here, we have allowed both the intrinsic birth rate and the carrying capacity K(t) to be functions of time. Eq. 9 can be analytically or numerically solved via generating function approaches, especially for simple functions and .

Since only describes the total number of individuals at time t, it cannot resolve the distribution of ages within the fluctuating population. Another shortcoming is the implicit assumption of exponentially distributed waiting times between birth and death events. The times since birth of individuals are not tracked. General waiting time distributions can be incorporated into a master equation approach by assuming an appropriate number of internal “hidden” states, such as the different phases in a cell division cycle [50]. After all internal states have been sequentially visited, the system makes a change to the external population-size state. The waiting time between population-size changes is then a multiple convolution of the exponential waiting-time distributions for transitions along each set of internal states. The resultant convolution can then be used to approximate an arbitrary waiting-time distribution for the effective transitions between external states. It is not clear, however, how to use such an approach to resolve the age structure of the population.

Bellman-Harris Fission Process

The Bellman-Harris process [4, 8, 25, 39, 41] describes fission of a particle into any number of identical daughters, such as events , , and in Fig. 1a. Unlike the master equation approach, the Bellman-Harris branching process approach allows interfission times to be arbitrarily distributed. However, it does not model the budding mode of birth indicated by process in Fig. 1a, nor does it capture interactions (such as carrying capacity effects) within the population. In such a noninteracting limit, the Bellman-Harris fission process is most easily analyzed using the generating function F(z, t) associated with the probability , defined asWe assume an initial condition consisting of a single, newly born parent particle, . If we also assume the first fission or death event occurs at time , we can define as the generating function conditioned on the first fission or death occurring at time and write F recursively [1, 2, 20] as:The function H encapsulates the probability that a particle splits into m identical particles upon fission, for each non-negative integer m. For binary fission, we have since . Since this overall process is semi-Markov [48], each daughter behaves as a new parent that issues its own progeny in a manner statistically equivalent to and independent from the original parent, giving rise to the compositional form in Eq. 12. We now weight over a general distribution of waiting times between splitting events, , to findThe Bellman-Harris branching process [2, 14] is thus defined by two parameter functions: , the vector of progeny number probabilities, and , the probability density function for waiting times between branching events for each particle. The probabilities can be recovered using a contour integral (or Taylor expanding) about the origin:Note that Eq. 13 incorporates an arbitrary waiting-time distribution between events, a feature that is difficult to implement in the master equation (Eq. 9). An advantage of the branching process approach is the ease with which general waiting-time distributions, multiple species, and immigration can be incorporated. However, it is limited in that an independent particle assumption was used to derive Eq. 13, where the statistical properties of the entire process starting from one parent were assumed to be equivalent to those started by each of the identical daughters born at time . This assumption of independence precludes treatment of interactions within the population, such as those giving rise to carrying capacity. More importantly, the Bellman-Harris equation is expressed purely in terms of the generating function for the total population size and cannot resolve age structure within the population.

Fig. 1

a A general branching process. I indicates a budding or simple birth process, where the parental individual produces a single offspring (a ‘singlet’) without death. II indicates binary fission, where a parent dies at the same moment two newborn twins occur (a ‘doublet’). III indicates a more general fission event with four offspring (a ‘quadruplet’). IV indicates death, which can be viewed as fission with zero offspring. b A binary fission process such as cell division. At time we have four individuals; two sets of twins. At time we have six individuals; two pairs of twins and two singlets

Leslie Matrices

Leslie matrices [31, 32] have been used to resolve the age structure in population models [9, 10, 12, 15, 31–33, 39, 44]. These methods essentially divide age into discrete bins and are implemented by assuming fixed birth and death rates within each age bin. Such approaches have been applied to models of stochastic harvesting [10, 15] and fluctuating environments [13, 30] and are based on the following linear construction, iterated over a single time step:The value indicates the population size in age group i; is the mean number of offspring arriving to age group 0 from a parent in age group i; and is the fraction of individuals surviving from age group i to . These models have the advantage of being based upon algebraic linearity, which enables many features of interest to be investigated analytically [6]. However, they are inherently deterministic (although they can be used to study extrinsic environmental noise) and the discretization within such models results in an approximation. Thus, a fully continuous stochastic model is desirable.

Martingale Approaches

Relatively recent investigations have used Martingale approaches to model age-structured stochastic processes. These methods stem from stochastic differential equations and Dynkin’s formula [37] and evolves general processes of the form . Here the vector represents the time dependent age-chart of the population with variable size n; f is a symmetric function of the individual ages; and F is a generic function of interest. A Martingale decomposition of the following form resultswhere the operator captures the mean behavior, and the stochastic behavior is encoded in the local Martingale process [26]. Such analyses have enabled several features of general birth-death processes, including both budding and fission forms of birth to be quantified. Specifically, the Malthusian growth parameter can be explicitly determined, along with the asymptotic behavior of the age structure. More recently there have been results related to coalescents and extinction of these processes [19, 22, 23]. However, we will show the utility of obtaining the probability density of the entire age chart of the population which allows efficient computations in transient regimes. The kinetic approach first developed in [16] introduces machinery to accomplish this.

Kinetic Theory

A brief introduction to the current formulation of our kinetic theory approach to age-structured populations can be found in [16]. The starting point is a derivation of a variable-dimension coupled set of partial differential equations for the complete probability density function describing a stochastic, interacting, age-structured population subject to simple birth and death. Variables in the theory include the population size n, time t, and the vector representing the complete age-chart for the n individuals. If we randomly label the individuals , then represents the probability that the individual has age in the interval . Since individuals are considered indistinguishable, is invariant under any permutation of the age-chart vector . These functions are analogous to those used in kinetic theories of gases [35]. Their analysis in the context of age-structured populations builds on the Boltzmann kinetic theory of Zanette [51] and results in the kinetic equationwhere and the age variables are separated from the time variable by a semicolon. The associated boundary condition is

Note that because is symmetric in the age arguments, the zero can be placed equivalently in any of the n age coordinates. The birth rate function can be quite general and can take forms such as for a simple birth process or to represent births arising from interactions between pairs of individuals.

Equation 17 applies only to the budding or simple mode of birth such as event in Fig. 1a. In [16] we derived analytic solutions for in pure death and pure birth processes, and showed that marginal densities obeyed a BBGKY-like (Bogoliubov-Born-Green-Kirkwood-Yvon) hierarchy of equations. Furthermore, when the birth and death rates are age-independent (but possibly number-dependent), the hierarchy of equations reduce to a single master equation for the total number of individuals n in the population. Characterizing all the remaining higher moments of the distribution remains an outstanding problem. Moreover, methods to tackle fission modes of birth such as those shown in Fig. 1b were not developed. These are the two contributions described in this paper. Before analyzing these problems, we summarize the pros and cons of the different approaches in Table 1.

Table 1

Advantages and disadvantages of different frameworks for stochastic age-structured populations. ‘Stochastic’ indicates that the model resolves probabilities of configurations of the population

Theory	Stochastic	Age-dependentrates	Age-structuredpopulations	Age-chart resolved	Interactions	Budding	Fission
Verhulst Eq.	×	×	×	×	✓	×	×
McKendrick Eq.	×	✓	✓	×	✓	✓ 1	×
Master Eq.	✓	×	×	×	✓	✓	✓
Bellman-Harris	✓	✓	×	×	×	×	✓
Leslie matrices	×	✓ 2	✓	×	✓	×	×
Martingale	✓	✓	× 3	×	✓	✓	✓
Kinetic theory	✓	✓	✓	✓	✓	✓	✓ 4

‘Age-dependent rates’ indicates whether or not a model takes into account birth, death, or fission rates that depend on an individuals age (time after its birth). ‘Age-structured Populations’ indicates whether or not the theory outputs the age structure of the ensemble population. ‘Age-chart resolved’ indicates whether or not a theory outputs the age distribution of all the individuals in the population. ‘Interactions’ indicates whether or not the approach can incorporate population-dependent dynamics such as that arising from a carrying capacity, or from birth processes involving multiple parents. ‘Budding’ and ‘Fission’ describes the model of birth and indicates whether the parent lives or dies after birth

1 Birth and death rates in the McKendrick-von Foerster equation can be made explicit functions of the total populations size, which must be self-consistently solved [17, 18]

2 Leslie matrices discretize age groups and are an approximate method

3 Martingale methods do not resolve the age structure explicitly, but utilize rigorous machinery

4 The kinetic approach for fission is addressed later in this work, but not in [16]

Analysis of Simple Birth-Death Processes

We now revisit the simple process of budding birth and death, and extend the kinetic framework introduced in [16]. We first show that the factorial moments for the density satisfy a generalized McKendrick-von Foerster equation. We also explicitly solve Eqs. 17 and 18, and derive for the first time an exact general solution for .

Moment Equations

The McKendrick-von Foerster equation has been shown to correspond to a mean-field theory of age-structured populations in which the birth and death rates and are population-independent [16]. This leaves open the problem of determining the age-structured variance (and higher-order moments) of the population size.

In [16], we derived the marginal k-dimensional distribution functions defined by integrating over age variables:The symmetry properties of indicate that it is immaterial which of the age variables are integrated out. From Eq. 17, we then obtainedSimilarly, integrating the boundary condition in Eq. 18 over of the (nonzero) variables, gives, for simple birth processes where ,We now show how to use these marginal density equation hierarchies and boundary conditions to derive an equation for the moment of the age density.

For , is the probability that we have n individuals and that if one of them is randomly chosen, it will have age in . Therefore, the probability that we have n individuals, and that there exists an individual with age in , is . Summing over all possible population sizes yields the probability that the system contains an individual with age in the interval . More generally, is the probability that there are n individuals, k of which can be labelled such that the has age within the interval . Summing over the possibilities , we thus introduce factorial moments and moment functions as:Here is the Pochhammer symbol, and and are Stirling numbers of the first and second kind, respectively [42, 43]. Although we are primarily interested in the functions , the factorial moments will prove to be analytically more tractable. One can then easily interchange between the two moment types by using the polynomial relationships involving Stirling numbers.

After multiplying Eq. 20 by and summing over all , we findwhere, for simplicity of notation, the arguments have been suppressed from and . In the case where the birth and death rates and are independent of the sample size, significant cancellation occurs and we find the simple equationWhen , one recovers the classical McKendrick-von Foerster equation describing the mean-field behavior after stochastic fluctuations are averaged out. Equation 24 is a natural generalization of the McKendrick-von Foerster equation and provides all the age-structured moments arising from the population size fluctuations. If the birth and death rates, and , depend on the population size, one has to analyze the complicated hierarchy given in Eq. 23.

To find the boundary conditions associated with Eq. 24, we combine the definition of with the boundary condition in Eq. 21 and obtainSetting , we recover the boundary condition associated with the classical McKendrick-von Foerster equation. For higher-order factorial moments, the full solution to the factorial moment is required for the boundary condition to the moment .

Specifically, consider the second factorial moments and assume the solution to the McKendrick-von Foerster equation is available (from e.g., Eq. 5). In the infinitesimal interval , the term is the Bernoulli variable for an individual having an age in the interval . Thus, in an extended age window , we heuristically obtain the expectationwhere is the stochastic random variable describing the number of individuals with an age in at time t. Using an analogous argument for the variance, we findThus, the second moment allows us to describe the variation of the population size within any age region of interest. Similar results apply for higher order correlations. We focus then on deriving a solution to and determining the variance of population-size-age-structured random variables. Eq. 24 for general k is readily solved using the method of characteristics leading towhere the propagator is defined as . We can now specify in terms of boundary conditions and initial conditions by selecting . Since is invariant to any permutation of its age arguments, we have only two conditions to consider. The initial condition encodes the initial correlations between the ages of the founder individuals and is assumed to be given. From Eq. 22, must be a symmetric function in the age arguments. A boundary condition of the form describes the fecundity of the population through time. This is not given but can be determined in much the same way that Eq. 8 was derived.

To be specific, consider a simple pure birth (Yule-Furry) process (, ) started by a single individual. The probability distribution of the initial age of the parent individual is assumed to be exponentially distributed with mean . Upon using transform methods similar to those used to derive Eq. 8, we obtain the following factorial moments (see Appendix 1 for more details):We have given for only since the region can be found by imposing symmetry of the age arguments in . After using Eq. 22 to convert and into and , we can use Eqs. 26 and 27 to find age-structured moments, particularly the mean and variance for the number of individuals that have age in the interval [a, b]:Note that in the limits and , we recover the expected exponential growth of the total population size for a Yule-Furry process. We also recover the known total population variance .

Full Solution

Equation 17 defines a set of coupled linear integro-differential equations in terms of the density . In [16], we derived explicit analytic expressions for in the limits of pure death and pure birth. Here, we outline the derivation of a formal expression for the full solution. To do so, it will prove useful to revert to the following representation for the density:If is restricted to the ordered region such that , can be interpreted as the probability density for age-ordered individuals (see [16] for more details). We will consider as a distribution over ; however, its total integral (n!) is not unity and it is not a probability density. We can use Eq. 32 to switch between the two representations, but simpler analytic expressions for solutions to Eq. 17 result when is used.

To find general solutions for expressed in terms of an initial distribution, we replace with in Eq. 17 and use the method of characteristics to find a solution. Examples of characteristics are the diagonal timelines portrayed in Fig. 2. So far, everything has been expressed in terms of the natural parameters of the system; the age of the individuals at time t. However, varies in time complicating the analytic expressions. If we index each characteristic by the time of birth (TOB) instead of age a, then b is fixed for any point (a, t) lying on a characteristic, resulting in further analytic simplicity. We use the following identity to interchange between TOB and age representations:We will abuse notation throughout our derivation by identifying . The method of characteristics then solves Eq. 17 to give a solution of the following form, for any This equation is defined in terms of a propagator that represents the survival probability over the time interval , for m individuals born at times , in a population of size n,where we have again used the definition . The propagator satisfies certain translational properties:The solution applies to any region of phase space where . If , say , then we must invoke the boundary conditions of Eq. 18 to replace with , where we have and will henceforth use the notationEq. 34 is then used to propagate backwards in time. To obtain a general solution, we need to repeatedly back-substitute Eq. 34 and the associated boundary condition, resulting in an infinite series of integrals. However, each term in the resultant sum can be represented by a realization of the birth-death process. We represent any such realization across time period [0, t], such as that given in Fig. 2, as follows.

Fig. 2

A sample birth death process over the time interval [0, t]. Red and white circles indicate births and deaths within [0, t]. The variables and denote TOBs of individuals present at time t, while , and indicate birth and death times of individuals who have died by time t. Terms arising from application of the recursion in Eq. 34 and boundary condition of Eq. 18 are given to the right

Let and denote the TOBs for m individuals born in the time interval [0, t], and n founder individuals, all alive at time t. Next, define and to be the TOBs of k individuals born in the time interval [0, t] and founder individuals, respectively. Here, all individuals are assumed to die in the time window [0, t]. Their corresponding times of death are defined as and , respectively. Thus, there will be individuals alive initially at time and individuals alive at the end of the interval [0, t].

Next, consider the realization in Fig. 2, where we start with the two individuals at time 0 with TOBs and . The individual with TOB survives until time t, while the individual with TOB dies at time . Within the time interval [0, t] there are three more births with TOBs , and , the last of which has a corresponding death time of , resulting in three individuals in total that exist at time t.

To express the distribution in terms of the initial distribution , conditional upon three birth and two death events ordered such that , we start with the distribution . Just prior to the first birth time , we have two individuals, so that and Eq. 34 yields (the death term does not contribute). To describe the birth at time , we use the boundary condition of Eq. 18 to construct .

Immediately after and before the next death occurs at time , three individuals exist and . Now, only the death term in Eq. 34 contributes andContinuing this counting, we find the product of terms displayed on the right-hand side of Fig. 2.

Next, we use the translational properties indicated in Eqs. 36 and 37 to combine the propagators associated with Fig. 2 into one term: . In other words, each birth-death pair (y, s) is propagated along the time interval it survives; from to . For example, the individual with TOB survives across the entire timespan [0, t], whereas the individual with TOB is born and dies at times and . These two individuals are propagated by the terms and , respectively. Provided the order is preserved and the values are negative, the form of the integral expressions in Fig. 2 are preserved.

After summing across all realizations (the configuration in Fig. 2 is one member of ) of the possible orderings of the birth and death times , , , and , we can write the general solution to Eq. 34 in the formThe terms and refer to the times below and above relative to the ordering of times , , , and . For example, in Fig. 2, and represent the lower and upper bounds of the vector found from the ordering . The term A(x) represents the vector of TOBs of the individuals alive just prior to time x. The term N(x) represents the number of individuals alive just prior to time x.

Although analytic and complete, the solution given in Eq. 40 is unwieldy and difficult to implement. One can truncate the sum to remove low probability contributions, such as realizations containing improbable numbers of intermediary births and deaths, and perform numerical integration. However, this approach also rapidly becomes infeasible as the dimensions increase. Therefore, we explore the general solution via event-based Monte-Carlo simulation. We initialize the process with a number of samples obtained from an initial distribution. Each sample is represented by a vector of birth times and is propagated forward in time. A timestep is chosen to be sufficiently small such that at most one birth or death event occurs within it, after which the vector is updated. This process is continued until the required time has been reached. Although the high dimensionality makes it difficult to sample enough realizations to sufficiently explore the distribution , lower dimensional marginal distributions such as , and , and their counterparts , can be sufficiently sampled.

Figure 3a, b show results from simulations of a pure death and a pure birth process, respectively. In Fig. 3a we assumed a population-independent linear death rate and initiated the pure death process with 10 individuals with initial ages drawn from a gamma distribution with unit mean and standard deviation . Fig. 3a(i) shows the simulated density which decreases in n with time. Figs. 3a(ii–iv) show that the weight of the reduced single-particle density function shifts to longer times and higher ages as the system size n is decreased. The sum over the population of the symmetric two-point correlation is shown in Fig. 3a(v). The observed structure indicates no correlations in the death only process and the peak at reflects the fact that older individuals die faster, shifting the mean age slightly below the initial age plus the elapsed time (1 + 2 = 3). Fig. 3b shows results from Monte-Carlo simulations of a pure birth process with growth rate and carrying capacity (Eq. 10). Here, we initiated the simulations with one individual with age drawn from the same gamma distribution . In this case, the reduced single-particle density exhibits peaks arising from both from the initial distribution and from birth (Fig. 3b(ii–iv)). The two-point correlation function exhibits a similar multimodal structure as shown in (v). In all simulations at least 400,000 trajectories were aggregated and the results are in good agreement with analytic solutions to Eq. 17.

Fig. 3

Monte–Carlo simulations of densities in age- and number-dependent birth-death processes. Column a shows results for a death-only process with a linear death rate function . We initiated all simulations from 10 individuals with initial age drawn from distribution . In column b, we consider a budding-only birth process with a carrying capacity (in Eq. 10). Here, simulations were initiated with a single parent individual with an initial age also drawn from the distribution P(a). In (i), we plot the total number density for both processes. We also plot the single-particle density function for the pure death process in a(ii–iv) and for the limited budding process in b(ii–iv). Finally, the population-summed two-point correlations functions for pure death and pure budding are shown in panels a(v) and b(v)

Age-Structured Fission-Death Processes

We now derive a kinetic theory for a binary fission-death process, as depicted in Fig. 1b. We find a hierarchy of kinetic equations, analogous to Eqs. 17 and 18, and determine the mean behavior.

Extended Liouville Equation for Fission-Death

The binary fission-death process is equivalent to a birth-death process except that parents are instantaneously replaced by two newborns. The process can also be thought of as a budding process in which the parent is instantaneously renewed. In order to describe both twinless individuals (singlets) and twins (a doublet), we have to double the dimensionality of our density functions. For example, in Fig. 1b at time , we have two pairs of distinct twins, with four individuals having two ages, whereas at time we have two singlets and two doublets. Thus, we define the ages of current singlets and twins by and , respectively, where m is the number of singlets and n the number of pairs of twins. Transforming to the time-of-birth (TOB) representation, we define the TOB of current singlets and twins as and , respectively. For simplicity, we will assume that no simple birth processes occur and that particles grow in number only through fission. The function is defined as the age-dependent fission rate of an individual (whether a singlet or a doublet) of age a when the system contains m singlets and n doublets. Similarly, we have death rate , and event rate . We suppose, for the moment, that the TOBs are ordered so that and . The quantity is then the probability of m singlets with ordered TOBs in and n twin pairs with ordered TOBs in . The density satisfies the following equation:where the partial age vectors are defined as and the singlet age vector, doublet age vector, and time arguments are separated by semicolons. The term represents the vector of all m singlet TOBs, except for the ith one. The first term on the right hand side of Eq. 41 represents the death of a singlet particle with an unknown TOB z in the interval , while the second term describes the death of any one of two individuals in a pair of twins (with TOB ).

The associated boundary conditions areThe first term on the right-hand side above represents the fission of one of a pair of twins, generating a new pair of twins of age zero (TOB t), and leaving behind a singlet with TOB . The second term represents the fission (and removal) of a singlet with unknown TOB z, giving rise to an additional pair of twins of age zero.

We now let and be unordered TOB vectors, and extend to the domain by defining , where is the ordering operator. Note that is not a probability distribution under this extension; however, can be interpreted as the probability that we have a population of m singlets and n pairs of twins, such that if we randomly label the singlets and the doublets , the singlet has age in and the doublet have age in . The density obeyswith associated boundary conditionEquations 44 and 45 provide a complete probabilistic description of the population of singlets and doublets undergoing fission and death.

Mean-Field Behavior

Here, we analyze the mean-field behavior of the fission-death process by first integrating out unwanted variables from the full density to construct marginal or “reduced” densities. Successive integrals over any number of the variables and can be performed, giving:For example, is the probability of finding at time t, m singlets and n doublets, regardless of age. After integrating Eq. 44 we find the double hierarchy of equationsSimilarly, integrating Eq. 45 yields boundary conditions for the marginal densities:We can now analyze the densities X(x, t) and Y(y, t), where is the probability that there exists at time t a singlet with TOB in and is the probability that at time t we have one doublet with TOB in . Analogous to Eq. 22, we defineUpon setting and , we multiply Eq. 47 by m and n, respectively, and sum both equations. If the fission and death rates and depend on population size, the resultant expressions are complex hierarchies which will be difficult to analyze. However, if and are size-independent, many cancellations occur and the resulting equations for X and Y simplify significantly, givingSimilarly, repeating the operation on the boundary conditions in Eq. 48, we find boundary conditions for X and Y:Note that if is the total population density, Eqs. 50 and 51 reduce to McKendrick-von Foerster-like equations:To solve Eqs. 50 and 51, we first defineand find solutions of the formprovided . For an initial time of , we find, upon setting ,We now substitute Eqs. 55 and 56 into Eqs. 51 to find a Volterra equation for B(t):Equation 57 along with Eqs. 55 and 56 constitute a complete solution for the mean density of singlets and doublets. Eqs. 55 and 56 also show that the total population density, , takes on a simple form in terms of B(t):while the total mean population is given byBefore analyzing a specific model of the fission-death process, we will first establish the equivalence of our noninteracting kinetic theory with the Bellman-Harris fission process (discussed in Sect. 2.3) in the mean-field limit.

Mean-Field Equivalence to the Bellman-Harris Process

Consider a Bellman-Harris fission process with an inter-branching time distributed according to the function and an associated cumulative density function defined by . Upon using the progeny distribution function given in Eq. 12, the Bellman-Harris model in Eq. 13 can be written equivalently asIf we restrict ourselves to a binary fission process, the progeny distribution function takes the form , where and are the death and binary fission probabilities, conditional on an event taking place. Thus, the mean population defined ashas the Laplace-transformed solutionWe now show that the same result arises from our full noninteracting (population-independent and ) kinetic approach. Since the fission and death rates can be expressed as and , Eq. 53 reduces to and . Starting from a single individual with age zero, Eq. 59 can be written aswhich has the Laplace-transformed solutionSimilarly, Eq. 57 becomeswith Laplace-transformed solutionSubstituting Eq. 66 in Eq. 64 results in Eq. 62 for , explicitly establishing the mean-field equivalence between the Bellman-Harris approach and our kinetic theory. Note that in the Bellman-Harris formulation, the waiting-time distributions of either fission or death have the same distribution g(a). In our kinetic theory, these rates can have distinct distributions, and , and can also depend on population size, providing much greater flexibility.

A Fission-Only Model of Cell Division

We now consider explicit results for a simple fission-only model () of cell division where cell cycle times are rescaled to be -distributed with unit mean and variance . This -distribution and its Laplace transform are explicitlyEquation 66 for B(t) can then be solved to yieldThe inverse Laplace transform is detailed in Appendix 2 and involves contour integration that yieldsSimilarly, from Eq. 62 we havewhich can also be evaluated via a similar Bromwich integral:For , is exponentially distributed, and we find the simple growth law , which is equivalent to the result found earlier in Sect. 3.1. This corresponds to a continuously compounded population. On the other hand, when is increased, the -distribution sharpens about unity. Figs. 4a–c show that as increases, the mean population size T(t) tends towards that given by the discrete-time Galton–Watson step process, as would be expected. In the limit, the population compounds at discrete, evenly timed intervals leading to an overall lower population compared to that of a process with more frequent branching (smaller ). In Figs. 4d–f, we have used the expression for B(t) in Eqs. 58 and 69 to give the mean age-time distribution T(a, t). Note that unlike the solution to the Bellman-Harris equation shown in Figs. 4a–c, the mean density T(a, t) (Eq. 58) resolves age structure.

Fig. 4

Plots of simulations and analytic results of a fission-only process with -distributed branching times. a, b, and c show mean populations as a function of time for dispersion values , , and , respectively. Red dotted trajectories are realizations of simulations, while the solid red line is the mean. The blue dashed curve is the mean population T(t) computed from Eq. 71 and is nearly indistinguishable from the red solid curve. The upper and lower black lines correspond to the continuous-time Markovian fission process and the discrete-time Galton-Watson process, respectively. d, e, and f depict the corresponding mean age-distributions T(x, t) computed from Eq. 58 but plotted as functions of time t and age a

Spatial Models

We now illustrate how our age-structured kinetic model can be generalized to include spatial motion such as diffusion and convection. We will follow the approaches described in Webb [49] for incorporating spatial effects in age-structured simple birth-death processes. Since these methods are adaptations of the McKendrick-von Foerster equation, they are deterministic and ignore stochastic fluctuations in population size. In a manner similar to how the McKendrick-von Foerster equation was extended to the stochastic domain using Eq. 17, here, we outline how to generalize the age-structured spatial process discussed in [49] to incorporate stochasticity.

Consider a simple budding-mode birth-death process such that is the probability density for a population containing n randomly labelled individuals with TOBs and positions . Although is again invariant under permutations of variables associated with different individuals, the relative orders of and must be preserved. For example, . For ease of presentation, we assume a one-dimensional system; generalizations to higher spatial dimensions are straightforward. We further suppose that individuals are undergoing identical, independent diffusion processes with diffusion constant D. Examples of other spatial processes that may be combined with stochastic age-structured kinetics can be found in [49]. We suppose that and are birth and death rates for an individual with age a and at spatial position q in a population of size n. Finally, the initial position of each newborn is determined by the position of the parent at the time of birth. The extended theory is described by the following kinetic equation for :The corresponding boundary condition capturing the influx of newborn individuals iswhich differs slightly from that in Eq. 18. In the original formulation, we do not track which individual is the parent of a newborn, whereas here the newborn has the same position () as the parent (), setting its identity as the individual. In addition to a boundary condition, Eq. 72 requires an initial condition to specify both the initial TOB and initial position of individuals.

As with our earlier analyses, we first express in terms of by introducing the propagator , which enables us to transform Eq. 72 to an inhomogeneous heat equation for the function ,whose solution can be expressed in the form [5]Here, denotes the identity matrix and is the multivariate normal density for the vector arising from a distribution with mean and covariance . This result expresses in terms of and is analogous to Eq. 34. This solution is valid provided ; for , we must invoke the boundary condition. One can then use Eq. 75 and the boundary condition to search for explicit solutions in much the same way as we did for our spatially independent kinetic theory.

Summary and Conclusions

We have developed a complete kinetic theory for age-structured birth-death and fission-death processes that allow for systematic and and self-consistent incorporation of interactions at the population level. Our overall result in [16], which we extend here, is the derivation of a kinetic theory for stochastic age-structured populations. The kinetic equations can be written in terms of a BBGKY-like hierarchy (or a double hierarchy in the case of fission). Methods of approximation and closure typically employed in gas/liquid kinetic theory, plasma physics, or fluid dynamics can then be applied.

The analysis presented in this paper provides three new results. First, in Eq. 24, we have shown that the factorial moments of the age structure can be described by an equation that naturally generalizes the McKendrick-von Foerster equation. In particular, for population-independent birth, death, and fission rates we can determine the variance of the population size for specific age groups in a population, something that was not previously feasible without some form of approximation.

Second, in Eqs. 17 and 18, we develop a complete probabilistic description of a population undergoing a binary fission and death process. Although a general analytic solution to these systems can be written down (Eq. 40), it is difficult to calculate and further work is needed to identify analytic techniques or numerical schemes that can more readily provide solutions. The methods we have introduced can also be viewed as a continuum limit of matrix population models.

Third, we also outlined how to incorporate spatial dependence of birth and death into our age-structured kinetic theory. We considered only the simplest model of free diffusion in which individuals to not interact spatially. Spatially-mediated interactions can be incorporated by way of a “collision operator” in a full theory that treats both age and space kinetically.

Finally, we note that the overall structure of our model is semi-Markov. That is, birth, death, and fission rates depend on only the time since birth of an individual and not on, for example, the number of generations removed from a founder. Such lineage aging processes are often important in cell biology e.g., the Hayflick limit [21]) and would require extension of our state space to include generational class [52]. These extensions will be explored in future work.