Population Extinction on a Random Fitness Seascape

The growth of spatially distributed populations is of importance to ecology, evolutionary biology, epidemiology, among many other fields. One of the simplest models of a species or pathogen propagating in space is the Fisher equation Fisher (1937); Aronson and Weinberger (1978); Gourley (2000); Kwapisz (2000); Hallatschek and Korolev (2009) given by

$$\frac{dy(x,t)}{dt}=\mu y-ay^{2}+D\nabla^{2}y~{}.$$

(1)

At each spatial coordinate $x$, the population size $y(x,t)$ initially increases exponentially in time $t$, at rate determined by the fitness parameter $\mu>0$. This (logistic type) growth is eventually slowed down by the nonlinear term, describing resource limitations, with the population asymptoting to $\mu/a$ (the capacity). Alternatively, for $\mu<0$, the population decays to zero. The diffusion term $D$ captures migration of the population to nearby points in space.

Such deterministic evolution is only valid for populations of infinite size for which fluctuations can be ignored. In particular, near extinction points where the population hits zero (such as when a species suffers ecological collapse or when a pathogen is without available hosts), the deterministic description does not capture the stochasticity of small number fluctuations. Transitioning from a model with a finite, discrete population into a continuum formulation results in a noise of amplitude proportional to the square root of the number of individuals. Including this so-called demographic noise Durrett and Levin (1994); Butler and Goldenfeld (2009); Traulsen et al. (2012); Constable et al. (2016); Weissmann et al. (2018) leads to the stochastic Fisher equation

$$dy=\left(\mu y-ay^{2}+D\nabla^{2}y\right)dt+\sigma_{d}\sqrt{y}dW~{},$$

(2)

with $\sigma_{d}$ as the noise amplitude, and $dW=dW(x,t)$ being a standard Wiener process. Models of this type have been studied in a number of contexts outside of population growth, thanks to connections to directed percolation Cardy and Sugar (1980); Täuber et al. (1998); Janssen and Täuber (2005); Peschanski (2009); Sipos and Goldenfeld (2011); Korolev et al. (2010); Horowitz and Kardar (2019) (with onset of laminar turbulence for $\mu=a=0$ Sipos and Goldenfeld (2011) as a recent application). Indeed, extinction transitions are expected to generically belong to the directed percolation universality class. Henkel et al. (2008); Taüber (2014)

The main focus of this paper is on a different form of noise. The implicit assumption of the Fisher equation is that the fitness $\mu$ is an intrinsic property, uniform in time and the same across all locations. For most populations, however, the growth rate is strongly influenced by external factors, such nutrient level, light, etc. at different times and places. Generalizing from a static but position dependent fitness landscape $\mu(x)$, a time-varying fitness $\mu(x,t)$ is referred to as a seascape Mustonen and Lässig (2009, 2010); Iram et al. (2019); Agarwala and Fisher (2019). Population growth on a random seascape is thus governed by

$$dy=\left(\mu y-ay^{2}+D\nabla^{2}y\right)dt+\sigma ydW,$$

(3)

with $\sigma^{2}$ being the variance of the fitness noise, and $dW=dW(x,t)$ a Wiener process as before. Note that for $a=0$, the linear equation describes the evolving weight $y(x,t)$ of a directed polymer in a quenched random energy landscape ($\sigma dW(x,t)$); with $\ln y(x,t)$ satisfying the Kardar-Parisi-Zhang equation, notable in the study of surface growth problems Kardar et al. (1986); Kardar and Zhang (1987); Derrida (1990); Chu and Kardar (2016); Borodin et al. (2014); Medina et al. (1989); Quastel and Spohn (2015); Sasamoto and Spohn (2010); Halpin-Healy (2013).

While much studied, the problem of directed polymers in random media is largely unsolved in other than one dimension. Such limitation will necessarily extend to characterizing the effect of a random seascape on a population in two or more dimensions. We thus resort to replacing near-neighbor migrations with jumps of arbitrary length. This “mean-field” limit is reasonable for a well mixed bacterial population, or for a pathogen infecting an active city. An alternative description is to replace continuous space with a series of nodes indexed by $k=1,2,\cdots,N$. If population $y_{k}$ of each node can jump to any other node in the network at rate $D/N$, then (for $N\gg 1$), evolution in a random seascape is described by

$$dy_{k}=\{\mu y_{k}-a{y_{k}}^{2}+D(\bar{y}-y_{k})\}dt+\sigma y_{k}dW_{k}~{},$$

