Extinction phase transitions in a model of ecological and evolutionary dynamics

TL;DR

This study investigates how environmental fluctuations influence the extinction transition in ecological models, revealing a novel infinite-noise universality class driven by temporal disorder.

Contribution

It demonstrates that temporal disorder causes an exotic infinite-noise phase transition, contrasting with the standard directed percolation class in static environments.

Findings

Extinction transition in static environments belongs to directed percolation universality class.

Temporal disorder induces a logarithmically slow decay and large fluctuations in populations.

Evidence supports the transition being of infinite-noise type, with temporal Griffiths phases.

Abstract

We study the non-equilibrium phase transition between survival and extinction of spatially extended biological populations using an agent-based model. We especially focus on the effects of global temporal fluctuations of the environmental conditions, i.e., temporal disorder. Using large-scale Monte-Carlo simulations of up to $3\times 10^7$ organisms and $10^5$ generations, we find the extinction transition in time-independent environments to be in the well-known directed percolation universality class. In contrast, temporal disorder leads to a highly unusual extinction transition characterized by logarithmically slow population decay and enormous fluctuations even for large populations. The simulations provide strong evidence for this transition to be of exotic infinite-noise type, as recently predicted by a renormalization group theory. The transition is accompanied by temporal Griffiths phases featuring a power-law dependence of the life time on the population size.

Tables

Table 1. Table 1: Critical exponents for the bacterial splitting model in a time-independent environment compared to high-accuracy results for the two-dimensional directed percolation (DP) universality class. The exponents above the horizontal line are directly measured in our simulations, the ones below are computed using exponent relations.

Exponent	this work	DP, Ref. Dickman99	DP, Ref. VojtaFarquharMast09
$\delta$	0.454(2)	0.4523(10)	0.4526(7)
$\Theta$	0.235(10)	0.2293(4)	0.233(6)
$z$	1.76(1)	1.767(1)	1.757(8)
$z{\nu }_{⟂}$	1.29(2)	1.292(4)	1.290(4)
${\nu }_{⟂}$	0.73(2)
$\beta$	0.586(12)

Equations

$$W(p)=c\delta (p-p_h)+(1-c)\delta (p-p_l)$$

$$\rho (\Delta ,t)=b^{-\beta /{\nu }_{\perp }}\rho (\Delta b^{1/{\nu }_{\perp }},t{b}^{-z}).$$

$$\rho (\Delta ,t)=t^{-\delta }{X}_{\rho }(\Delta t^{1/(z{\nu }_{\perp })})$$

$$P_s(t)\sim t^{-\delta },N_s(t)\sim t^{\Theta },R(t)\sim t^{1/z}.$$

$$\rho (t)\sim [\ln (t/t_0){]}^{-\bar{\delta }}$$

$$\rho (\Delta ,t)=(\ln b{)}^{-\beta /{\bar{\nu }}_{\perp }}\rho [\Delta (\ln b{)}^{1/{\bar{\nu }}_{\perp }},t{b}^{-z}],$$

$$\rho (\Delta ,t)=[\ln (t/t_0){]}^{-\beta /{\bar{\nu }}_{\perp }}{X}_{\rho }\{\Delta [\ln (t/t_0){]}^{1/{\bar{\nu }}_{\perp }}\}$$

$$P_s(t)$$

$$N_s(t)$$

$$R(t)$$

Full text

1 Introduction

The scientific study of population growth and extinction has a long history. It dates back, at least, to the efforts of Euler and Bernoulli in the 18th century Euler1748 ; Bernoulli1760 . Early work was based on simple deterministic equations in which any spatial dependence was averaged out. Later, stochastic versions of these models were also considered. Recently, attention has turned to approaches that provide more realistic descriptions by incorporating fluctuations in space and time as well as features such as heterogeneity and mobility (see, e.g., Refs. Bartlett61 ; Britton10 ; OvaskainenMeerson10 and references therein).

The behavior of a biological population close to the transition between survival and extinction is of great conceptual and practical importance. On the one hand, extinction transitions are prime examples of non-equilibrium phase transitions between active (fluctuating) and inactive (absorbing) states, a topic of considerable current interest in statistical physics MarroDickman99 ; Hinrichsen00 ; Odor04 ; HenkelHinrichsenLuebeck_book08 . On the other hand, efforts to predict and perhaps even prevent the extinction of biological species on earth require a thorough understanding of the underlying mechanisms. The same also holds for the opposite problem, viz., efforts to predict, control and eradicate epidemics.

