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 with individuals that range in size from single cells to multicellular organisms.Age-dependent population dynamics, where birth and death rates depend on an organism's age, are important in quantitative models of demography [33], biofilm formation [3], stem cell differentiation [45,49], and lymphocyte proliferation and death [56].In cell-based applications, replication is governed by a time-dependent cell cycle [40,43,54], while for higher organisms, the ability to give birth depends on their maturation time.For applications involving small numbers of individuals, a stochastic description of the age-structured population is also desirable.Finding 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,35,36], which discretizes ages into discrete bins and the continuous-age McKendrick-von Foerster equation and first studied by McKendrick [32,38] and subsequently von Foerster [16], Gurtin and MacCamy [21,22], and others [28,53].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,31].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 have also been developed [4,8,44,46].Although this approach, the Bellman-Harris equation, allows for age-dependent event rates, an assumption of independence precludes population-dependent event rates.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 t1 we have four individuals; two sets of twins.At time t2 we have six individuals; two pairs of twins and two singlets.
More recent methods [23,26,27,30] have utilized Martingale approaches, which have been used mainly to investigate the asymptotics of age structure, coalescents, and estimation of Malthusian growth rate parameters.
Thus, it is desirable to have 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 [50,51].In a recent publication [20], 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 [20] 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 McKendrickvon Foerster equation.Second, we show how the methods in [20] 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.We also draw attention to the companion paper [19], where quantum field theory techniques developed by Doi and Peliti [13,14,42] are used to address the same problem, providing alternative machinery for age-structured modeling.
In the next section, we give a detailed overview of the different techniques currently employed in age-structured population modeling.In Section 3, we use previous results [20] to show how the moments of age-structured population size obey a generalized McKendrick-von Foerster equation.In Section 4, we develop the kinetic theory for branching processes involving fission.In Section 5, we demonstrate how our theory can be applied to a microscopic model of cell growth.In Section 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 McKendrickvon Foerster equation, the master equation, the Bellman-Harris equation, Leslie matrices, Martingale methods, and our recently introduced kinetic approach [20].

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 ρ(a, t) such that ρ(a, t)da is the expected number of individuals with an age within the interval [a, a + da].The total number of organisms at time t is thus n(t) = ∞ 0 ρ(a, t)da.Suppose each individual has a rate of giving birth β(a) that is a function of its age a.For example, β(a) 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, its rate of dying µ(a) can also depend on its age, and will typically increase with age.
The McKendrick-von Foerster equation is most straightforwardly derived by considering the total number of individuals with age in [0, a]: N (a, t) = a 0 ρ(y, t)dy.The number of births per unit time from all individuals into the population of individuals with age in [0, a] is B(t) = ∞ 0 β(y)ρ(y, t)dy, whilst the number of deaths per unit time within this cohort is D(a, t) = a 0 µ(y)ρ(y, t)dy.Within a small time window ε, the change in N (a, t) is In the ε → 0 limit, we find Upon taking ∂ ∂a of Eq.
The associated boundary condition arises from setting a = 0 in Eq. 2: Finally, an initial condition ρ(a, t = 0) = g(a) 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 is To explicitly identify the solution, we need to calculate the fecundity function B(t).By substituting the solution in Eq. 5 into the boundary condition of Eq. 4 and defining the propagator U (a 1 , a 2 ) ≡ exp − a2 a1 µ(s)ds , we obtain the following Volterra integral equation: After Laplace-transforming with respect to time, we find Solving the above for B(s) and inverse Laplace-transforming, we find the explicit expression which 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 0 ≤ a < ∞ and uses the boundary conditions, the rate equation for the total population is ṅ(t) = ∞ 0 (β(a) − µ(a))ρ(a, t)da.Thus, n(t) will in general eventually diverge or vanish depending on the details of β(a) and µ(a).In the special case β(a) = µ(a), the population is constant.
What is missing are interactions that stabilize the total population.Eqs. 3 and 4 assume no higherorder 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 heuristically incorporated by assuming that the birth and death rates (β(a; n(t)) and µ(a; n(t))) can depend on the total population n(t) [11,21,22].However, as shown in [20], 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 ρ n (t) 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,31] where β n (t) and µ n (t) 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 β 0 (t) 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 β n and µ n .Since ρ n (t) 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 [54].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,29,44,46] describes fission of a particle into any number of identical daughters, such as events II, III, and IV 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 I 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 ρ n (t), defined as We assume an initial condition consisting of a single, newly born parent particle, ρ n (0) = δ n,1 .If we also assume the first fission or death event occurs at time τ , we can define F (z, t|τ ) as the generating function conditioned on the first fission or death occurring at time τ and write F recursively [1,2,24] as: The function A encapsulates the probability a m that a particle splits into m identical particles upon fission, for each non-negative integer m.For binary fission, we have Since this overall process is semi-Markov [52], 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 F (z, t|τ ) over a general distribution of waiting times between splitting events, g(τ ), to find The Bellman-Harris branching process [2,17] is thus defined by two parameter functions: a m , the vector of progeny number probabilities, and g(τ ), the probability density function for waiting times between branching events for each particle.The probabilities ρ n (t) can be recovered using a contour integral surrounding or a Taylor expansion about the origin: Note that Eq. 13 allows one to incorporate 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.