(4)

and similarly for demographic noise we have

$$dy_{k}=\{\mu y_{k}-a{y_{k}}^{2}+D(\bar{y}-y_{k})\}dt+\sigma_{d}\sqrt{y_{k}}dW_{k}~{},$$

(5)

where $\bar{y}=\frac{1}{N}\sum_{k}y_{k}$ is the spatial average, with $dW_{k}=dW(k,t)$ being spatiotemporal noise. In a sense, the diffusion coefficient $D$ now behaves like a mean-field coupling strength. The problem can be solved exactly in this case, as for $N\to\infty$, the probability distribution for $\bar{y}$ converges to a $\delta$-function, reducing the above equations to those of independent variables whose distributions can be obtained by standard methods. The important population average $\bar{y}$ is then computed self-consistently.

We carry our this program below initially for seascape noise, in the three cases of: neural evolution ($\mu=a=0$) in Sec. II, decaying population ($\mu<0$) in Sec. III, and saturating growth ($\mu>0$, $a>0$) in Sec. IV. For comparison, the cases of demographic noise, as well as mixed demographic and seascape noise, are studied Sec. V.

Despite the commonality in their deterministic parts, the two types of stochastic noise relate to two quite distinct universality classes; that of directed polymers in random media for Eq. (4), versus directed percolation for Eq. (5). As such, we might expect them to exhibit different behaviors in the all-important extinction regime. This is indeed the case: Near the extinction threshold, we show that seascape noise leads to a broadly distributed population, characterized by a tail falling off as a power-law. The power-law exponent varies continuously with $\sigma^{2}/D$ such that, when the noise is large, fluctuations in the population much exceed the mean. The extinction transition no longer belongs to the (“mean-field”) directed percolation universality class. However, as we show in Sec. V, demographic noise restores the directed percolation universality class.

Many subtleties of seascape noise can be gleamed by considering Eq. (4) in the ‘neutral’ case, given by

$$dy_{k}=D\left(\bar{y}-y_{k}\right)dt+\sigma y_{k}dW_{k}~{}.$$

(6)

The above equation can be interpreted as describing a form of synchronization: each $y_{k}$ is being pulled towards the spatial mean $\bar{y}$ with strength $D$, but they are forced apart by the noise $\sigma y_{k}dW_{k}$. The interplay between the two tendencies manifests in the population variance, with one limiting form being a “synchronized” population described by a tight distribution, and the other being an “incoherent” population that is broadly distributed Strogatz (2003); Pikovsky et al. (2003).

While it would be interesting to analyze the case of a finite $N$, computations simplify in the thermodynamic limit of $N\to\infty$. If starting with a uniform initial condition, say $y_{k}(t=0)\equiv y_{0}$ for all $k$, then the evolving probability distributions of all individual members will be identical. Thus, if there exists some long-time stationary distribution with finite mean, then this mean will be the same for each member. So by simply applying the law of large numbers and taking expectations, we obtain $\bar{y}=\langle y\rangle=y_{0}$. The second equality follows from the conservation of average population in the neutral limit, remaining at $Ny_{0}$ at all times.

Therefore, the large-$N$ dynamics are captured by the single stochastic differential equation

$$dy=D(y_{0}-y)dt+\sigma ydW~{},$$

(7)

independent of the spatial index $k$. From here, we can construct a Fokker-Plank equation for the evolution of the probability density function $\rho(y,t)$, namely

$$\partial_{t}\rho=-\partial_{y}\left[D(y_{0}-y)\rho-\frac{\sigma^{2}}{2}{\partial_{y}}\left(y^{2}\rho\right)\right]~{}.$$

(8)

A stationary distribution ($\partial_{t}\rho\equiv 0$) is obtained by setting the probability current in the above square bracket to zero, leading to

$$0=y^{2}\rho^{\prime}+(c_{D}+2)y\rho-c_{D}y_{0}\rho~{},$$

(9)

where we define $c_{D}=2D/\sigma^{2}$. This is a solvable ordinary differential equation, with a properly normalized solution given by

$$\rho(y)=\frac{\left(c_{D}y_{0}\right)^{c_{D}+1}}{\Gamma(c_{D}+1)}y^{-2-c_{D}}e^{-c_{D}y_{0}/y}~{},$$

(10)

which is a scaled inverse Chi-squared distribution with $\nu=2(c_{D}+1)$ and $\tau=y_{0}c_{D}/(1+c_{D})$. This prediction agrees with simulated infinite-$N$ dynamics depicted in Fig. 1.