(Random) temporal fluctuations of the environment play a crucial role for the behavior of a biological population close to extinction. In contrast to intrinsic demographic noise, environmental noise causes strong population fluctuations even for large populations which makes extinction easier. As a result, the mean time to extinction for uncorrelated environmental noise only grows as a power of the population size rather than exponentially as it would for time-independent environments Leigh81 ; Lande93 . Recent activities have also analyzed the effects of noise correlations in time on these results (for a recent review, see, e.g., Ref. OvaskainenMeerson10 ).

Most of the above work focused on space-independent (single-variable or mean-field) models of population dynamics. Significantly less is known about the effects of temporal environmental fluctuations on the dynamics of spatially extended populations and, in particular, on their extinction transition. Kinzel Kinzel85 established a general criterion for the stability of a non-equilibrium phase transition against such temporal disorder. Jensen Jensen96 ; Jensen05 studied directed bond percolation with temporal disorder and reported that the extinction transition is characterized nonuniversal critical exponents that change continuously with disorder strength. Vazquez et al. VBLM11 demonstrated that the power-law dependence between population size and mean time to extinction also holds for spatially extended systems in some parameter region around the extinction transition (which they called the temporal Griffiths phase).

More recently, Vojta and Hoyos VojtaHoyos15 used renormalization group arguments to predict a highly unconventional scenario, dubbed infinite-noise criticality, for the extinction transition in the presence of temporal environmental fluctuations. This prediction was confirmed for the one-dimensional and two-dimensional contact process by extensive Monte-Carlo simulations BarghathiVojtaHoyos16 .

Here, we study the transition between survival and extinction of a spatially extended biological population by performing large-scale Monte-Carlo simulations of an agent-based off-lattice model. In the absence of temporal environmental fluctuations, we observe a conventional extinction transition in the well-known directed percolation universality class of non-equilibrium statistical physics. Environmental fluctuations are found to destabilize the directed percolation behavior. The resulting extinction transition is characterized by logarithmically slow population decay and enormous fluctuations even for large populations. Numerical data for the time-dependence of the average population size, the spreading of the population from a single organism, as well as the mean time to extinction are well-described within the theory of infinite-noise critical behavior VojtaHoyos15 .

Our paper is organized as follows. In Sec. 2, we introduce our model and the Monte-Carlo simulations. Section 3 reports on the extinction transition in time-independent environments, while section 4 is devoted to the case of temporally fluctuating environments. In the concluding Sec. 5, we discuss how general these findings are, we compare them to earlier simulations, and we discuss generalizations to correlated disorder as well as spatial inhomogeneities.

2 Model and simulations

2.1 Definition of the model

The model was originally conceived DeesBahar10 ; SKMB13 as a model of evolutionary dynamics in phenotype space but it can be interpreted as a model for population dynamics in real space as well. Organisms reside in a continuous two-dimensional space of size $L\times L$. In the evolutionary interpretation, the two coordinates correspond to two independent phenotype characteristics (in arbitrary units) while they simply represent the real space position of the organism for population dynamics.

The time evolution of the population from generation to generation consists in three steps, (i) reproduction, (ii) competition death, and (iii) random death. In the reproduction step, each organism produces $N_{fit}$ offspring. We consider two different reproduction schemes, as sketched in Fig. 1.

In the asexual fission (or “bacterial splitting”) scheme, each organism produces offspring independently of the other organisms. The offspring coordinates are chosen at random from a square of side $2\mu$ centered at the parent position. In the evolutionary context, $\mu$ represents the mutability. In the assortative mating scheme, each organism mates with its nearest neighbor. More specifically, if organisms $a$ and $b$ are the nearest neighbors of each other, the pair produce $2\,N_{fit}$ offspring. If, on the other hand, organism $b$ is the nearest neighbor of $a$, but $c$ is the nearest neighbor of $b$, organism $b$ will produce $N_{fit}$ offspring with $a$ and another $N_{fit}$ offspring with $c$. The offspring coordinates are randomly chosen from a rectangle around both parents, extended by $\mu$ in each direction.

