Solving the stochastic dynamics of population growth

Authors: Loïc Marrec, Claudia Bank, Thibault Bertrand

Abstract

Population growth is a fundamental process in ecology and evolution. The population size dynamics during growth are often described by deterministic equations derived from kinetic models. Here, we simulate several population growth models and compare the size averaged over many stochastic realizations with the deterministic predictions. We show that these deterministic equations are generically bad predictors of the average stochastic population dynamics. Specifically, deterministic predictions overestimate the simulated population sizes, especially those of populations starting with a small number of individuals. Describing population growth as a stochastic birth process, we prove that the discrepancy between deterministic predictions and simulated data is due to unclosed‐moment dynamics. In other words, the deterministic approach does not consider the variability of birth times, which is particularly important with small population sizes. We show that some moment‐closure approximations describe the growth dynamics better than the deterministic prediction. However, they do not reduce the error satisfactorily and only apply to some population growth models. We explicitly solve the stochastic growth dynamics, and our solution applies to any population growth model. We show that our solution exactly quantifies the dynamics of a community composed of different strains and correctly predicts the fixation probability of a strain in a serial dilution experiment. Our work sets the foundations for a more faithful modeling of community and population dynamics. It will allow the development of new tools for a more accurate analysis of experimental and empirical results, including the inference of important growth parameters.

INTRODUCTION

Population growth is at the heart of fundamental processes in cell biology, evolution, and ecology, from the expansion of bacteria colonies and large‐scale animal populations to the propagation of an advantageous mutation. Predicting population growth dynamics is thus paramount in multiple fields. The outspoken goal of population modeling is to accurately describe the variation in the number of individuals in a population.

Historically, deterministic models are most commonly used to describe population dynamics (Anderson, 2014; Brauer & Castillo‐Chavez, 2013; Hannon et al., 1997). In these models, the population size is generically described by a continuous variable whose temporal dynamics are governed by an ordinary differential equation. Whereas most of these models are nonlinear—which means that analytical progress is not impossible but limited in some cases (Tsoularis & Wallace, 2002)—it is often relatively simple and computationally fast to obtain accurate numerical solutions, possibly explaining their widespread use. A paradigmatic model of population growth in theoretical ecology is the well‐known logistic equation whose study traces back to as early as the middle of the nineteenth century (Verhulst, 1838). The logistic differential equation was initially derived from introducing a self‐limiting property in the growth of a biological population to the unconstrained Malthusian exponential growth model (Malthus, 1798). It was rediscovered independently later on (Lotka, 1925; McKendrick & Pai, 1912; Pearl & Reed, 1920). Verhulst's logistic growth model derivation stemmed from the observation that unhindered exponential population growth is unrealistic. Even in the absence of predation relations, intraspecies competition for environmental resources such as food or habitat will lead to a characteristic saturation level, an upper bound on the population size known as the carrying capacity. Owing to its ease of use, the simplest logistic growth was used to model biological systems at all scales, from the population growth of micro‐organisms (Carlson, 1913; Pearl, 1927) to that of large mammal herds (Morgan, 1976) and fish schools (Krebs, 1985).

Further refinements to the logistic growth function led to the development of a generalized logistic growth model (Tsoularis & Wallace, 2002), which captures several commonly used population growth models including Blumberg (Blumberg, 1968), Richards (Richards, 1959) and Gompertz (Gompertz, 1825) growth models. Whereas amenable to easy progress and qualitative predictions, these deterministic models are not entirely faithful to the growth of a natural population, which is inherently stochastic (Bartlett, 1960; McKendrick, 1925). This stochasticity results from both intrinsic (e.g., demography) and extrinsic (e.g., environmental change) noise (Haefner, 2012; Lande et al., 2003). More recent studies have shown that deterministic and stochastic approaches yield critically different results (Allen & Burgin, 2000; Baker & Simpson, 2010; Wakano & Iwasa, 2013; Wilson, 1998).

Although it is often assumed that a deterministic equation can describe the dynamics of a large‐volume stochastic system, it is clear that this criterion alone is not sufficient (Gustafsson & Sternad, 2013; Kurtz, 1972). The range of validity of deterministic models is put in question. Even if new conditions for a deterministic equation to describe well the stochastic dynamics of a population have been outlined (Gustafsson & Sternad, 2013), they are not exhaustive and quantitative methods to overcome this discrepancy are missing. Recognizing these limitations, stochastic models have proved helpful in ecology for the past decades (Black & McKane, 2012; Bolker & Pacala, 1997; Keeling, 2000b; Marion et al., 1998; Matis & Kiffe, 1999).

Many studies use deterministic equations to fit experimental and empirical data, allowing the estimation of essential biological parameters. For example, logistic growth models have been used in microbiology to estimate microbial division rates (Kahm et al., 2010; Sprouffske & Wagner, 2016). A prediction based on deterministic models carries the risk of poorly estimating parameters of interest, such as the division rate of antimicrobial‐resistant bacteria, crucial to implementing political measures to slow the spread of antimicrobial resistance. Deterministic models were also used recently to estimate the basic reproductive number during viral outbreaks; those used logistic‐like equations (Aviv‐Sharon & Aharoni, 2020; Pelinovsky et al., 2020; Shen, 2020; Wu et al., 2020), as well as compartmental models such as SEIR (Shen et al., 2020; Sunhwa & Moran, 2020; Wan et al., 2020). Identifying when a deterministic equation does not correctly describe the average dynamics of stochastic population growth, understanding the reasons for this disagreement, and proposing solutions to remedy it, is thus of paramount importance.