An important characteristic of this solution is the behavior in the tail, where for large populations $y$, the density behaves as

$$\rho\propto y^{-2-c_{D}}~{};$$

(11)

confirmed in Fig. 2. Even though the mean of the distribution stays at the initial value of $y_{0}$, the seascape fluctuations lead to a power-law tail which results in a variance (and higher moments) which may either be finite or infinite! Specifically, we have

$$\frac{\langle y^{2}\rangle}{\langle y\rangle^{2}}=\frac{2D}{2D-\sigma^{2}}~{},$$

(12)

such that the variance becomes infinite when $c_{D}$ crosses 1, and the population distribution becomes broad compared to the mean. Parenthetically, we note that such a “transition” does not occur in case of demographic noise as described later in the text.

Now consider Eq. (4) in the presence of negative fitness $\mu$. The nonlinear term is now asymptotically irrelevant as can be seen by the transformation to $z=ye^{-\mu t}$, such that the dynamics in this ‘moving frame’ becomes

$$dz=\left[-az^{2}e^{\mu t}+D\left(\bar{z}-z\right)\right]dt+\sigma zdW~{}.$$

(13)

Since $\mu<0$, the nonlinear term drops out at large times, leaving us with

$$dz=D\left(\bar{z}-z\right)dt+\sigma zdW~{},$$

(14)

which is simply the exact same dynamics as for $\mu=0$ in Eq. (6). Therefore, the same long-time results hold in these transformed coordinates. This includes the exact distribution predicted in Eq. (10) as well as the moment ratio in Eq. (12), as confirmed by the simulations depicted in Fig. 3. Therefore, while the mean population size decays exponentially, the ratio

$$\frac{\langle y^{2}\rangle}{\langle y\rangle^{2}}=\frac{2D}{2D-\sigma^{2}},$$

(15)

indicates that the mean ceases to be a good measure of the distribution for $2D/\sigma^{2}<1$.

The most intricate case corresponds to when $\mu$, $a$, $D$, and $\sigma$ are all strictly positive:

$$dy_{k}=\{\mu y_{k}-a{y_{k}}^{2}+D\left(\bar{y}-y_{k}\right)\}dt+\sigma y_{k}dW_{k}~{}.$$

(16)

The intricacy arises as the noise modifies the mean value of $y$ from the bare value of $\mu/a$. In reality, the dynamics of the mean depend on the value of every higher moment. More intuitively, note that the restoring force of the logistic growth term is asymmetric, in that it “punishes” overly large values of $y$ harder than small ones. Under the effect of noise, this biases $y$ to smaller values, and we should expect the true value of $\langle y\rangle$ to be smaller than $\mu/a$.

To properly analyze this case, we need to construct a self-consistent solution in the long-time, large–$N$ limit. Assuming a long-time stationary distribution exists, and given that every node has the same initial condition $y_{k}(t=0)=y_{0}$, then the uniformity of their stochastic dynamics indicates that all nodes arrive to the same distribution $\rho$, for all $k$, i.e. $\bar{y}=\int y\rho(y)dy$. The stationary distribution $\rho(y)$ is obtained by examining the Fokker-Plank equation associated with the stochastic differential equation

$$dy=[\mu y-a{y}^{2}+D\left(\bar{y}-y\right)]dt+\sigma ydW~{},$$

(17)

given by

$$\partial_{t}\rho=-\partial_{y}\left[\left(\mu y-a{y}^{2}+D(\bar{y}-y)\right)\rho-\frac{\sigma^{2}}{2}{\partial_{y}}\left(y^{2}\rho\right)\right]~{}.$$

(18)

Setting the probability current to zero leads to the ordinary differential equation

$$0=y^{2}\rho^{\prime}+(2+c_{D}-c_{\mu})y\rho-c_{D}y_{0}\rho+c_{a}\rho y^{2}~{},$$

(19)

where $c_{D}=2D/\sigma^{2}$, $c_{\mu}=2\mu/\sigma^{2}$, and $c_{a}=2a/\sigma^{2}$. This admits an unnormalized solution of the form

$$\hat{\rho}(y)=e^{-c_{D}\bar{y}/y}e^{-c_{a}y}y^{-2-c_{D}+c_{\mu}}~{}.$$

(20)

Note that finite $\bar{y}$ and $a$ cut off power-law tails of the distribution for small and large $y$, respectively, so that all moments of the distribution are now finite.