After the reproduction step, the parent organisms are removed (die). The offspring then undergo two consecutive death processes. First, if the distance between two offspring organisms is below the competition radius $\kappa$, one of them is removed at random. This process simulates, e.g., the competition for limited resources. Second, each surviving offspring dies with death probability $p$. These random deaths are statistically independent of each other and model predation, accidents, diseases, etc. Offspring that survive both death processes form the parent population for the next generation. Note that the maximum number of organisms that can survive without competition-induced deaths, $N_{max}$, corresponds to a hexagonal close packing of circles of radius $\kappa/2$. Its value reads $N_{max}=(2/\sqrt{3})(L/\kappa)^{2}$.

If the parameters of this model, i.e., the number $N_{fit}$ of offspring per parent as well as mutability $\mu$, competition radius $\kappa$, and death probability $p$ do not depend on the position in space, no organism is given a fitness advantage. This corresponds to neutral selection in the evolutionary context. In the present paper, we only consider this case, but we will briefly discuss the effects of a nontrivial fitness landscape in the concluding section.

While we assume the model parameters to be uniform in space, we do investigate variations in time (temporal disorder). They model temporal fluctuations in the population’s environment, e.g., climate fluctuations. These fluctuations are global because they affect all organisms in a given generation in the same way. Specifically, we contrast the case of a constant, time-independent death probability $p$ and the case in which $p$ varies randomly from generation to generation.

2.2 Monte Carlo simulations

Close to the extinction transition, the population is expected to fluctuate strongly on large length and time scales. This implies that large systems and long simulation times are necessary to obtain accurate results. We have therefore simulated large systems with sizes of up to $1500\times 1500$. For a competition radius of $\kappa=0.25$, this allows populations of up to 42 million organisms to exist without competition-caused death. These systems are at least three orders of magnitude larger than those used in earlier simulations DeesBahar10 ; SKMB13 ; SKOB17 .

To tune the population between survival and extinction, we vary the mutability $\mu$ while all other parameters are held fixed, including the competition radius $\kappa=0.25$ and the number of offspring per parent, $N_{fit}=2$. In simulations with time-independent death probability, we fix its value at $p=0.3$. Temporal disorder is introduced by making the value of $p$ for each generation an independent random variable drawn from a binary distribution

[TABLE]

with $p_{h}=0.5$, $p_{l}=0.1$, and $c=0.5$.

We perform two different types of simulations. The first kind starts from large populations of about 80% of the competition-limited maximum $N_{max}$ defined above (33 million organisms for the largest systems). These organisms are initially placed at random positions, and we follow their time evolution for up to 100,000 generations. To obtain good statistics, we average over up to 40,000 individual runs for each parameter set. The second kind of simulation starts from a single organism and observes the growth of the population and their spreading over space. Here, we average over more than $10^{6}$ individual runs.

3 Extinction transition in time-independent environments

3.1 Asexual fission (bacterial splitting)

In order to find the extinction transition, we analyze the time dependence of the number $N$ of organisms, starting from large populations (80% of $N_{max}$) for many different values of the mutability $\mu$, as illustrated in Fig. 2.

This figure shows a double-logarithmic plot of population density $\rho=N/N_{max}$ vs. time $t$ (measured in generations) for parameters $L=1500$, $\kappa=0.25$, $N_{fit}=2$, and $p=0.3$,

Three different regimes are clearly visible. For $\mu$ larger than some critical value $\mu_{c}$, the population density initially decays but then settles onto a constant, time-independent value. This is the active, surviving phase of the model. For $\mu<\mu_{c}$, the population density decays to zero faster than a power law; replotting the data in semilogarithmic form confirms an exponential decay. This is the inactive phase in which even large populations go extinct quickly. The two phases are separated by the critical point $\mu=\mu_{c}$ of the extinction transition for which the population density decays like a power law.

To determine the value of $\mu_{c}$, we fit the population density to the power-law form $\rho(t)\sim t^{-\delta}$. For $\mu=0.32810$, we obtain a high-quality fit (reduced $\chi^{2}\approx 1.0$) over the time interval from $t=200$ to the longest times. The data for $\mu=0.32809$ and 0.32811 yield significantly worse fits. We thus conclude $\mu_{c}=0.32810(1)$ where the number in brackets indicates the error of the last digit. We have compared the results of different landscape sizes to ensure that this value is not affected by finite-size effects. The power-law fit of the critical curve ($\mu=0.32810$) yields an exponent $\delta=0.4537$ with a very small statistical error of about $10^{-4}$. A larger contribution to the error stems from the remaining uncertainty in $\mu_{c}$. By comparing the fits for $\mu=0.32809$ and $\mu=0.32811$, we estimate this error to be 0.002. Our final result for the decay exponent therefore reads $\delta=0.454(2)$.