In this work, we leap forward by solving the stochastic dynamics of population growth in the absence of deaths analytically. This resolution allows us to identify the extent to which a deterministic approach is a good approximation of the growth dynamics and to lay the foundation for future inference methods of growth parameters based, for example, on the likelihood function calculation. We consider several classical population growth models. First, we model the population growth as a stochastic birth process and simulate stochastic realizations of these kinetics. We compare their ensemble average to the predictions of the respective deterministic models. We show that the deterministic approach generically overestimates the average population size. This prediction error is larger when the initial number of individuals is very low. To explain the reason behind this discrepancy, we derive a master equation formalism describing the stochastic population growth dynamics and the moment equations. We find that the difference between the population size estimated by the deterministic equation and the mean of the simulated data is due to unclosed moment dynamics. We show that some moment‐closure approximations reduce the difference, although it remains globally substantial. Instead, we derive an exact solution to the stochastic population growth. Finally, we apply our results and show that our solution leads to a better prediction of the dynamics of two competing strains and the probability of fixation of a mutant in a serial passage experiment.

BIAS OF DETERMINISTIC APPROACHES

Pure‐birth models

Given the inherent stochasticity of population growth processes, we first establish whether a deterministic equation correctly describes the mean trajectory of stochastic growth. We consider four distinct growth models belonging to the generalized logistic growth models: Blumberg, Gompertz, Logistic, and Richards models (Tsoularis & Wallace, 2002). Our choice was motivated by their widespread use to fit experimental or empirical data to estimate growth parameters in microbiology and ecological communities (Ghenu et al., 2022; López et al., 2004; Sprouffske & Wagner, 2016). These kinetic models differ by their per capita growth rates, bN.

Under Malthusian growth, the per capita growth rate is constant and independent of the population size N; we denote as b this intrinsic birth rate—also called exponential birth rate or Malthusian parameter. Note that simple unbounded exponential growth would occur if no restrictions were imposed on the population size (e.g., nutrients, available space). In the growth models considered here, the per capita growth rate is explicitly dependent on the population size. To model environmental constraints such as availability of space or food, one then generically introduces a carrying capacity K that limits the population size N to the range N0≤N≤K assuming no deaths, where N0 is the initial population size. Specifically, the per capita birth rate bN decreases as the population size increases and vanishes when N=K.

For a general Logistic growth model, the deterministic equation describing the dynamics of N reads (1) dNdt=BN=bNα1−NKβγ, where BN denotes the population growth rate. The population growth rates for each of the four nonlinear models considered here are provided in Table 1 along with those of the exponential growth model.

As seen in Figure 1, the four population growth models chosen here display very different birth rate curves. Birth rates in all models are non‐monotonic and vanish when N→0 and N→K by construction. In other words, these models impose that a population of size zero cannot grow, and no population can grow beyond the carrying capacity. We note that the Logistic model displays a symmetry around the population size N=K/2, whereas in both the Blumberg and Richards models, the exponents β and γ, respectively, offer an extra degree of freedom to tune the shape of the growth rate curve and in particular, its asymmetry. We note that the population size at inflection, that is, when the population growth rate is maximum, is given by (2) N*=K1+βγα−1/β, and is thus dependent on exponents β and γ for the Blumberg and Richards models, respectively. The Gompertz model shows the fastest growth of all at small population sizes. As shown in Figure 2a, the population size in all deterministic models follows a sigmoid curve—also called an S‐shape curve—with its characteristic initial phase with slow growth, exponential growth phase, finally followed by a stabilization of the population size at a finite steady‐state population size (3) limt→∞Nt=K.

For the Logistic, Richards, and Gompertz models, full analytical solutions are available (Tsoularis & Wallace, 2002). These are given by (4a) NLt=K1+K/N0−1e−bt,

(4b) NRt=K1+K/N0β−1e−βbt1/β,

(4c) NGt=KN0/Ke−bt.

Those are represented in Figure 2a. Note that no analytical (i.e., closed‐form) solution exists for the Blumberg model; in this last case, we proceeded with a direct numerical solution of Equation (1).

However, as pointed out earlier, population growth is inherently stochastic. We, therefore, evaluated the validity range of the above deterministic descriptions by simulating individual stochastic trajectories for the four growth models introduced above using a Gillespie algorithm (Gillespie, 1976, 1977) (see Appendix A). To this end, we recast our problem into a pure‐birth process for which population growth results from an individual A reproducing at a size‐dependent rate bN following the elementary reaction (5) A→bN2A, with BN=bN×N. Naturally, multiple stochastic models may lead to the same deterministic model under a mean‐field approximation. Here, we focus on one particular microscopic scenario. Still, as we will argue in the next section, other formulations, including those where birth and death processes are taken into account explicitly, lead to even more drastic disagreement. We average all our results over 10^5^ independent stochastic trajectories to obtain the time‐dependent average population size. Note that we only consider the case of α=1 to ensure that we can faithfully match the stochastic models to their deterministic limits. Indeed, these are the only cases for which we obtain a well‐defined deterministic limit (see Appendix B for details).

As shown in Figure 2a, we observe a substantial difference between the deterministic predictions and the stochastic mean population size starting from a single individual. To quantify this disagreement, we calculated the relative error η defined as (6) ηt=Ndt−NstNst, where Ndt and Nst are, respectively, the time‐dependent population sizes predicted by the deterministic model and measured in our stochastic simulations. We observed relative errors as large as 30% independent of the carrying capacity (see Appendix D, Figure A2). Figure 2b shows that the largest error for the parameters chosen is obtained for the Richards model, whereas the smallest is for the Gompertz model (still at a substantial maxη≈12%). We note that the quantitative value of the relative errors measured for the Blumberg and Richards models depend on the choice of exponents β and γ, respectively. Nevertheless, the measured errors remain substantial over a wide range of exponents (see Appendix C, Figure A1). For all the growth models studied here, the deterministic dynamics given by Equation (1) overestimate the population size at all times. It is interesting to note that the maximum error is located around the inflection point (i.e., at t such that d2N/dt2=0) predicted by the deterministic equation. The inflection point corresponds to the population size where the population birth rate reaches its maximum and starts decreasing.