Leslie Matrices
Leslie matrices [35,36] have been used to resolve the age structure in population models [9, 10, 12, 18, 35-37, 44, 49].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,18] and fluctuating environments [15,34] and are based on the following linear construction, iterated over a single time step: The value n i indicates the population size in age group i; f i is the mean number of offspring arriving to age group 0 from a parent in age group i; and s i is the fraction of individuals surviving from age group i to i + 1.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 of 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 [41] and considers general processes of the form F (f (a n (t))), where the vector a n (t) 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 results, where the operator G captures the mean behavior, and the stochastic behavior is encoded in the local Martingale process M (f,F ) t [30].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 [23,26,27].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 [20] 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 [20].The starting point is a derivation of a variable-dimension Liouville equation for the complete probability density function ρ n (a n ; t) 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 a n = (a 1 , a 2 , . . ., a n ) representing the complete age chart for the n individuals.If we randomly label the individuals 1, 2, . . ., n, then ρ n (a n ; t)da n represents the probability that the i th individual has age in the interval [a i , a i + da i ].Since individuals are considered indistinguishable, ρ n (a n ; t) is invariant under any permutation of the age-chart vector a n .These functions are analogous to those used in kinetic theories of gases [39].Their analysis in the context of age-structured populations builds on the Boltzmann kinetic theory of Zanette [55] and results in a BBGKY-like (Bogoliubov-Born-Green-Kirkwood-Yvon) hierarchy of equations: where γ n (a) = β n (a)+µ n (a) and the age variables are separated from the time variable by the semicolon.
The associated boundary condition is given by Note that because ρ n (a n−1 , 0; t) is symmetric in the age arguments, the zero can be placed equivalently in any of the n age coordinates.The birth rate function is quite general and can take forms such as for a simple birth process or 1≤i<j≤n−1 β n−1 (a i , a j ) 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 I in Fig. 1A.In [20] we derived analytic solutions for ρ n (a n ; t) in pure death and pure birth processes and showed that when the birth and death rates are constant, the hierarchy of equations reduces to a single master equation.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 main areas of focus of the present work.Before analyzing these problems, we summarize the pros and cons of the different techniques in Table 1.

Analysis of Simple Birth-Death Processes
We now revisit the simple budding birth and death process and extend the the kinetic framework introduced in [20].We first show that the factorial moments for the density ρ n (a n ; t) 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 ρ n (a n ; t).
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 '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 [21,22]. 2 Leslie matrices discretize age groups and are an approximate method. 3Martingale methods do not resolve the age structure explicitly, but utilize rigorous machinery. 4The kinetic approach for fission is addressed later in this work, but not in [20].Kinetic Theory 4