To convert Eq. (20) to a proper PDF, we need the normalization factor

	$\displaystyle Z$	$\displaystyle=\int_{0}^{\infty}\hat{\rho}(y)dy$
		$\displaystyle=2\left(\frac{c_{a}}{c_{D}\bar{y}}\right)^{\frac{1+c_{D}-c_{\mu}}{2}}\mathcal{K}_{1+c_{D}-c_{\mu}}\left(2\sqrt{c_{D}c_{a}\bar{y}}\right)~{},$		(21)

where $\mathcal{K}_{\gamma}$ is a modified Bessel function of the second kind. The mean of the distribution is then given by

	$\displaystyle\langle y\rangle Z$	$\displaystyle=\int_{0}^{\infty}y\hat{\rho}(y)dy$
		$\displaystyle=2\left(\frac{c_{a}}{c_{D}\bar{y}}\right)^{\frac{c_{D}-c_{\mu}}{2}}\mathcal{K}_{c_{D}-c_{\mu}}\left(2\sqrt{c_{D}c_{a}\bar{y}}\right)~{},$		(22)

leading to the self-consistency equation,

$$\bar{y}=\frac{c_{D}}{c_{a}}\left(\frac{\mathcal{K}_{\beta}(2x)}{\mathcal{K}_{\beta+1}(2x)}\right)^{2}~{},$$

(23)

where we define $\beta=c_{D}-c_{\mu}$ and $x=\sqrt{c_{D}c_{a}\bar{y}}$ for the sake of convenience. The above equation can be solved numerically; as depicted in Fig. 4, there is excellent agreement between the thus analytically constructed PDF and the result from of numerical simulation.

We may well ask if the extinction transition described above persists with the addition of demographic noise, i.e. for the equation

$$dy_{k}=\mu y_{k}-a{y_{k}}^{2}+D(\bar{y}-y_{k})+\sigma y_{k}dW_{k}+\sigma_{d}\sqrt{y_{k}}dW^{\prime}_{k}~{}.$$

(31)

We first consider demographic noise by itself (with $\sigma=0$), and then the mixed case ($\sigma_{d}\neq 0$ and $\sigma\neq 0$).

Both seascape fluctuations and demographic noise have well defined origins. The former arises from natural variations in population fitness and access to resources, whereas the latter is the result of moving from a discrete to a continuum description. Both are unavoidable in real world systems, and neglecting one or the other can cause drastically incomplete predictions of the breadth of a population. In particular, seascape fluctuations broaden the tail of the distribution at large values, while demographic stochasticity controls its behavior close to zero.

To understand the importance of seascape fluctuations, it is instructive to consider Eq. (3) in the limit of $a=0$ and $D=0$. At each point $\ln y$ now performs a random walk, resulting in broad log-normal distributions for which characterizing the population in terms of mean $\langle y\rangle$ and variance $\langle y^{2}\rangle_{c}$ is inadequate. A finite diffusion coefficient (still with $a=0$) is expected to modify the log-normal distribution. The resulting behavior, in the universality class of directed polymers in random media (DPRM), is highly dependent on dimensionality: In one dimension it leads to the celebrated Tracy–Widom distribution Dotsenko (2010); Quastel and Spohn (2015); Chu and Kardar (2016), while in dimensions larger than 2, a critical value of $D_{c}$ separates broad and narrow distributions Halpin-Healy (2012); Kardar and Zhang (1987); Halpin-Healy (2013). To circumvent difficulties associated with spatial dimensionality, in this paper we consider a mean-field limit (Eq. (4)) in which the parameter $D$ accounts for migration between any pair of sites, irrespective of their separation. In this case, seascape noise of strength $\sigma$ generates a power-law tail in the distribution for large populations which falls off with exponent $2+2D/\sigma^{2}$. Much like the case of DPRM at high dimensions, this signals a transition from broad to a narrow distributions for $D>D_{c}=\sigma^{2}/2$.

For $a=0$, a finite $\mu$ simply makes the mean, $\bar{y}$, time dependent (proportional to $e^{\mu t}$) without modifying the shape of the distribution. For $\mu>0$, a finite $a$ is needed to curtail the exponential growth of the population; a self-consistency requirement now enables computing a stationary $\bar{y}$. The finite $a$ also leads to an exponential decay of the distribution at large values, resulting in well-behaved moments. The most interesting aspect of the problem is then the nature of the extinction transition as $\mu\to 0$: For $D>D_{c}$, the mean (and all other moments) vanish linearly with $\mu$, as expected for a mean-field extinction transition. However, for $D<D_{c}$ the mean vanishes as $\mu^{\sigma^{2}/(2D)}$, with higher moments vanishing with other exponents (harking back to the broad distribution expected for $\mu=0$).