Furthermore, we observe that the error uniformly decreases as N0, the initial number of individuals in the population, increases (see Figure 2c). This discrepancy limits the range of validity of the deterministic models as the initial number of individuals in the population is assumed to be small in many applications, for example, patient zero in a disease spreading scenario, single cell mutation in a mutation fixation experiment, small number of cells in a bacterial colony expansion, etc. To summarize, the deterministic equation fails to describe the average stochastic trajectory. Although the relative error depends on the specific growth model, and thus on the per capita birth rate, it remains substantial in all cases tested and increases with decreasing initial population size. Importantly, the relative error between deterministic predictions and measured mean stochastic population size is independent of the carrying capacity. Thus the discrepancy does not vanish in the limit of large but finite population sizes (see Appendix D, Figures A2 and A3).

Stochastic population growth is a Markovian jump process. Therefore, the times between jumps from size N to N+1 are exponentially distributed with parameter BN. For small initial population sizes N0, the rates at which the population initially grows are low (see Figure 1), e.g., the rate of the first reproduction is given by bN01−N0/Kβγ for the generalized logistic model. In turn, this implies that early in the process, the distributions of reproduction times lead to a large variance in the population size. We postulate that this large variance accumulated over the growth process is responsible for the disagreement between the deterministic and the mean stochastic trajectories. This postulate is consistent with our observations (1) the relative error increases when the initial population size decreases, and (2) the relative error is maximal around the inflection point where the exponential distribution of reproduction times is the tightest. Interestingly, when K→∞ (large volume limit), the rate of first reproduction converges to bN0 and so is entirely controlled by the initial population size confirming that the observed disagreement remains valid in this limit. Finally, we note that the above deterministic models are mean‐field models which intrinsically assume an underlying population size distribution peaked around its mean in contrast to the wide population size distributions observed in the stochastic models.

Birth‐death processes

In the previous section, we focused on populations that can only increase in size over time. This assumption, which entails neglecting the death of individuals, is predominant in microbiology, where the models used to fit population growth data do not explicitly include death rates (Kahm et al., 2010; Ram et al., 2019; Sprouffske & Wagner, 2016). Similarly, pharmacodynamic models, which aim at quantifying how antibiotics inhibit growth or kill cells, commonly replace the birth rate with a net birth rate (i.e., birth rate minus death rate) in Equation (1) (Czock & Keller, 2007; Regoes et al., 2004). In this way, the population grows if the net birth rate is positive, decreases if it is negative, and remains constant for zero net growth rates. However, stochastic population growth can also be modeled as a Markovian jump process where births and deaths are distinct events leading to distinct changes in population N→N+1 and N→N−1, respectively.

We do not expect deterministic models to fare better in the presence of explicit death events with small initial population sizes. For the sake of simplicity and without loss of generality, we add to Equation (1) a linear death term leading to the modified differential equation dN/dt=bNα1−N/Kβγ−dN. We further simulate a stochastic birth‐death process known to lead to this deterministic equation in the mean‐field limit using once again a Gillespie algorithm (Gillespie, 1976, 1977). We show that a deterministic approach fails to describe the dynamics of the average population growth over many stochastic realizations, as in the death‐free case. Strikingly, we demonstrate that it also fails to predict the correct steady‐state population size.

Indeed, Figure 3a shows that deterministic models do not predict quantitatively either time‐dependent population sizes or steady‐state population sizes averaged over many stochastic realizations. We argue that this difference is due to rapid extinctions, which are frequent occurrences when considering low initial population sizes and large ratios of death rate to birth rate (see Figures 3b,c). Specifically, demographic stochasticity leads to extinction with probability p0=d/bN0 (or p0=d/blogKN0 for the Gompertz model), which is non‐zero as long as the death rate is strictly positive. Note that the probability of extinction is obtained from the linear birth‐death process (Kendall, 1948), assuming density independence at the beginning of the growth when the population size is very small compared to the carrying capacity. For nonlinear models, these early extinction events are not taken into account in deterministic approaches.

One may wonder whether a fairer approach to matching the stochastic models with the deterministic limits should be to first condition the master equation on survival. One can show that this does not lead to a substantial reduction in the error in predicting the transient population size but allows one to correctly predict the quasi‐steady‐state population size (see Appendix D for details). Furthermore, the linear birth‐death process, whose analytical solution is known (Kendall, 1948), has a deterministic limit giving the population size averaged over the survival and early extinction trajectories (see Appendix D for details).

In summary, we conclude that deterministic formalism is not a good predictor of the average population growth dynamics, even for large carrying capacities in the presence or the absence of explicit death events. We also note that the discrepancy between deterministic and average stochastic population sizes worsens as the initial population size decreases. We further conclude that deterministic approaches fail at predicting the steady‐state population size when death events are explicitly introduced. In the following, we focus on the pure‐birth process, which is already of great interest in microbiology, as we pointed out earlier. In the next section, we adopt a stochastic approach to describe the population growth and obtain an exact analytical solution for the population size distribution at all times.

ERROR REDUCTION BY MOMENT‐CLOSURE APPROXIMATIONS

To identify the reasons behind the poor performance of the deterministic equation, we return to a stochastic formalism. Generically, any population growth in the absence of death may be described by a stochastic birth process whose rates are defined by the model (Kendall, 1949) (see Table 1 for examples). Let us consider a population whose number of individuals at time t is denoted by N, whereas its initial population size is N0. As in Equation (5), we consider that each individual in the population replicates with the same per capita rate bN. Here, the population size increases from N to N+1 individuals at a total rate BN, where BN was defined in Table 1 for several population growth models.