Moment Equations
The McKendrick-von Foerster equation has been shown to correspond to a mean-field theory of agestructured populations [20].Specifically, if we construct the expected population density X(a, t) n (a; t), we obtain the McKendrick-von Foerster equation for X(a, t).This leaves open the problem of determining the age-structured variance (and higher-order moments) of the population size.
In [20], we derived the marginal k−dimensional distribution functions defined by integrating ρ n (a n ; t) over n − k age variables: The symmetry properties of ρ n (a n ; t) indicate that it is immaterial which of the n − k age variables are integrated out.From Eq. 17, we then obtained n+1 (a k , y; t)dy.
Similarly, integrating the boundary condition in Eq. 18 over n − k 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 k th moment of the age density.For k = 1, ρ n (a; t)da is the probability that we have n individuals and that if one of them is randomly chosen, it will have age in [a, a + da].Therefore, the probability that we have n individuals, and there exists an individual with age in [a, a + da], is nρ n (a; t)da.Summing over all possible population sizes n ≥ 1 yields the probability ρ(a, t)da = n nρ (1) n (a; t)da that the system contains an individual with age in the interval [a, a + da].More generally, n k ρ (k) n (a k ; t)da k is the probability that there are n individuals, k of which can be labelled such that the i th has age within the interval [a i , a i + da i ].Summing over the possibilities n ≥ k, we thus introduce factorial moments X (k) (a k ; t) and moment functions Y (k) (a k ; t) as: Here (n is the Pochhammer symbol, and s(k, ), S(k, ) are Stirling numbers of the first and second kind [47,48].Although we are primarily interested in the functions Y (k) (a k ; t), the factorial moments X (k) (a k ; t) 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 (n) k and summing over all n ≥ k, we find where, for simplicity of notation, the arguments (a k ; t) have been suppressed from ρ (k) n and X (k) .In the case where the birth and death rates β n (a) = β(a) and µ n (a) = µ(a) are independent of the sample size, significant cancellation occurs and we obtain the elegant equation To find the boundary conditions associated with Eq. 24, we combine the definition of X (k) with the boundary condition in Eq.21 and obtain Setting X (0) ≡ 0, we recover the boundary condition associated with the classical McKendrick-von Foerster equation.For higher-order factorial moments, the full solution to the (k − 1) st factorial moment X (k−1) (a k−1 ; t) is required for the boundary condition to the k th moment X (k) (a k−1 , 0; t).Specifically, consider the second factorial moments and assume the solution X (1) ≡ Y (1) to the McKendrick-von Foerster equation is available (from e.g., Eq. 5).In the small interval da, the term Y (1) da is the Bernoulli variable for an individual having an age in the interval [a, a + da].Thus, in an extended age window Ω, we heuristically obtain the expectation where Y Ω (t) 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 find Thus, the second moment Y (2) 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 Y (2) and determining the variance of population-size-age-structured random variables.Eq. 24 for general k is readily solved using the method of characteristics where the propagator is defined as U (a, b) ≡ exp − b a µ(s)ds .We can now specify X (k) in terms of boundary conditions and initial conditions by selecting m = min {a k , t}.Since X (k) (a k ; t) ≡ X (k) (π(a k ); t) is invariant to any permutation π of its age arguments, we have only two conditions to consider.The initial condition X (k) (a k ; 0) = g(a k ) encodes the initial correlations between the ages of the founder individuals and is assumed to be given.From Eq. 22, X (k) (a k ; 0) must be a symmetric function in the age arguments.A boundary condition of the form X (k) (a k−1 , 0; t) ≡ B(a k−1 ; t) 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 (β(a) = β, µ(a) = 0) 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 A for more details): We have given X (2) (a, b; t) for only a < b since the region a > b can be found by imposing symmetry of the age arguments in X (2) .After using Eq.22 to convert X (1) and X (2) into Y (1) and Y (2) , 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 a → 0 and b → ∞, we recover the expected exponential growth of the total population size E(Y [0,∞] ) = e βt for a Yule-Furry process.We also recover the known total population variance Var(Y [0,∞] ) = e βt (e βt − 1).