While the above unconventional critical behavior is an accurate description of Eq. (4) as $\mu\to 0$, its relevance to an actual species extinction must also account for demographic stochasticity. For $D=0$, demographic noise in Eq. (5) is incompatible with a continuous probability distribution, with fluctuations causing local extinctions resulting in a growing delta-function in the PDF at the absorbing point $y=0$. A finite $D$ removes the delta-function by seeding new population at all sites through migration from other sites, resulting in the steady state distribution of Eq. (33), with a normalizable divergence at the origin. (For $|\mu|\ll D$, the same distribution, albeit with a time dependent $\bar{y}(t)$ can be used to describe growing or shrinking populations.) However, the boundary between growth and extinction is no longer at $\mu=0$ in the presence of demographic noise, and the self-consistency condition should be used to identify a critical value of $\mu_{c}$. The nature of the extinction transition at $\mu_{c}$ is now conventional, belonging to the expected directed percolation universality class. Finally, we find that when both seascape and demographic noise are present, the resulting extinction transition (at a finite $\mu_{c}$) is again conventional. This is not too surprising as extinction characterizes behavior when $y\to 0$, where demographic fluctuations are paramount.

The overarching theme of this investigation is the emphasis on the need to characterize the full distribution of a population and not just its mean. In the course of this investigation, in the mean-field limit, we have encountered a few distributions of some repute: For seascape noise, a scaled inverse chi-squared distribution appears; for demographic noise an Erlang shows up, and when both noises are mixed (for $\mu=a=0$), a rescaled beta prime distribution manifests. The special functions characterizing these distributions are typically solutions of relatively simple differential equations. In our context, these equations appear naturally in computing steady states of Fokker-Planck equations, providing examples of their applicability.

The covid-19 epidemic has motivated myriad current studies of spread and control of epidemics Giordano et al. (2020); Tuite et al. (2020); Chatterjee et al. (2020); Kucharski et al. (2020). With notable exceptions (e.g. Ref. Bittihn and Golestanian (2020)), most studies follow a number of time-dependent quantities that are implicitly assumed to characterize the whole distribution. The approach of this paper can be used to examine actual distributions via the self-consistency condition in the mean-field limit, for multiple stochastic equations, generalizing Eq. (V.4). Even for a single quantity, we may generalize the exponent of the stabilizing nonlinear term from 2 to an arbitrary $\alpha$. This is not without motivation, since instead of an embedded logistic equation at every point in space, we would now have a Richards’ equation. The Richards’ equation acts as a generalization of logistic growth, and has been used in a variety of ecological models Richards (1959); Fekedulegn et al. (1999). Indeed, most results presented in this paper have natural generalizations for $\alpha>1$, including corresponding distributions. The most notable complexity arises in the self-consistency equation for the mean, which requires numerical analysis and lacks closed forms.

Another interesting query is the applicability of the mean-field results to the finite dimensional extinction problem. There is a certain richness of dimensional dependence in directed percolation models, but we may nonetheless expect some qualitative behavior in the mean-field case to carry over to finite dimensions Korolev et al. (2010). For example, in the important cases of one and two dimensions, the limiting distribution for $a=0$ arising from Eq. (3) is always broad. It would be interesting to investigate the nature of the extinction transition with seascape noise (with and without demographic noise) in these dimensions.

Finally, we note that certain predictions in this model may be experimentally replicable. Well-mixed microbial populations are often represented via a mean-field migration model, which we employ here. While tracking entire distributions is cumbersome (albeit not impossible), the mean and standard deviation of the population should be relatively easy to track. Our model predicts that the ratio of these quantities makes a dramatic change, passing from finite to very large as the strength of seascape noise is increased above a threshold value. The fitness of a population can be controlled via the resources afforded to the microbes, and local fluctuations in this quantity can be constructed by randomizing access to resources (e.g. light levels). While the results presented in this paper may not be robust to various perturbations, an experiment probing the effects of seascape fluctuations in fitness should yield highly informative results.

This research was supported by a McDonnel fellowship to Bertrand Ottino-Löffler, as well as by NSF through grants # DMR-1708280 and # PHY-2026995 (MK).