We focus on finite‐sized populations that grow in a constant environment with a carrying capacity K. To fully account for the stochasticity inherent to demographic noise, we use a microscopic and probabilistic description in continuous time of the birth events within the population. More specifically, we write a system of differential equations describing the probabilities PN0,Nt that a population has a given size N at a given time t knowing that it started with N0 individuals. Because the growth rates vanish when N→K, the size of our population is at most K, with the state N=K being an absorbing state, that is, the population indefinitely remains in this state once reaching it for the first time. Put simply, our stochastic process, while continuous in time, has a finite discrete number of possible states. Here, we assume that the population jumps between successive sizes with a rate dependent on its current size leading to a fully coupled system of equations.

This system of differential equations, formally called the master equation, governs the time evolution of the probabilities PN0,Nt. Writing the master equation for a stochastic jump process requires one to think about gain and loss terms to the probabilities PN0,Nt; for our system, it reads (Gardiner, 2009; Van Kampen, 2011) (7) dPN0,Ntdt=BN−1PN0,N−1t−BNPN0,Nt, when N0≤N≤K−1 and the probability normalization condition imposes PN0,Nt=1−∑N=N0K−1PN0,Nt, when N=K. The first term on the right‐hand side of Equation (7) is a gain term corresponding to an increase in the population size from N−1 to N individuals via a birth event, whereas the second term is a loss term corresponding to the population size transitioning from N to N+1 individuals. Although writing down the master equation for a stochastic jump process is often easy, computing the formal solution of master equations is arduous and has been an active field of investigation for decades (Gardiner, 2009; Van Kampen, 2011; Weber & Frey, 2017).

Importantly, the master Equation (7) contains all information about the growth dynamics; in particular, as it governs the probability distribution PN0,Nt, it contains all information to compute the averaged population size trajectory over time. Rather than solving directly Equation (7), we derive an equation governing the moments of PN0,Nt, where the mth moment of the population size is defined as (8) Nm=∑N=N0KNmPN0,Nt, with m a positive integer. For instance, the equation governing the time‐evolution of the first moment (i.e., m=1), which corresponds to the mean population size, reads (9) dNdt=BN, where as the equation for the second moment (i.e., m=2) satisfies (10) dN2dt=2N+1BN, where again averages are defined as ⋅=∑N=N0K⋅PN0,Nt.

The moment equations are closed for linear population birth rates as in the Malthusian (or exponential) growth model. Taking for the sake of simplicity and without loss of generality b=1, the population birth rate of the exponential model is given by BN=N. Writing out the moment equations, we obtain, for instance, (11a) dNdt=N,

(11b) dN2dt=N+2N2,

(11c) dN3dt=N+3N2+3N3.

First, the average population size given by Equation (11a) shows that the deterministic model accurately describes the average stochastic population size dynamics. Secondly, the equation for the second moment (11b) is closed as it only depends on the second moment itself and the first moment, which can be obtained analytically by solving Equation (11a). Similarly, the third moment (11c) is closed as well; the equation for the mth moment only depends on the m first moments, and so by solving the moment equations sequentially, we obtain a closed equation for any moment of the distribution. This way, important distribution properties such as its variance or skewness can be studied analytically. Note that the master equation (7) is analytically solvable for linear population birth rates, also called Malthusian growth (Taylor et al., 1998), which allows one to obtain distribution properties without solving the ordinary differential equations satisfied by the moments.

Meanwhile, the moment dynamics are unclosed for nonlinear population birth rates, which is often the case for finite populations. In other words, the equation for mth moment may involve higher‐order moments, leading to an infinite hierarchy of moment equations that is not solvable. Approximation techniques to get around this problem exist, and here, we discuss the accuracy of the most commonly used.

The most basic way to get around this problem is to apply a so‐called mean‐field approximation. The mean‐field approximation relies on the crucial assumption that the distribution of populations sizes is well‐peaked so that PN0,N≈δN−N, where δ is the Kronecker‐delta function (i.e., δN−N=1 if N=N and 0 otherwise). This approximation naturally leads to the approximations N≈N, and BN≈BN, which enable us to recover the deterministic limit given by Equation (1). More formally, a Kramers‐Moyal expansion in combination with a diffusion approximation of the master equation (7) and a mean‐field approximation leads to the same resulting deterministic equation (Gardiner, 2009) (see Appendix B for details). As we argued earlier, this deterministic model has been very popular, for instance, in population genetics theory focusing on the evolutionary effects of natural selection and mutation, as it simplifies calculations and circumvents the need for the master equation framework (Crow & Kimura, 1970; Ewens, 2004; Kingman, 1979; Nagylaki, 2013; Nei, 1987; Rice, 2004). Importantly, we argued in Section 2 that the distributions of population sizes were wide (see also Figure 5). So unsurprisingly, the mean‐field approximation fails to satisfactorily describe the population size averaged over stochastic realizations.

Going beyond the mean‐field approximation requires us to close the hierarchy of moment equations; these methods are called moment‐closure approximations. They have been extensively used to provide analytic approximations to nonlinear stochastic population growth models (Krishnarajah et al., 2005; Nåsell, 2003b; Singh & Hespanha, 2007). In the following, we focus on the Logistic growth model and proceed to several common moment‐closure approximations. Writing the first few moment equations for the Logistic growth model, we obtain (12a) dNdt=N−N2K,

(12b) dN2dt=N+2−1KN2−2KN3,

(12c) dN3dt=N−3−1KN2+3−3KN3−3KN4.