Full Solution
Equation 17 defines a set of coupled linear integro-differential equations in terms of the density ρ n (a n ; t).
In [20], we derived explicit analytic expressions for ρ n (a n ; t) in the limits of pure death and pure birth.
Here, we outline the derivation of a formal expression for the full solution.It will prove useful to revert to the following representation for the density: If a n is restricted to the ordered region such that a 1 ≤ a 2 ≤ . . .≤ a n , f n can be interpreted as the probability density for age-ordered individuals (see [20] for more details).We will consider f n as a distribution over R n ; 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 f n (a n ; t) is used.
To find general solutions for f n (a n ; t) expressed in terms of an initial distribution f m (•; 0), we replace ρ n (a n ; t) with f n (a n ; t)/n! 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 a n of the individuals at time t.However, a n varies in time complicating the analytic expressions.If we index each characteristic by the time of birth (TOB) b = t − a 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 t − a n ≡ [t − a 1 , t − a 2 , . . ., t − a n ].The method of characteristics then solves Eq. 17 This equation is defined in terms of a propagator Ûn (b m ; t 0 ; t) ≡ U n (a m ; t 0 ; t) that represents the survival probability over the time interval [t 0 , t], for m individuals born at times b m , in a population of size n, where we have again used the definition γ n (a) = β n (a) + µ n (a).The propagator Û satisfies certain translational properties: The solution fn applies to any region of phase space where t 0 ≥ max(b n ).If t 0 = max(b n ), say t 0 = b n , then we must invoke the boundary conditions of Eq. 18 to replace fn , where we have used and will henceforth use the notation Eq. 34 is then used to propagate fn−1 (b n−1 ; b n ) 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.
Let b m ∈ [0, t] and b n < 0 denote the TOBs for m individuals born in the time interval [0, t], and n founder individuals, all alive at time t.Next, define y k ∈ [0, t] and y < 0 to be the TOBs of k individuals born in the time interval [0, t] and founder individuals, respectively.Here, all k + individuals are assumed to die in the time window [0, t].Their corresponding times of death are defined as s k and s , respectively.Thus, there will be n + individuals alive initially at time t = 0 and m + n 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 b 1 and y 1 .The individual with TOB b 1 survives until time t, while the individual with TOB y 1 dies at time s 1 .Within the time interval [0, t] there are three more births with TOBs b 1 , b 2 and y 1 , the last of which has a corresponding death time of s 1 , resulting in three individuals in total that exist at time t.
To express the distribution f3 (b 2 , b 1 ; t) in terms of the initial distribution f2 (b 1 , y 1 ; 0), conditional upon three birth and two death events ordered such that 0 < y 1 < s 1 < b 1 < b 2 < s 1 < t, we start with the distribution f2 (b 1 , y 1 ; 0).Just prior to the first birth time y 1 , we have two individuals, so that f3 (•; y − 1 ) ≡ 0 and Eq.34 yields f2 (b 1 , y 1 ; y − 1 ) = f2 (b 1 , y 1 ; 0) Û (b 1 , y 1 ; 0, y 1 ) (the death term does not contribute).To describe a birth at time y 1 , we use the boundary condition of Eq. 18 to construct f3 (b 1 , y 1 , y 1 ; ). Immediately after y 1 and before the next death occurs at time s 1 , three individuals exist and f2 (•; y + 1 ) ≡ 0. Now, only the death term in Eq. 34 contributes and Continuing 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: Û (y 1 ; 0, s 1 ) Û (b 1 ; 0, t) Û (y 1 ; y 1 , s 1 ) Û (b 1 ; b 1 , t) Û (b 2 ; b 2 , t).In other words, each birth-death pair (y, s) is propagated along the time interval it survives; from max{y, 0} to min{s, t}.For example, the individual with TOB b 1 < 0 survives across the entire timespan [0, t], whereas the individual with TOB y 1 is born and dies at times y 1 and s 1 .These two individuals are propagated by the terms U (b 1 ; 0, t) and U (y 1 ; y 1 , s 1 ), respectively.Provided the order 0 < y 1 < s 1 < b 1 < b 2 < s 1 < t is preserved and the values b 1 , y 1 < 0 are negative, the form of the integral expressions in Fig. 2   Fig. 3 Monte-Carlo simulations of densities in age-and capacity-dependent birth-death processes.Row A are results for a death-only process a linear death rate function µ(a) = a.We initiated all simulations from 10 individuals with initial age drawn from distribution P = 128a 3 e −4a /3.In row B, we consider the a buddingonly birth process with a carrying capacity K = 5 (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 ρ (0) n (t) = da ρn(a; t) for both processes.We also plot the single-particle density function ρn=1,5,9(a; t = 2) for the pure death process in A(ii-iv) and ρn=1,3,5(a; t = 5) for the limited budding process in B(ii-iv).Finally, the population-summed two-point correlations functions n ρ (2) n (a1, a2; t) for pure death and pure budding are shown in panels A(v) and B(v).
After summing across all realizations C m,k, (the configuration in Fig. 2 is one member of C 2,1,1 ) of the possible orderings of the birth and death times b m , y k , y , s k and s , we can write the general solution to Eq. 34 in the form The terms t − (x) and t + (x) refer to the times below and above x relative to the ordering of times b m , y k , y , s k and s k .For example, in Fig. 2 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 b n 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 b n 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 f n (a n ; t), lower dimensional marginal distributions such as f n (a 1 ; t) and f (2) n (a 1 , a 2 ; t), and their counterparts ρ n , can be sufficiently sampled.Figures 3A and 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 µ(a) = a and initiated the pure death process with 10 individuals with initial ages drawn from a gamma distribution with unit mean and standard deviation 1  2 .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 ρ (2) n (a 1 , a 2 ; t = 2) is shown in Fig. 3A(v).The observed structure indicates no correlations in the death only process and the peak at a 1 = a 2 ≈ 2.6 reflects the fact that older individuals die faster, shifting the mean age slightly below the initial age plus the elapsed time (1+2=3).Results from Monte-Carlo simulations of a pure birth process with growth rate β 0 = 1 and carrying capacity K = 5 (Eq.10).Here, we initiated the simulations with one individual with age drawn from the same gamma distribution P (a) = 128a 3 e −4a /3.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 n (a 1 , a 2 ; t = 5) 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.Similar analytic results can be obtained using Doi-Peliti second quantization methods, as is demonstrated in the companion paper [19].In particular, the age-structured population-size function ρ n (t) is expanded into a similar sum, where each term can be interpreted two ways: as an element in a perturbative expansion and also represented as a Feynman diagram in a path integral expansion.The moment equations from Section 3.1 that generalize the McKendrick equation can also be derived using second quantization.

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 our density functions.For example, in Fig. 1B at time t 1 , we have two pairs of distinct twins, with four individuals having two ages, whereas at time t 2 we have two singlets and two doublets.Thus, we define the ages of current singlets and twins by a m and a n , 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 x m = t − a m and y n = t−a n , respectively.For simplicity, we will assume that no simple birth processes occur and that particles grow in number only through fission.The function β m,n (a) 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 µ m,n (a), and event rate γ m,n (a) = β m,n (a)+µ m,n (a).We suppose, for the moment, that the TOBs are ordered so that x 1 ≤ x 2 ≤ . . .≤ x m and y 1 ≤ y 2 ≤ . . .≤ y m .The quantity f m,n (x m ; y n )dx m dy n is then the probability of m singlets with ordered TOBs in [x m , x m +dx m ] and n twin pairs with ordered TOBs in [y n , y n + dy n ].The density f m,n satisfies the following equation: where the partial age vectors are defined as x i,j = (x i , . . ., x j ) and the singlet age vector, doublet age vector, and time arguments are separated by semicolons.The term x . ., x m ) represents the vector of all m singlet TOBs, except for the i th 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 [x i , x i+1 ], while the second term describes the death of any one of two individuals in a pair of twins (with TOB x i ).
The associated boundary conditions are: The 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 x i .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 x m and y n be unordered TOB vectors, and extend f m,n to the domain R m+n by defining f m,n (x m ; y n ; t) = f m,n (T(x m ); T(y n ); t), where T is the ordering operator.Note that f m,n is not a probability distribution under this extension; however, ρ m,n (x m ; y n ; t)dx m dy n = 1 m!n! f m,n (x m ; y n ; t)dx m dy n 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 1, 2, . . ., m and the doublets 1, 2, . . ., n, the i th singlet has age in [x i , x i + dx i ] and the j th doublet have age in [x j , x j + dx j ].The density ρ m,n obeys with associated boundary condition Equations 44 and 45 provide a complete probabilistic description of the population of singlets and doublets undergoing fission and death.We draw attention to the parallel paper [19], where we derive an equivalent hierarchy using methods used in quantum field theory developed by Doi and Peliti [13,14,42].

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 ρ m,n (x m ; y n ; t) to construct marginal or "reduced" densities.Successive integrals over any number of the variables x m and y n can be performed, giving: For example, ρ (0,0) m,n (; ; t) 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 equations To solve Eqs.50 and 51, we first define and find solutions of the form X(x, t) = X(x, t 0 )U (x; t 0 , t) + 2Y (x, t 0 )U (x; t 0 , t)(1 − U (x; t 0 , t)), provided t 0 ≥ x.For an initial time of t = 0, we find, upon setting t 0 = max{0, x}, We now substitute Eqs.55 and 56 into Eqs.51 to find a Volterra equation for B: 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, T (x, t) = X(x, t) + 2Y (x, t), takes on a simple form in terms of B(t): while the total mean population T (t) = ∞ 0 T (x, t)dx is given by Before 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 Subsection 2.3) in the mean-field limit.