To analyze the off-critical behavior of the population density (values of $\mu$ close to but different from $\mu_{c}$), we employ a scaling ansatz appropriate for a nonequilibrium continuous phase transition (see, e.g., Ref. Hinrichsen00 )

[TABLE]

Here, $\Delta=\mu-\mu_{c}$ denotes the distance from the transition point, and $b$ is an arbitrary length scale factor. $\beta$, $\nu_{\perp}$, and $z$ are the order parameter, correlation length, and dynamical critical exponents, respectively. Setting the arbitrary scale factor to $b=t^{1/z}$ yields the scaling form

[TABLE]

of the density, with $\delta=\beta/(z\nu_{\perp})$ and scaling function $X_{\rho}$. (At criticality, $\Delta=0$, eq. (3) turns into the power-law decay $\rho\sim t^{-\delta}$ discussed earlier.) The scaling form (3) can be used to find the exponent combination $z\nu_{\perp}$: If one plots $\rho\,t^{\delta}$ vs. $\Delta t^{1/(z\nu_{\perp})}$, the data for all values of $\mu$ and $t$ should collapse onto a single master curve. To perform this analysis, we fix the decay exponent at the value $\delta=0.454$ found earlier and vary $z\nu_{\perp}$ until the best data collapse is obtained. This yields $z\nu_{\perp}=1.29$. The collapse for $z\nu_{\perp}=1.27$ or 1.31 is of visibly lower quality. We therefore arrive at the final estimate $z\nu_{\perp}=1.29(2)$. The resulting scaling plot is shown in Fig. 3.

In addition to the simulations that follow the time evolution of large populations, we also perform runs that start from a single organism in the center of the landscape and observe the spreading of the population through space. Specifically, we analyze the survival probability $P_{s}(t)$, i.e., the probability that the population has not gone extinct at time $t$. We also measure the average number $N_{s}(t)$ of organism in the population at time $t$ as well as the mean-square radius $R(t)$ of the cloud of organisms in space. Right at the extinction transition, these quantities are expected to follow the power laws

[TABLE]

From the quality of the fits of our data to these power laws, we determine the critical point to be $\mu_{c}=0.32809(2)$, in good agreement with the estimate obtained from the density $\rho(t)$ above. The resulting exponent values read $\delta=0.449(7)$, $\Theta=0.235(10)$, and $z=1.76(1)$. The value of $\delta$ is slightly less precise but agrees with the value obtained from $\rho(t)$. Figure 4 shows $P_{s}(t)$, $N_{s}(t)$, and $R(t)$ at the critical point.

The values of the critical exponents resulting from our simulations are summarized in Table 1.

They fulfill the hyperscaling relation GrassbergerdelaTorre79 $\Theta+2\delta=d/z$ within their error bars. For comparison, Table 1 also shows the exponents obtained in two high-accuracy studies of the contact process Dickman99 ; VojtaFarquharMast09 which is a paradigmatic model in the directed percolation GrassbergerdelaTorre79 universality class. Our exponents agree very well with the directed percolation values, providing strong evidence for the extinction transition in our model to belong to this universality class. This is in agreement with a conjecture by Janssen and Grassberger Janssen81 ; Grassberger82 , according to which all absorbing state transitions with a scalar order parameter, short-range interactions, and no extra symmetries or conservation laws belong to the directed percolation universality class.

3.2 Assortative mating

In addition to the bacterial splitting case, we also perform simulations using the assortative mating reproduction scheme for which the offspring is located in a rectangle around the organism and its nearest neighbor, expanded by the mutability $\mu$ (see Fig. 1). It is important to point out that assortative mating introduces a kind of long-range interaction between the organisms because the distance between an organism and its nearest neighbor can become arbitrarily large. It is thus not clear whether or not one should expect the extinction transition to belong to the directed percolation universality class.

Figure 5 shows the time evolution of the population density $\rho$ starting from a large population (80% of $N_{max}$) for many different values of $\mu$. The parameters $L=1500$, $\kappa=0.25$, $N_{fit}=2$, and $p=0.3$ are chosen to be identical to the bacterial splitting simulations discussed in the last section.