Notably, the equation for the first moment N depends on the second moment N2, whereas the equation for the second moment N2 depends on the third moment N3, and so on. To close this hierarchy of moment equations, two routes are often (i) closure methods can rely on a cumulant truncation procedure, in which the k first cumulant equations are approximated by setting all cumulants of order higher than k to 0 (Matis & Kiffe, 1996, 1999), (ii) closure methods can also be based on assumptions on the form of the underlying distribution of population sizes PN0,Nt (Isham, 1991; Krishnarajah et al., 2005; Marion et al., 1998; Nåsell, 2003b; Whittle, 1957). Most recently, the latter method has been extensively used; in these latter approximations, one often focuses only on the first two moments.

Here, we test several common moment‐closure approximations that express N3 as a function of the first two moments N and N2, allowing us to close the first two moment equations (12a) and (12b). We report all moment‐closure approximations tested here in Table 2. As shown in Figure 4, all moment‐closure approximations tested here show a substantial disagreement with the simulated stochastic mean population. Whereas the Binomial, separable derivative‐matching (SDM), and mean‐field approximations overestimate the mean population size, the New‐Poisson, Nåsell‐Poisson, and Normal moment‐closure approximations underestimate it. We report absolute relative errors ranging from ∼6% for the Binomial approximation to ∼26% for the normal approximation. Indeed, classical moment‐closure approximations fail in problems with very skewed underlying probability distributions, for which accurate knowledge of the higher order moments is needed (Krishnarajah et al., 2005). For instance, the third moment of the distribution (the one being approximated in the methods presented here) is directly related with the lopsidedness of the distribution. Figure 5 shows that our population growth dynamics are plagued by large skewness. In other words, Figure 5 shows that the population size can have many different values at a given time during stochastic growth.

Notably, some moment‐closure approximations perform better than the mean‐field approximation, although the latter is the most widely used. While these moment‐closure approximations generally reduce the error in predicting the dynamics of population growth, we note that they cannot generically be applied to the other growth models considered in this paper. Indeed, the right‐hand side of Equation (9) depends on the terms 1−N/KγN, logNN and Nβ+1 for the Blumberg, Gompertz, and Richards models, respectively. These are not easily expressed in terms of combinations of moments of the distribution. Since these nonlinear population growth models are widely used, we must devise a generically applicable method, which we do in what follows. However, we recognize that moment‐closure approximations allow growth dynamics to be described from one or two ordinary differential equations, which may be easier than a fully probabilistic description.

TOWARD AN EXACT SOLUTION TO STOCHASTIC POPULATION GROWTH

First Transition rate matrix

A first approach to try and solve the master equation directly is to recast it in the language of Markov chains. Namely, stochastic population growth can be interpreted as a continuous‐time, K‐state Markov jump process with transition rate matrix R=Rij1≤i≤K,1≤j≤K, where the matrix element Rij with i≠j is understood to be the rate at which the population switches from value N=j to N=i (note that other conventions exist). The diagonal elements of the transition rate matrix are generically fixed by enforcing conservation of total probability, which implies that the columns of R sum up to zero, that is, ∑i=1KRij=0, so that Rjj=−∑i=1,i≠jKRij.

As the population sizes potentially vary from 1 to K, we first define Pt as the column vector of dimension K

(13) Pt=PN0,1tPN0,2t⋮PN0,Kt, with the ith‐component corresponding to the time‐dependent probability of having a population of size i. We then rewrite Equation (7) in matrix form as (14) dPtdt=R⋅Pt, where the K×K transition rate matrix R is defined as (15) Rij=−Bjifi=j,Bjifi=j+1,0otherwise. leading to (16) R=−B100…0B1−B20…00B2−B3…0⋮⋮⋮⋮⋮0000−BK

In the case where the transition rate matrix does not have repeated entries, the solution to Equation (14) can be written as a weighted superposition of the K eigenvectors vk of the transition rate matrix multiplied by an exponential function whose rate is given by the associated eigenvalue μk. Namely, we write (17) Pt=∑k=1Kvkckeμkt, a vector whose components give us the time‐dependent probability that the population size is N given that it was N0 at t=0, (18) PN0,Nt=∑k=1Kvk,Nckeμkt.

Note that the coefficients ck are obtained by the imposition of the initial conditions and here must satisfy PN0,N0=∑k=1Kvk,Nck=δN,N0, where δ is the Kronecker delta (i.e., δN,N0=1 if N=N0 and 0 otherwise). Furthermore, ck=0 for 1≤N<N0 since a population size lower than N0 cannot be reached. Indeed, death events are not considered here, so the population size can only increase over time.

As the transition rate matrix R is lower triangular, the eigenvalues are equal to the diagonal entries of the matrix, and we obtain μk=−Bk, 1≤k≤K. Finally, we compute the eigenvectors as (19) vk,N=0if1≤N≤k−1,∏q=kN−1BqBq+1−Bkotherwise.

Notably, the eigenvectors vk are ill‐defined when the birth rates degenerate (i.e., when Bk=Bk′ for k≠k′). Indeed, Equation (19) shows that some of the components of the eigenvectors diverge when birth rates degenerate. The degeneracy typically happens when the population birth rate curve displays particular symmetries (see Figure 1). For instance, the Logistic growth model presents a mirror symmetry with respect to K/2; for all 1≤k≤K, we obtain Bk=BK−k. Diagonalization of the transition rate matrix would thus not apply to the Logistic growth model.

In the case where every birth rate is distinct (i.e., Bk≠Bk′ if k≠k′), we compute the constants ck, and we finally get (20) PN0,Nt=1BN∑k=N0NBke−Bkt∏q=N0q≠kNBqBq−Bk.