Mean-field Equivalence to the Bellman-Harris Process
Consider a Bellman-Harris fission process with a branching time distribution g(a) and a cumulative density function G(a) = a 0 g(a )da , along with the progeny distribution function A(•) given in Eq. 12.The Bellman-Harris model in Eq. 13 can be written in the form If we restrict ourselves to a binary fission process, the progeny distribution function takes the form A(y) = a 0 + a 2 y 2 , where a 0 and a 2 = 1 − a 0 are the death and binary fission probabilities, conditional on an event taking place.Thus, the mean population defined as has the Laplace-transformed solution We now show that the same result arises from our full noninteracting (population-independent β(a) and µ(a)) kinetic approach.Since the fission and death rates can be expressed as β(y) = a2g(y) 1−G(y) and µ(y) = a0g(y) 1−G(y) , Eq. 53 reduces to U (x; x, t) = 1 − G(t − x) and U (0; 0, t) = 1 − G(t).Starting from a single individual at age zero, Eq. 59 can be written as which has the Laplace-transformed solution Similarly, Eq. 57 becomes with Laplace-transformed solution Substituting Eq. 66 in Eq. 64 results in Eq. 62 for T (s), 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, β n (a) and µ n (a), and can also depend on population size, providing much greater flexibility.

A Cell Division Model
We now consider explicit calculations for a simple fission-only model of cell division where cell cycle times are rescaled to be Γ -distributed with unit mean and variance 1 α .This Γ -distribution and its Laplace transform g(s) are explicitly Equation 66 for B(t) can then be solved to yield The inverse Laplace transform is detailed in Appendix B and involves contour integration that yields (69) Similarly, from Eq. 62 we have which can also be evaluated via a similar Bromwich integral: For α = 1, g(t) = e −t is exponentially distributed, and we find the simple growth law T (t) = e t , which is equivalent to the result E(Y [0,∞] ) = e βt found earlier in Subsection 3.1.As α → ∞, the gammadistribution sharpens about unity and the process becomes more discrete-like in time.Figs.4A,B,C show that as α is increased, the mean population size T (t) tends towards that given by the discrete-time Galton-Watson step process, as would be expected.In Figs.4D,E,F, we have used the expression for B(t) in Eqs.58 and 69 to give the mean age-time distribution T (x, t).Note that unlike the solution to the Bellman-Harris equation shown in Figs.4A,B,C, the mean density T (x, t) (Eq.58) resolves age structure.

