Kinetic theory of age-structured stochastic birth-death processes

Classical age-structured mass-action models such as the McKendrick-von Foerster equation have been extensively studied but are unable to describe stochastic fluctuations or population-size-dependent birth and death rates. Stochastic theories that treat semi-Markov age-dependent processes using, e.g., the Bellman-Harris equation do not resolve a population's age structure and are unable to quantify population-size dependencies. Conversely, current theories that include size-dependent population dynamics (e.g., mathematical models that include carrying capacity such as the logistic equation) cannot be easily extended to take into account age-dependent birth and death rates. In this paper, we present a systematic derivation of a new, fully stochastic kinetic theory for interacting age-structured populations. By defining multiparticle probability density functions, we derive a hierarchy of kinetic equations for the stochastic evolution of an aging population undergoing birth and death. We show that the fully stochastic age-dependent birth-death process precludes factorization of the corresponding probability densities, which then must be solved by using a Bogoliubov-–Born–-Green–-Kirkwood-–Yvon-like hierarchy. Explicit solutions are derived in three limits: no birth, no death, and steady state. These are then compared with their corresponding mean-field results. Our results generalize both deterministic models and existing master equation approaches by providing an intuitive and efficient way to simultaneously model age- and population-dependent stochastic dynamics applicable to the study of demography, stem cell dynamics, and disease evolution.

Figure 1

(a) A simple age-dependent birth-death process. Each parent gives birth with an age-dependent rate ${\beta }_n\left(a\right)$, which may also depend on the total population size $n$. Individuals can also die (open circles) at an age- and population-dependent rate ${\mu }_n\left(a\right)$. (b) Age trajectories in the upper ($a>t$) octant are connected to those in the lower one ($a<t$) through the birth processes. Individuals that exist at time $t=0$ can be traced back and defined by their time of birth $b_i$. Here the labeling is ordered according to increasing age. The pictured trajectories define characteristics $a_i\left(t\right)$ that can be used to solve Eq. (12).

Figure 2

Comparison of $P(a,t)$ [Eq. (26)] derived from the McKendrick–von Foerster equation with $P_m(a,t)$ of a fully stochastic pure birth process with constant $\beta =0.1$. We start with $N=10$ individuals and analyze our quantities at time $t=10$ for ages $a<t$. (a) Each of the 100 gray lines counts the number of individuals younger than age $a$ in one simulation. The solid black curve indicates the deterministic (McKendrick–von Foerster) solution $P(a,t)={\int }_0^a\rho (y,t)dy$, which can also be obtained through $P(a,t)={\sum }_{m=1}^{\infty }m{P}_m(a,t)$. The shaded region represents the inter-quartile range of $P_m(a,t)$, (the central count of individuals occupied by 50% of simulations). (b) Distribution constructed from 1000 simulations (bars) and theoretical distribution $P_m(a=5,t=10)$ (black curve).