To determine the critical point, we fit the population density to the power law $\rho(t)\sim t^{-\delta}$. The best power law is found for $\mu=0.39736$ (high-quality fit with reduced $\chi^{2}\approx 1$ from $t=500$ to the longest times) while the fits for $\mu=0.39734$ and $0.39738$ are of lower quality. We therefore conclude that the critical point is at $\mu_{c}=0.39736(2)$. The exponent resulting from the power-law fit of the critical curve is $\delta=0.459(5)$. The majority of the error stems from the uncertainty in $\mu_{c}$; the statistical error is only about $2\times 10^{-4}$.

To study the off-critical behavior, we again perform a scaling analysis of the population density according to eq. (3). We fix $\delta$ at the value found above, $\delta=0.459$, and vary $z\nu_{\perp}$ until we find the best data collapse. This yields $z\nu_{\perp}=1.31(2)$. The resulting high-quality scaling plot is shown in Fig. 6.

The critical exponents of the extinction transition for the assortative mating reproduction scheme, $\delta=0.459(5)$ and $z\nu_{\perp}=1.31(2)$ are very close to the directed percolation exponents listed in Table 1. Specifically, $z\nu_{\perp}$ agrees within the error bars with the directed percolation value while $\delta$ is very slightly too high (the errors bars almost touch, though). We believe the small deviation stems from finite-size effects which are significantly stronger than in the bacterial splitting case. Simulations for a smaller landscape of $L=500$ yield an exponent value of $\delta=0.466$ which suggests that the true (infinite system) exponent will be slightly below 0.459(5).

We conclude that the extinction transition in the assortative mating case also belongs to the directed percolation universality class despite the potential long-range interactions introduced by the assortative mating. This may be caused by the fact that almost all mating pairs have short distances. A long mating distance occurs only if a single organism is far way from all others.111In contrast, for anomalous directed percolation with long-range spreading, every organism can produce offspring far away from its own position JOWH99 ; HinrichsenHoward99 . In this case, the positions of its offspring (and offspring of offspring) move quickly towards those of the other organisms. The net effect is not that different from the isolated organism simply dying.

Alternatively, the directed percolation exponents could hold in a transient time regime only while the long-range interactions cause the true asymptotic exponents to differ from directed percolation. However, our simulations of large populations for fairly long times do not show any indications of a crossover away from directed percolation behavior. Determining rigorously whether or not the assortative mating model belongs to the directed percolation universality class therefore remains a task for the future.

4 Extinction transition fluctuating environments

We now turn to the effects of temporal environmental fluctuations, i.e., temporal disorder, on the extinction transition. To this end, we perform simulations (using the bacterial splitting reproduction scheme) in which the death probability $p$ varies randomly from generation to generation. Specifically, the value of $p$ for each generation is an independent random variable drawn from the binary distribution (1) with $p_{h}=0.5$, $p_{l}=0.1$, and $c=0.5$.

The environmental noise induces enormous population fluctuations even for large populations. This is illustrated in Fig. 7 which shows the time evolutions of three individual populations subject to different realizations of the temporal disorder.

All three populations initially have the same large population of $N_{0}=0.8\,N_{max}\approx 1.5\times 10^{7}$ organisms. Nonetheless, after fewer than 100 generations, the populations already differ by more than two orders of magnitude. For comparison, the figure also shows three populations not subject to environmental noise but only to the demographic (Monte Carlo) noise of the stochastic time evolution. In this case, the populations curves are smooth and differ from each other by only a few percent. They are thus indistinguishable in Fig. 7. Environmental noise is much more efficient in inducing population fluctuations than intrinsic demographic noise because it affects the entire population in the same way while the demographic noise acting on one organism is independent from that acting on other organisms.

To investigate the extinction transition quantitatively, we now analyze the time evolution of the average population density. Figure 8 presents a double logarithmic plot of $\rho$ vs. $t$, averaged over 20,000 to 40,000 temporal disorder realizations.

As in the clean case, Fig. 2, we can identify the active, surviving phase in which the average population density $\rho$ approaches a nonzero constant for long times as well as the inactive phase, in which the density decays faster than a power law with time. However, the curve separating the two phases is clearly not a power law. In fact, none of the curves can be described by a power law over any appreciable time interval. This suggests that the extinction transition in the presence of temporal disorder, i.e., environmental noise, is unconventional.