Note that the condition on the transition rate matrix not having repeated entries translates to imposing on the growth model not to lead to degenerate birth rates. The solution to the stochastic dynamics of population growth models with degenerate birth rates cannot be directly obtained via Equation (17); indeed the degeneracy in the birth rates leads to degenerate eigenvalues μii∈1,M, with M<K. This means that the solution would then be written as a linear combination of exponential terms with parameter μi, whose coefficients are polynomial factors. In general, identifying these polynomial factors can be a very tedious task. To circumvent this issue, and thus the lack of universality of Equation (17), we develop in the next section a new resolution method that applies to any population growth model.

Exact distribution of waiting times between birth events

As we just saw, the method based on the transition rate matrix is unsatisfactory because it does not yield a closed‐form solution for Equation (14) if birth rates are degenerate. Here, we suggest an approach based on waiting times between birth events. The underlying idea is that in the absence of deaths, population growth can be interpreted as a succession of events (births) happening in a well‐defined order separated by waiting times, which are random variables. To reach size K from its initial size N0, the population has to grow one individual at a time and go from N0 to N0+1, then from N0+1 to N0+2, etc. in this precise order. Therefore, all the information needed to derive the time‐dependent probability PN0,Nt should be contained in the distribution of the time between two birth events. In other words, in the approach based on the master equation and the transition rate matrix, the reasoning is based on population sizes. In contrast, in this approach, our reasoning is based on waiting times between successive events.

We denote by τN the time elapsed between two births where the population size increases from N to N+1 individuals. Owing to the Markovian nature of the process, τN is a stochastic variable exponentially distributed with mean BN−1. Then, the probability of having a given number N of individuals at time t must be equal to the probability that N−N0 births occurred by t and not N−N0+1 yet. Quantitatively speaking, (24) PN0,Nt=Prob∑k=N0N−1τk<t<∑k=N0Nτk.

The sum of n exponentially distributed random variables with rates λi1≤i≤n follows a hypoexponential distribution (Amari & Misra, 1997; Ross, 2007), which we denote Htλi1≤i≤n. Hypoexponential distributions were previously studied in the context of population genetics (Strimmer & Pybus, 2001) but also cell and systems biology (Chao et al., 2019; Gavagnin et al., 2019; Golubev, 2016; Yates et al., 2017).

Using the expression for the probability density function for the hypoexponential distribution (see Table 3 for the three cases to consider), we write the exact solution to Equation (7) in the form (25) PN0,Nt=1BNHtBN0BN0+1…BN.

Notably, in the case where all population growth rates are distinct, the hypoexponential distribution takes the form (22), and we recover exactly the solution introduced in the previous section (see Equation (20)). In Figure 5, we compare our exact solution (25) to the mean population size measured in simulations of the stochastic process for our four nonlinear growth models and show perfect agreement in all cases. We also confirm in Figure 5 that the full probability distributions measured from simulations agree with our exact solution.

APPLICATIONS

Finally, we cover two examples of applications in which we show that using a deterministic model instead of an exact solution to stochastic population growth leads to quantitatively very different results and may misinterpret important experimental results.

Population growth dynamics within a community

First, we extend our results to the study of a community composed of multiple strains evolving in the same environment. Community dynamics are at the heart of theoretical ecology and have been recently applied to the growing field of microbiome studies, where predicting the relative abundance of each microbial strain in the gut microbiota represents an opportunity for medical diagnosis and treatment (Gebrayel et al., 2022). For simplicity and without loss of generality, we focus on the case of two competing strains. Consider, for instance, the population dynamics of a wild‐type (W) strain and a mutant strain (M) competing in a batch culture environment; we denote their intrinsic birth rates bW and bM, respectively. As before, we define N as the size of the community, whereas n (resp. N−n) denotes the number of M (resp. W) individuals. We assume that the size of the community is limited by a single carrying capacity K. Note that our approach is easily generalizable to cases with multiple strains or with different carrying capacities and with explicit interaction parameters.

Furthermore, we introduce the relative fitness of the two strains, defined as r=bM/bW; this ratio indicates which strain is favored by natural selection. Specifically, if r>1 (resp. r<1), then strain M is beneficial (resp. deleterious), with r=1 corresponding to the neutral case. Each time an individual reproduces, the probability that this individual is of strain M is then given by (26) ΓN,n=rnr−1n+N.

Here, we start with an initial community size N0, composed of n0 individuals from strain M and N0−n0 individuals from strain W. The probability PN0,Nt that the community has a total size N at time t, knowing that the initial size of the community was N0, is given directly by Equation (25) with population reproduction rates (27) BN=∑n=n0N−N0+n0PN,n∣N0,n0BN,nM+BN,nW, where BN,nM and BN,nW are the rates at which each population increases. For the Logistic model, these rates are for instance given by (28a) BN,nM=bM1−N/Kn,

(28b) BN,nW=bW1−N/KN−n.

Furthermore, the probability PN,n∣N0,n0 of finding n individuals of type M when the total number of individuals is N satisfies (29) PN+1,n∣N0,n0=1−ΓN,nPN,n∣N0,n0+ΓN,n−1PN,n−1∣N0,n0, subject to the initial conditions PN0nN0n0=δn,n0. Note that Equation (29) has been extensively studied in (Houchmandzadeh, 2018). By definition, the probability to observe n individuals of type M at time t, knowing that we had initially n0 such individuals, is (30) Pn0,nt=∑N=N0KPN0,NtPN,n∣N0,n0, which we exactly compute using Equations (25), and (29). Then, using Equations (25), (27), and (29), we compute the average stochastic community and population sizes (31a) Nt=∑N=N0KNPN0,Nt,