Spatial Models
We now illustrate how our kinetic (in age) theory can be generalized to include spatial motion such as diffusion and convection.We will follow the approaches described in Webb [53] 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, in this section we briefly outline how to generalize the age-structured spatial process discussed in [53] to incorporate stochasticity.
Consider a simple budding-mode birth-death process such that ρn (b n ; q n ; t) is the probability density for a population containing n randomly labelled individuals with TOBs b n and positions q n .Although ρn (b n ; q n ; t) is again invariant under permutations of the elements, the relative orders of b n and q n must be preserved: ρ2 (b 1 , b 2 ; q 1 , q 2 ; t) = ρ2 (b 2 , b 1 ; q 2 , q 1 ; t), for example.The TOB arguments in this density can be readily transformed to age (rather than TOB) by a transformation of the type given in Eq. 33.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 incorporated can be found in [53].We suppose that β n (a; q) and µ n (a; q) are birth and death rates for any 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 ρn (b n ; q n ; t): The corresponding boundary condition capturing the influx of newborn individuals is which 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 (q n ) as the parent (q i ), setting its identity as the i th individual.In addition to a boundary condition, Eq. 72 requires an initial condition ρ n (b n ; q n ; 0) to specify both the initial TOB and initial position of individuals.
As with our earlier analyses, we first express ρ n in terms of ρ n+1 by introducing the propagator ds , which enables us to transform Eq. 72 to an inhomogeneous heat equation for the function whose solution can be expressed as [5] ρ n (b n ; q n ; t) =U n (b n ; q n ; t 0 , t) Here, I n denotes the n × n identity matrix and N q (x, Σ) is the multivariate normal density for the vector q arising from a distribution with mean x and covariance Σ.This result expresses ρ n in terms of ρ n+1 and is analogous to Eq. 34.This solution is valid provided t 0 > max{x}; for t 0 = max{x}, 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.In the companion paper, we derive the mean-field equations for this spatial kinetic theory using quantum field theoretic methods developed by Doi and Peliti [19].

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 [20] which we extend here is the derivation of a kinetic theory for stochastic agestructured 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 additional specific mathematical results.Firstly, 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 populationindependent 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.
Secondly, 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.
Thirdly, we also outlined how to incorporate spatial dependence of birth and death into our agestructured 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.
All of our results can also be derived using from quantum field theoretical approaches [13,14,42], which are described in detail in a parallel paper [19].Such methods provide alternative machinery to analyze the statistics of age-and space-structured populations and may provide new avenues for calculation.
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 [25]) and would require extension of our state space to include generational class [56].These extensions will be explored in future work.involves a branch point at s = 0, we construct a branch cut along the negative real axis and define s = re iθ where θ ∈ (−π, π).The denominator s α − 2 also produces poles at s = 2 where n is an integer with |n| ≤ α 2 .The contour required for the Bromwich integral is shown in Fig. 5 and is evaluated using Cauchy's residue theorem.
The integrals around the outer perimeter and the origin contribute zero in the limit as R → ∞ and ε → 0. The branch cuts and poles provide the nonzero contributions.First, consider the integrals along the branch cut.Writing the variable s as re iθ , for θ = ±π, we integrate 1 Combining the contributions from the branch cut and the residues results in L −1 (t) 1 s α −2 , which, when substituted into Eq.68, gives the final result in Eq. 69.
The derivation for the Laplace inversion in Eq. 70 is similar.Note that the value s = 1 is a removable singularity and the same set of poles and paths for branch cut integrals needs to be considered.The details are left to the reader.

)
When k = 1, 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, β n and µ n , depend on the population size, one has to analyze the complicated hierarchy given in Eq. 23.

Fig. 2 A
Fig.2A sample birth death process over the time interval [0, t].Red and white circles indicate births and deaths within [0, t].The variables bi > 0 and b j < 0 denote TOBs of individuals present at time t, while yi > 0, y j < 0, and si, s j ∈ [0, t] 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.
, t − (b 2 ) = [s 1 , b 1 ] and t + (b 2 ) = [b 2 , s 1 ] represent the lower and upper bounds of the vector b 2 = [b 1 , b 2 ] found from the ordering 0 < y 1 < s 1 < b 1 < b 2 < s 1 .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.

Fig. 4
Fig. 4 Mean population and age-time distributions for a fission-only process with Γ -distributed branching times.A, B, and C show mean populations as a function of time for dispersion values α = 1, α = 10, and α = 100, 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.

Fig. 5
Fig.5Bromwich integral for calculating the inverse Laplace transform in Eq. 80.The integral along γ is evaluated using the residues at the poles and the integrals along the branch cut in Cauchy's theorem.