Recently, Vojta and Hoyos VojtaHoyos15 developed a real-time (“strong-noise”) renormalization group approach to absorbing state transitions with temporal disorder. This theory predicts an exotic “infinite-noise” critical point. It is characterized by a slow, logarithmic decay of the average population density at the transition point,

[TABLE]

where $\bar{\delta}=1$, and $t_{0}$ is a non-universal time scale. To test this prediction, we plot in Fig. 9 the density data in the form $\rho^{-1}$ vs. $\ln(t)$. In this plot, the predicted critical behavior (5) yields a straight line, independent of the value of $t_{0}$.

The figure shows that the data for $\mu=0.4129$ indeed follow a straight line for more than two orders of magnitude in time. Accordingly, a fit of these data to eq. (5) from $t=50$ to the longest times ($t=10,000$) is of high quality (reduced $\chi^{2}\approx 1$), significantly better than the fits for $\mu=0.4127$ and $\mu=0.4130$. (The resulting $t_{0}$ value is $t_{0}=2.572$.) Consequently, we conclude that the critical point is located at $\mu=0.4129(2)$ and fulfills the predicted infinite-noise behavior.

Based on the renormalization group approach VojtaHoyos15 , a heuristic scaling theory was developed in Ref. BarghathiVojtaHoyos16 . The scaling ansatz for the average population density reads

[TABLE]

and the predicted exponent values are $\beta=0.5$, $\bar{\nu}_{\perp}=0.5$, and $z=1$. Setting the length scale factor $b=(t/t_{0})^{1/z}$ leads to the scaling form

[TABLE]

where the non-universal time scale $t_{0}$ is necessary because logarithms are not scale free. Interestingly, due to the logarithmic time dependence, the value of the dynamical exponent $z$ does not appear in (7). For $\Delta=0$, we recover the logarithmic decay (5) with $\bar{\delta}=\beta/\bar{\nu}_{\perp}=1$.

The scaling form implies that the data for all $\mu$ and $t$ should collapse onto a master curve if one graphs $\rho\,\ln(t/t_{0})$ vs. $\Delta\,[\ln(t/t_{0})]^{2}$. The corresponding scaling plot, using the value $t_{0}=2.572$ found earlier is presented in Fig. 10.

Clearly, the data collapse is not particularly good. A better collapse can be achieved if we allow the microscopic time scale on the $x$-axis to differ from $t_{0}$. Within the scaling description of the critical point this corresponds to a subleading correction to the leading scaling behavior. Figure 11 shows a scaling plot of $\rho\,\ln(t/t_{0})$ vs. $\Delta\,[\ln(t/t_{1})]^{2}$ where $t_{1}$ is independent of $t_{0}$.

The value $t_{1}=13.0$ leads to nearly perfect data collapse. Alternatively, a good collapse can be obtained by varying the exponent in the scaling combination $\Delta\ln(t/t_{0})]^{1/\bar{\nu}_{\perp}}$ on the $x$-axis. Figure 12 demonstrates that $1/\bar{\nu}_{\perp}=2.8$ results in a good collapse.

Distinguishing between these two scenarios, viz., (i) critical exponents as predicted by the renormalization group theory, but with corrections to scaling and (ii) an exponent value that differs from the renormalization group theory, requires data over significantly longer times that are presently beyond our numerical means.

In addition to the simulations studying the time evolution of large populations, we also perform runs that start from a single organism and observe the spreading of the population through space. The scaling theory of Ref. BarghathiVojtaHoyos16 . leads to the following predictions for the time dependencies of the survival probability $P_{s}$, the number of organisms $N_{s}$, and the radius $R$ of the organism cloud at criticality:

[TABLE]

Interestingly, except for the subleading logarithmic corrections, these time dependencies are the same as in the active, surviving phase where the population spreads ballistically, $R\sim t$ and $N_{s}\sim t^{2}$. The theory does not predict the values of the exponents $y_{N}$ and $y_{R}$ governing the logarithmic corrections.

Figure 13 presents the results of our spreading simulations at the critical point, $\mu=0.4129$. The data are averages over 150,000 realizations of the temporal disorder with 10 attempts of growing the population per realization. Each attempt starts from a single organism in the center of a landscape of size $L=1500$. The other parameters are identical to those used earlier, $\kappa=0.25$, $N_{fit}=2$, and $p=0.1$ or 0.5 with 50% probability.