(31b) nt=∑n=n0KnPn0,nt.

On the contrary, a deterministic description of the community dynamics leads to the system of differential equations (32a) dndt=bMn1−NK,

(32b) dNdt=bMn+bWN−n1−NK.

We once again compare the results of stochastic simulations quantitatively to the predictions of Equation (31) on the one hand and Equation (32) on the other hand. We show in Figure 6a,b that the deterministic model grossly overestimates the size of the community, whereas our stochastic solution perfectly matches the simulated data. Strikingly, the deterministic model is shown to overestimate the equilibrium population size of strain M. Our stochastic approach provides an exact prediction of the average community and population sizes and, more importantly, yields the full time‐dependent probability distributions of community and population sizes, which is not possible with a deterministic approach (see Figure 6c,d).

Fixation probability in a serial passage experiment

Finally, we show that our exact stochastic solution yields the fixation probability of a strain in a serial passage experiment. For this, let us assume that we start an experiment with the same number of individuals M and W. As usual, the initial size of the community formed by strains M and W is given by N0; we introduce the dilution rate D defined as the ratio of the initial size of the community to the carrying capacity, D=N0/K.

In a serial passage experiment, the community grows for a time τ before one applies a bottleneck by taking a random sample of N0 individuals; mathematically, this corresponds to selecting N0 individuals from the community following a binomial law. One then proceeds with a new growth phase of length τ before applying a new bottleneck. This process is repeated until only a single strain is left in the community. In these experiments, a quantity of interest is the probability pfix that the strain M fixes and that the strain W goes extinct. In particular, optimizing the fixation probability as a function of the dilution ratio D and the waiting time τ has attracted much attention, especially in the context of directed evolution (LeClair & Wahl, 2017; Wahl et al., 2002; Wahl & Gerrish, 2001). To our knowledge, existing studies answered this question using only deterministic approaches for modeling population dynamics.

To calculate pfix, we model the system as a Markov chain on the number of individuals M after each bottleneck event. At the moment of applying a bottleneck, the population contains Nτ individuals, including nτ individuals of strain M. If picking a single individual randomly from the community, the probability that this individual is of type M is thus given by nτ/Nτ. The probability Πk→l that the number of individuals M goes from k to l when a bottleneck is applied follows the binomial distribution (33) Πk→l=N0lnτNτl1−nτNτN0−l, where nτ depends on its starting value k. Note that here we follow the convention used in (Wahl et al., 2002; Wahl & Gerrish, 2001) and model bottlenecks as draws from a population with replacement leading to binomial distributions for Πk→l; one could also model bottlenecks as draws without replacement, naturally leading to hypergeometric distributions. Under a deterministic approach, nτ and Nτ are obtained by solving Equation (32) with the initial conditions n0=k and N0=DK. However, in a stochastic approach, we write the Πk→l as a sum over all possible pairs (N, n) at time τ weighted by their respective probabilities, (34) Πk→l=∑N=N0K∑n=kNPN0,NτPN,n∣N0,k×N0lnNl1−nNN0−l, where PN,n∣N0,n0 is governed by Equation (29). Finally, we note that conservation of probabilities imposes that Πl→l=1−∑k≠lΠk→l.

We define P=P0P1…PN0 as the column vector of probabilities Pi to have i individuals of strain M in the random sample of N0 individuals from the community following a bottleneck event. As such, the serial passage experiment defines a discrete‐time Markov process, in which time is measured in units of τ with t thus corresponding to the number of bottlenecks. Here we follow reference Marrec and Bitbol (2018) and take the limit of continuous time to write the master equation governing P as (35) dPtdt=R⋅Pt.

Here, the elements of the transition rate matrix R are given by (36a) Rij=Πj→i,

(36b) Rjj=−∑i≠jΠj→i, in which (36b) ensures conservation of probability.

The Markov process thus defined possesses two absorbing states, namely n=0 and n=N0, which correspond to the extinction and fixation of strain M, respectively. By definition, once the system reaches one of these states, it remains there indefinitely. Mathematically, this implies that the first and last columns of the transition rate matrix are filled with zeros as these columns contain the transition rates out of the n=0 and n=N0 states, respectively. As the matrix R is not invertible, we introduce the reduced transition rate matrix R˜ in which the rows and the columns corresponding to the absorbing states are removed. According to this definition, R˜ is invertible and we can write the fixation probability as (37) pfix=−∑i=1N0−1Πi→N0R˜−1i,n0, where n0 is the number of individuals M at the beginning of the experiment. As shown in Figure 7, the deterministic approach grossly overestimates the fixation probability for all dilution rates D. Here, we confirm this result for all nonlinear growth models studied. In contrast, we show that our stochastic solution exactly matches the results of stochastic simulations. This result has substantial consequences as we show in particular that the dilution ratio predicted by a deterministic approach to optimize the fixation probability is far from being the actual optimal dilution ratio, although it is commonly used in the literature, as was already argued.

DISCUSSION

In conclusion, we showed that a deterministic approach to population growth leads to biased predictions of the average behavior of this inherently stochastic process. More precisely, deterministic models overestimate the population size averaged over large numbers of stochastic realizations. This overestimation increases with decreasing initial population size. Qualitatively, the bias of the deterministic approach is due to the variability of the waiting times between reproduction events, which is particularly important at small population sizes. Quantitatively, the bias of the deterministic approach is due to unclosed‐moment dynamics. Importantly, we showed that moment‐closure approximations are insufficient to substantially reduce the relative difference between analytical predictions and average population sizes and do not apply to all population growth models.

It should be noted that methods other than moment‐closure approximations have been used to describe the average stochastic population dynamics by a simple ordinary differential equation. For example, the results of Kurtz (1972) on the convergence of average stochastic trajectories to a deterministic model in the case of large carrying capacity and large initial numbers of individuals have been extended more recently in this context of the so‐called fluid limit. This limit theorem shows that the dynamics of a population with a small initial size compared to the carrying capacity can be accurately described by a deterministic equation with a random initial size (Baker et al., 2018; Barbour et al., 2016). However, it is important to note that while Barbour et al. (2016) explicitly show that the whole distribution of the population size is converging to that of the fluid limit with random initial conditions, this often involves a non‐explicit function. In practice, as obtaining the full time‐dependent population size distribution is essential for the inference of growth parameters, the fluid limit has limited applicability. Nonetheless, for some growth models, it can provide explicit results for the average population size. Furthermore, by construction, the fluid limit provides an accurate description of the long‐time population growth dynamics but fails at describing the short‐time transient growth accurately. This is particularly limiting, for instance in the context of epidemics where the goal is to infer important parameters such as the reproduction rate as early as possible. On a more technical note, the random initial size distribution can be tedious to derive for some population growth models. In this context, we would be remiss not to mention a method based on a theorem by Cohn and Jagers (1994), which also allows one to describe the average stochastic population dynamics, but in this case, at any time (Uecker & Hermisson, 2011). However, as in the fluid limit, this method does not provide a way to obtain the population size distribution.

In contrast, we proposed two methods to derive exact solutions to the stochastic population growth either by solving the master equation directly, which requires the diagonalization of the transition rate matrix or by tracking reproduction times instead of population sizes. The first method was shown to be valid only in situations where the reproduction rates are distinct, whereas the second is generic. Our solution has revealed that the temporal distribution of population sizes is proportional to a hypoexponential distribution.

Note that as argued earlier, stochasticity in population growth can generically stem from intrinsic and extrinsic noises. Here, we focused on the effect of intrinsic noise and more specifically on the case where the only source of stochasticity resides in the waiting times between births, in an otherwise homogeneous population. A natural extension of our work would be to consider heterogeneous populations in which individuals display growth variability, that is, different intrinsic birth rates b. In this context, it would be interesting to study the effect of this intrapopulation variability on the growth dynamics both computationally and analytically by extending the present work.

It is also interesting to note that the role of extrinsic noise on population growth dynamics has also attracted attention in the past. A case in point is the study of the effect of stochastic temporal environmental fluctuations on growth (and other vital) rates and eventually population growth dynamics Steiner (2020). While early works focused on homogeneous populations (Lande et al., 2003), the works of Tuljapurkar and Orzack in the early 1980s extended these results to structured populations, with a particular focus on age‐structured matrix population models (Tuljapurkar & Orzack, 1980). In a series of seminal papers, they provided analytical insights into population structure and long‐run growth rates under stochastic environmental fluctuations showing that the distribution of population size is lognormal, and thus skewed (Tuljapurkar, 1982a). In particular, Tuljapurkar's small noise approximation, a perturbation formula used to calculate the population growth rate in stochastic environments, showed that population growth in constant and fluctuating environments differ, even when the mean growth rates are designed to be equal (Tuljapurkar, 1982b). Formally, it quantifies the decrease in the population growth rate in structured populations in response to the increased variance in vital rates due to randomly varying environments. In an interesting parallel with our own study, this result disproves the then commonly accepted assumption that averaging growth rates across realizations of the fluctuating environment would lead to an accurate estimate of the effective growth rates. Following on from Tuljapurkar's original work on sensitivity analyses to variance in vital rates, recent studies performed an age‐specific sensitivity analysis of stable, stochastic, and transient growth in structured populations (Giaimo & Traulsen, 2022).

Finally, we showed that our solution provides a more accurate description (than a deterministic approach) of the time‐dependent and steady‐state population sizes in a community composed of competing strains and the fixation probability of a mutation in a serial passage experiment. Thus, our theory offers an opportunity to quantify the dynamics of microbial communities from colonization to coexistence and thus contributes to the growing field of microbial eco‐evolutionary dynamics. For example, the gut of C. elegans worms colonized by two neutrally‐competing strains was shown to transition from a single‐strain composition at a low colonization rate to coexistence at a high colonization rate (Vega & Gore, 2017). In previous work, a deterministic approach was used to predict this transition (Vega & Gore, 2017). Our exact solution enables the expansion of this work by quantifying the abundance distribution of either strain within each worm's gut.

Our work opens new perspectives on population dynamics, ecology, and evolution. Importantly, our theory yields the full time‐dependent probability distribution of population sizes (see Equation (25)). Based on this result, an interesting future research direction would be to improve inference methods for growth parameters, as current methods suffer from substantial limitations (Ghenu et al., 2022). We propose that the present work constitutes the first step toward an exact inference method since it allows for the exact calculation of the likelihood function (Cranmer et al., 2020).

AUTHOR CONTRIBUTIONS

Loïc Marrec: Conceptualization (lead); data curation (lead); formal analysis (lead); investigation (lead); methodology (lead); project administration (lead); resources (lead); software (lead); validation (lead); visualization (lead); writing – original draft (equal); writing – review and editing (equal). Claudia Bank: Funding acquisition (lead); supervision (lead); validation (equal); writing – review and editing (equal). Thibault Bertrand: Formal analysis (supporting); investigation (supporting); methodology (supporting); resources (supporting); software (supporting); validation (supporting); visualization (supporting); writing – original draft (equal); writing – review and editing (equal).

FUNDING INFORMATION

CB is grateful for funding from ERC Starting Grant 804569829 (FIT2GO) and SNSF Project Grant 315230_204838/1 (MiCo4Sys).

CONFLICT OF INTEREST STATEMENT

We declare we have no competing interests.

Supporting information