The figure shows that the data for times larger than about 100 are well described by the theoretical predictions. Accordingly, fits to eqs. (8), (9), and (10) yield reduced $\chi^{2}\approx 1$. The resulting values of the exponents $y_{N}$ and $y_{R}$ have small statistical errors of about $3\times 10^{-2}$. However, as they stem from subleading terms, they are sensitive towards small changes of the fit interval. Taking this into account, we arrive at the estimates $y_{N}=3.9(4)$ and $y_{R}=1.1(2)$. These exponents fulfill the equality BarghathiVojtaHoyos16 $y_{N}=2\bar{\delta}+dy_{R}$ with $\bar{\delta}=1$. However, they differ from the values found in simulations of the two-dimensional contact process BarghathiVojtaHoyos16 . This can either mean that these subleading exponents are non-universal or that our values have not yet reached the asymptotic long-time regime. Distinguishing these possibilities will require simulations over much longer time intervals.

Finally, we determine how the mean time to extinction (or average life time) of a population depends on its size. In time-independent environments, a finite population that is nominally on the active, surviving side of the extinction transition can decay only via a rare fluctuation of the demographic (Monte-Carlo) noise. The probability of such a rare event decreases exponentially with population size, leading to an exponential dependence of the lifetime on the population size.

In the presence of (uncorrelated) temporal disorder, the extinction probability is enhanced due to rare, unfavorable fluctuations of the environment. Within space-independent (mean-field) models of population dynamics, this leads to a power-law relation (rather than an exponential one) between life time and population size Leigh81 ; Lande93 . Vazquez et al. VBLM11 pointed out that this power-law behavior can be understood as the temporal analog of the Griffiths singularities known in spatially disordered systems. They dubbed the parameter region where the power-law behavior occurs the temporal Griffiths phase.

In order to address this question in our model, we determine the survival probability $P(t)$ (the probability that the population has not gone extinct at time $t$) of finite populations close to the extinction transition. Figure 14 shows $P(t)$ for $\mu=0.420$ and several population sizes.

Fitting the survival probability to the exponential $P(t)\sim\exp(-t/\tau)$ yields the life time $\tau$ for each population size. We perform analogous simulations for a range of mutabilities from $\mu=0.4$ on the inactive, extinct side of the extinction transition to $\mu=0.44$ on the active, surviving side. The resulting graph of life time $\tau$ vs. system size $L$ (or, equivalently, initial population $N_{0}$) is presented in Fig. 15.

For mutabilities $\mu>\mu_{c}=0.4129$ (populations on the active side of the extinction transition), the figure indeed shows a power-law dependence $\tau\sim N_{0}^{1/\kappa}$ of the life time $\tau$ on the initial population $N_{0}$. The exponent $\kappa$ is non-universal; power-law fits of our data give $\kappa=0.53$, 0.66, and 0.89 for mutabilities $\mu=0.440$, 0.430, and 0.420, respectively. The real-time renormalization group of Ref. VojtaHoyos15 predicts that $\kappa$ increases as $\mu$ approaches the extinction transition and (in two space dimensions) reaches the value $\kappa_{c}=2$ right at $\mu_{c}$. We also expect a logarithmic correction to the leading power law analogous to those in eqs. (8), (9), and (10). We therefore fit the lifetime curve for $\mu=\mu_{c}=0.4129$ to the function $\tau\sim N_{0}^{1/2}[\ln(N_{0})]^{y_{\tau}}$ which leads to a high-quality fit with $y_{\tau}\approx 2.5$. On the inactive side of the extinction transition, the theory predicts a slow logarithmic increase, $\tau\sim\ln(N_{0})$, of the life time with population size. This behavior is indeed found for the $\mu=0.400$ curve in Fig.  15.

5 Conclusion

In summary, we have studied the extinction transition of a spatially extended biological population by performing Monte Carlo simulations of an agent-based model. For time-independent environments, i.e., in the absence of temporal disorder, we have found the extinction transition to belong to the well-know directed percolation universality class for both the asexual fission (bacterial splitting) and the assortative mating reproduction schemes. In the bacterial splitting case, this behavior is expected from the conjecture by Janssen and Grassberger Janssen81 ; Grassberger82 , according to which all absorbing state transitions with a scalar order parameter, short-range interactions, and no extra symmetries or conservation laws belong to the directed percolation universality class. The assortative mating case is more interesting because each organism mates with its nearest neighbor which can, in principle, be arbitrarily far away. Even though this introduces a long-range interaction in the model, our simulations still yield directed percolation critical behavior.

The main part of our paper has been devoted to the effects of temporal environmental fluctuations, i.e., temporal disorder, on populations close to the extinction transition. The question of whether or not a given universality class is stable against (weak) temporal disorder is addressed by Kinzel’s generalization Kinzel85 of the Harris criterion (see also Ref. VojtaDickman16 for version of the criterion that applies to arbitrary spatio-temporal disorder). According to the criterion, a critical point is stable against temporal disorder if the correlation time exponent $\nu_{\parallel}=z\nu_{\perp}$ fulfills the inequality $z\nu_{\perp}>2$. The directed percolation universality class in two dimensions features a value $z\nu_{\perp}\approx 1.29$ (see Table 1) which violates Kinzel’s criterion, predicting that temporal disorder must qualitatively change the extinction transition.

Simulations of our model in which the death probability $p$ varies randomly from generation to generation confirm this expectation. They yield an unconventional extinction transition characterized by logarithmically slow population decay and enormous fluctuations even for large populations. The time-dependence of the average population size, the spreading of the population from a single organism, and the mean time to extinction are well-described by the theory of infinite-noise critical behavior VojtaHoyos15 . This is also important from a conceptual statistical mechanics point of view because it confirms the universality of the infinite-noise scenario by showing that it holds for off-lattice problems as well as for simple lattice models such as the contact process with temporal disorder BarghathiVojtaHoyos16 .

It is worth pointing out that the description of the off-critical behavior required us to include a correction-to-scaling term. This likely stems from the fact that the logarithmically slow dynamics leads to a slow crossover to the true asymptotic regime. This slow crossover may also explain why recent simulations of a similar model SKOB17 appear to be compatible with directed percolation behavior even though that model does contain temporal disorder. For weak temporal disorder, the population is expected to show directed percolation behavior in a transient time regime before the true asymptotic critical behavior is reached. As the systems in Ref. SKOB17 are much smaller than ours (up to 30,000 organisms vs. $3.3\times 10^{7}$), their simulations probably do not reach the asymptotic regime.

The slow crossover to the asymptotic regime also suggests that the power-law scaling with nonuniversal exponents observed by Jensen in a directed percolation model with temporal disorder Jensen96 ; Jensen05 does not represent the true asymptotic behavior. Instead, the nonuniversal power laws hold in a transient time interval before the crossover to the logarithmic infinite-noise behavior.

The temporal disorder considered in the present paper is uncorrelated, i.e., “white noise”. The effects of long-range noise correlations (leading to “red noise”) on population extinction have attracted considerable interest in recent years (see, e.g., Ref. OvaskainenMeerson10 and references therein). Our simulations can be easily generalized to this case. From the analogy with long-range correlated spatial disorder IbrahimBarghathiVojta14 , we expect that long-range correlations of the temporal disorder further increase its effects and lead to a change of the universality class.

In the present paper, we have investigated the extinction transition of populations in a homogeneous landscape. In the evolutionary context of the model this corresponds to neutral selection. What about the effects of spatial inhomogeneities? According to the Harris criterion Harris74 $d\nu_{\perp}>2$ (where $d$ is the space dimensionality), (uncorrelated) spatial disorder is a relevant perturbation and destabilizes the directed percolation universality class. Renormalization group calculations HooyberghsIgloiVanderzande03 ; HooyberghsIgloiVanderzande04 and extensive Monte Carlo simulations of the contact process DickisonVojta05 ; VojtaFarquharMast09 ; Vojta12 have established that the resulting critical point is of exotic infinite-randomness kind. We expect similar behavior for spatially disordered versions of the present model.

If the system contains (uncorrelated) disorder in both space and time, the generalized Harris criterion governing the stability of the pure critical behavior reads $(d+z)\nu>2$ VojtaDickman16 ; AlonsoMunoz01 . The directed percolation universality fulfills this criterion in all dimensions suggesting that such spatio-temporal disorder is an irrelevant perturbation, at least if it is sufficiently weak (see also discussion in Ref. Hinrichsen00 ).

This work was supported by the NSF under Grant Nos. DMR-1205803 and DMR-1506152. We acknowledge valuable discussions with S. Bahar.