key: cord-0010883-iwu022jh
authors: Spanio, Tommaso; Hidalgo, Jorge; Muñoz, Miguel A.
title: Impact of environmental colored noise in single-species population dynamics
date: 2017-10-02
journal: Phys Rev E
DOI: 10.1103/physreve.96.042301
sha: 8758b0e7492fdfd29fe3b4069fd42a29216abab4
doc_id: 10883
cord_uid: iwu022jh

Variability on external conditions has important consequences for the dynamics and the organization of biological systems. In many cases, the characteristic timescale of environmental changes as well as their correlations play a fundamental role in the way living systems adapt and respond to it. A proper mathematical approach to understand population dynamics, thus, requires approaches more refined than, e.g., simple white-noise approximations. To shed further light onto this problem, in this paper we propose a unifying framework based on different analytical and numerical tools available to deal with “colored” environmental noise. In particular, we employ a “unified colored noise approximation” to map the original problem into an effective one with white noise, and then we apply a standard path integral approach to gain analytical understanding. For the sake of specificity, we present our approach using as a guideline a variation of the contact process—which can also be seen as a birth-death process of the Malthus-Verhulst class—where the propagation or birth rate varies stochastically in time. Our approach allows us to tackle in a systematic manner some of the relevant questions concerning population dynamics under environmental variability, such as determining the stationary population density, establishing the conditions under which a population may become extinct, and estimating extinction times. We focus on the emerging phase diagram and its possible phase transitions, underlying how these are affected by the presence of environmental noise time-correlations.

Variability on external conditions has important consequences for the dynamics and the organization of biological systems. In many cases, the characteristic timescale of environmental changes as well as their correlations play a fundamental role in the way living systems adapt and respond to it. A proper mathematical approach to understand population dynamics, thus, requires approaches more refined than, e.g., simple white-noise approximations. To shed further light onto this problem, in this paper we propose a unifying framework based on different analytical and numerical tools available to deal with "colored" environmental noise. In particular, we employ a "unified colored noise approximation" to map the original problem into an effective one with white noise, and then we apply a standard path integral approach to gain analytical understanding. For the sake of specificity, we present our approach using as a guideline a variation of the contact process-which can also be seen as a birth-death process of the Malthus-Verhulst class-where the propagation or birth rate varies stochastically in time. Our approach allows us to tackle in a systematic manner some of the relevant questions concerning population dynamics under environmental variability, such as determining the stationary population density, establishing the conditions under which a population may become extinct, and estimating extinction times. We focus on the emerging phase diagram and its possible phase transitions, underlying how these are affected by the presence of environmental noise time-correlations. DOI: 10.1103/PhysRevE.96.042301

Population dynamics is a core matter in the modeling of ecological communities, genetics, and epidemics [1] [2] [3] . Combined with the increasing volume of available experimental (big) data, it constitutes a fundamental tool to shed light on the laws governing complex communities of living systems [4] . The traditional approach to population dynamics consists in the analysis of coupled deterministic equations describing the evolution of species abundances in a given community [5] . This procedure-whose outcome is not necessarily simple [6, 7] -is adequate in many cases. However, deterministic approaches neglect the effect of fluctuations, and these are now acknowledged to be both inherent and essential to the organization of communities of living systems [8] . On the one hand, the discreteness and finiteness of populations lead to demographic noise, which has been shown to be responsible for a wealth of nontrivial phenomena, such as the emergence of complex statistical patterns in neutral communities [9, 10] , quasiperiodic oscillations in prey-predator systems [11] , species formation [12] , and others [13] [14] [15] . On the other hand, populations are strongly affected by fluctuations in external conditions [16, 17] , which in most of the cases are highly unpredictable. This source of stochasticity, usually called environmental noise, can have important consequences for, e.g., ecosystem stability [18, 19] and evolutionary dynamics [20] [21] [22] [23] [24] , and fosters species coexistence [25] [26] [27] [28] .

Remarkably, theoretical and empirical evidence reveals that these phenomena strongly rely on a specific interplay between the characteristic timescale of environmental variations and the intrinsic timescale of the dynamics [27] [28] [29] [30] [31] . Owing to this, * hidalgo@pd.infn.it theoretical approaches have to be constructed beyond simple white-noise approximations, i.e., including "colored" (timecorrelated) noise [32, 33] . The question of how environmental colored noise affects population extinction has been widely studied in the literature, and there are contrasting positions on whether environmental fluctuations increase or decrease the risk of extinction, as this may actually depend on subtle differences of the underlying dynamics as well as the actual "color" of the fluctuations [34] [35] [36] [37] [38] [39] [40] . Of particular interest for our analyses here is the remarkable work of Kamenev et al. [40] , who analyzed a logistic growth population-dynamic model in which birth and death rates fluctuate in time, showing that, depending on the interplay between the system size and the temporal scale of the environment, the model exhibits qualitatively different functional dependencies of the mean extinction time with the system size.

In this paper, we bring the question of how time-correlated environmental noise affects population dynamics to the context of phase transitions and analyze in detail one of the most standard models in the study of population dynamics, the contact process [41, 42] -which can also be described as a birth-death process of the Malthus-Verhulst class-in the presence of time-correlated environmental variability [43] [44] [45] [46] [47] . To study this model we employ the so-called "unified colorednoise approximation," which is exact in the limits of very large and very short correlation times [48, 49] , and study the resulting effective (white-noise) problem employing a standard path integral approach; analytical results are tested against direct computational simulations obtained employing a very careful numerical analysis scheme [50] . Using this combined approach, we scrutinize the model phase diagram and identify the parameters for which the population becomes extinct with certainty and those for which the population survives, as well as the threshold separating them, and how the resulting phase transition depends on environmental-noise time correlations. As we will show, the phase diagram becomes much richer in this case than in its noiseless counterpart.

From a broader perspective, our study provides a simple and general framework for the analysis of population dynamics with colored noise, blending together several analytical [48, 49] and numerical [50] tools already available in the literature, which can be straightforwardly implemented to other similar scenarios beyond the case presented here.

The contact process (CP) [41, 42] is a prototypical model for the study of population dynamics with extinction, with applications in different fields such as in epidemic spreading [51] , ecology [52] , and propagation of neural activity [53] .

We use this simple model as a guideline here, but results are easily generalizable to other similar models for population or spreading dynamics. In the CP (see Fig. 1 ), nodes in a FIG. 1. Single-species dynamics with extinction. (a) A community of individuals grow under limited conditions with the dynamics of the contact process running on any given network. Each of the N nodes in the network can be either occupied by up to one individual (active node) or remain empty (inactive node). Individuals can reproduce at empty neighboring nodes at rate λ and also die and be removed from the community at rate μ. (b) The model exhibits different behaviors depending on parameter values: for low values of the reproduction rate, birth processes do not compensate deaths, and any population becomes extinct in the long term (absorbing phase). After a threshold value, population density exhibits a nonzero stationary value (active phase). Both regimes are separated by a critical point of a continuous phase transition. (c) However, any finite population eventually becomes extinct due to demographic fluctuations. Phases can be then distinguished looking at the scaling of the mean-extinction time with the system size, which is logarithmic in the absorbing phase and exponential in the active phase, while it becomes a power-law just at the critical point. (d) In the simplest scenario, the rate at which individuals reproduce (and similarly for the death rate) can be considered as a constant parameter (blue line). However, variability on the external conditions such as the temperature, humidity, pH,... may strongly influence rate parameters (red line). What is the impact of environmental variability in the previous panels?

given network (e.g., a square lattice) can be either occupied and active or vacant and inactive. Active nodes produce new offspring at neighboring empty sites at rate λ and can also die and be removed from the community at rate μ. The total system size N is fixed, representing limitation of space or resources, imposing an upper bound on the active population size. For the sake of simplicity, we neglect spatial effects and restrict our analysis to the simplest case of a well-mixed community (or, equivalently, a fully connected network). At each time t, the state of the system is determined by the total number of active sites, n(t), or equivalently, the population density, ρ(t) = n(t)/N. For very large populations, demographic fluctuations can be neglected, and the dynamics of ρ becomes deterministic [41, 42] :

The stationary density, ρ * [see Fig. 1 (b)] is either ρ * = 0 (the so-called "absorbing" state) if λ < λ c = μ, i.e., if births do not balance deaths and the population progressively shrinks, leading to extinction, or ρ * = 1 − μ/λ ("active" phase) if λ > λ c and the population survives indefinitely. A "critical" point, λ = λ c , separates the absorbing from the active phase; this value represents the extinction threshold and is a fundamental parameter of our forthcoming analysis.

In finite systems, demographic fluctuations can drive population extinction even when parameters correspond to the active phase [54] ; as a matter of fact, the only "truly" stable solution in the long term is the absorbing state in such a case. Still, it is possible to characterize the phases of a finite population by the mean time to reach extinction, T , as a function of the system size, N . Different functional dependencies emerge for each of the phases of the CP [43, 55] (see Fig. 1 ): T scales logarithmically with the system size in the absorbing phase (λ < λ c ), while it increases exponentially in the active phase (λ > λ c )-meaning that extinctions become extremely rare for sufficiently large populations in the active phase-and scales as a power-law right at the critical point, λ = λ c . Thus, in a nutshell, the CP represents a prototypical paradigm to analyze single-species communities with extinction. Such a dynamics can be characterized by means of (i) the phase diagram, which describes the stationary state of the system as a function of the parameters, (ii) the critical point, representing the extinction threshold, and (iii) the scaling of the mean-extinction time with the system size, as a proxy for population stability. In the next section we analyze how time-correlated environmental noise changes each of these elements and how this depends on the environmental auto-correlation time.

For the sake of simplicity, we assume that the environment influences homogeneously the population (i.e., demographic rates are global variables), and it does so by affecting only one parameter, which here we take to be the birth rate, leaving all other parameters unchanged (other choices are possible but they do not significantly affect the forthcoming results). The state of the environment is encoded in a time-dependent variable, (t), assumed to be independent of the state of system, so that λ → λ(t) =λ + σ (t), whereλ is the mean value, σ is a constant, and (t) follows an Ornstein-Uhlenbeck process [54] (see Fig. 1 

where ξ (t) is a zero-mean Gaussian white noise with

(2), it follows that λ is distributed as a Gaussian variable with meanλ and autocorrelation function (λ(t) −λ)(λ(t ) −λ) = σ 2 e |t−t |/τ [54] . As we are interested in the interplay between the timescale of the dynamics and the environment, we keep the correlation time of the environment, τ , as a control parameter throughout our analysis. Let us also note that, for sufficiently large values of σ , it may occur that λ(t) < 0, so we restrict our analysis to the regime of small variability, σ λ . Numerical analyses are performed keeping the constraint that λ(t) = 0 if a negative value is reached, but such events are extremely rare for σ λ . We have verified that other forms with bounded colored noise, as for instance dichotomous Markov noise [56] , do not change qualitatively our main results.

In the well-mixed scenario, substituting λ →λ + σ (t) in Eq. (1), one readily finds the following stochastic differential equation for the averaged density ρ in the infinite-size limit:

The set of stochastic equations formed by Eqs. (2) and (3) constitutes the starting point of our analysis. In this work, we exclusively focus on the impact of environmental fluctuations, as here we consider the limit of very large populations, where demographic fluctuations can be safely neglected. The important case with the combined effects of demographic and environmental stochasticity-relevant for finite systems-has also been explored in the literature [23, 27, 40, 57, 58] .

We analyze how environmental colored noise changes the phase diagram and mean-extinction times (as sketched in Fig. 1 ) by employing both analytical and computational approaches. To this end, we combine two different analytical tools: the unified colored noise approximation (UCNA) [48] (see also Appendix A) and a path-integral approach to calculate extinction times in finite populations [40] . Numerical simulations of the stochastic particle ("individual-based") model have been implemented by means of Anderson's next reaction method [50] , which can be adapted to the case in which rates vary stochastically in time (see Appendix B).

We first compute the stationary density as a function of parameter values. The process defined by Eqs. (2) and (3) is Markovian and, thus, the theory of Markovian processes applies [54] . The standard approach to solve it consists in finding the steady-state distribution P st (ρ, ) by solving the corresponding Fokker-Planck equation and then computing its associated averaged density ρ * = 1 0 dρρ ∞ −∞ d P st (ρ, ). However, this program cannot be completed analytically in the present case, as an exact integral does not exist.

The UCNA allows us to construct an approximate Markovian process for just one variable-much more susceptible of analytical understanding-describing the evolution of the population density with white noise [48] . In a nutshell, the UCNA method consists in the adiabatic elimination of the environmental variable; this can be safely done when the intrinsic dynamics and the environmental one operate at very different timescales. As a matter of fact, the method provides an exact equation for τ → 0 and τ → ∞, whereas it is only approximate for intermediate timescales τ . Thus, the UCNA can be understood as an "interpolation" between the dynamics for rapidly and slowly varying environments, respectively [59] (see Appendix A).

For simplicity, it is convenient to rewrite Eq. (3) in terms of a new variable with additive rather than multiplicative noise. In particular, defining

where, to ease the notation we have introduced the drift term f (x) =λ − μ − μe x . Unless explicitly stated, the forthcoming expressions remain valid for other choices of f (x). In the case of Eq. (4), after the elimination of the environmental variable , the UCNA approximation leads to a Langevin equation with multiplicative noise of the form (see Appendix A for a more detailed presentation of the UCNA method):

where η(t) is a Gaussian white noise with zero mean and η(t)η(t ) = δ(t − t ), which has to be understood in the Stratonovich sense [48] , and where

Equation (5) is equivalent to the following Fokker-Planck equation for the probability distribution P (x,t) [54] :

whose stationary solution, P st (x), can be found analytically imposing the zero-flux condition:

where we have introduced the potential V (x) = − x f (y)dy and Z is a normalization constant. Introducing the expression of f (x), i.e., Eq. (4), and reverting the change of variables, P st (ρ) = P st [x(ρ)]dx/dρ, one finally obtains the stationary distribution for the dynamics of Eq. (3) under the UCNA approximation:

Although Eq. (10) has a complicated form, its shape is chiefly controlled by the factor ρ¯λ −μ σ 2 τ −1 , as illustrated in Fig. 2 . In particular, ifλ < λ c = μ, then a nonintegrable singularity appears at ρ = 0, i.e., the distribution is not normalizable. This means that Eq. (10) is not a truly stationary solution, and the only stationary solution corresponds to the absorbing state, P st (ρ) = δ(ρ) [61] . On the other hand, forλ > λ c = μ, Eq. (10) can be safely normalized. Consequently, our first important result is that environmental variability does not shift the critical point for the process described by Eq. (3) .

Remarkably, there is an important difference with respect to the case with constant rates, as the active phase splits into two regions (see Fig. 2 ): the first one, spanning in λ c <λ < λ c + σ 2 τ ≡ λ c , for which P st (ρ = 0) exhibits an integrable singularity, and the second one, forλ > λ c , for which P (ρ = 0) = 0. Therefore, we find a regionλ ∈ [λ c ,λ c ], where the system can be found arbitrarily close to the absorbing state. We call this region the "weakly active" phase. Moreover, such a region can be itself divided into two subregions; after a certain valueλ = λ b the distribution becomes bimodal, only to recover its monomodality whenλ > λ c . The value of λ b does not have a simple analytical form and needs to be numerically determined.

We have compared the prediction of the PDF, Eq. (10), with simulations of the individual-based model. To capture the form of the quasistationary distribution (i.e., the distribution conditioned to the fact that the system is not in the absorbing state), it suffices to instantaneously introduce one active particle in the system if the population becomes extinct (more sophisticated methods provide slightly more accurate results of the quasistationary distribution [62] ). Figure 2(b) illustrates the time series of the population density for a system size N = 1000, for different values ofλ taken in the weakly active and active phases, respectively. Histograms are represented in Fig. 2(c) , together with the corresponding theoretical prediction, Eq. (10) (solid curves), illustrating a rather good agreement between them.

These results (summarized in Fig. 2 ) have been obtained neglecting the effects of demographic noise, which may play a fundamental role when the system approaches to the absorbing state. In fact, one may have doubts about the physical meaning of ρ * in the weakly active region, as the distribution exhibits a singularity at ρ = 0 and the system may become extinct with a relatively large probability. To shed light on this issue, in the next section we analyze the mean extinction time in the weakly active and active regions, elucidating the meaning of the different phases in the context of finite populations.

Mean-first passage times of a stochastic process [63] can be computed using the framework of path integrals [40, 49, 64] . Our strategy here is to apply such a method to the "effective" process obtained from the UCNA approximation method, Eq. (5).

The idea behind the path integral approach is that one can express the probability of a particular realization of the process, i.e., of a path {(x(t),ẋ(t))} t ≡ (x,ẋ), as

FIG. 2. Phases of the contact process with environmental noise (well-mixed scenario). Individual-based dynamics are implemented using Anderson's next reaction method [50] , where the death rate is fixed to μ = 1 and the birth rate is a stochastic Gaussian variable with meanλ, variance σ 2 , and temporal correlation τ [modeled as an Ornstein-Uhlenbeck process, Eq. (2)]. We study the behavior of the model for different values ofλ, and we take μ = 1, τ = 5, σ 2 = 0.1, N = 1000. Extinction is avoided by introducing an active particle at a random location to measure quasistationary distributions. (a) Phases of the model: (i) in the absorbing phase,λ <λ c = μ, the only stationary solution is extinction; (ii) in the weakly active phase,λ c <λ <λ c = μ + τ σ 2 , there is a positive quasi stationary density but the system can approach arbitrarily close to extinction due to fluctuations of the environment; this phase can be divided into two phases, depending on whether the corresponding probability distribution function (PDF) is uni-or bimodal; (iii) in the active phase, excursions close to the absorbing state become extremely rare. where S in the action associated with such a path, i.e., the time-integral of the Lagrangian, S[(x,ẋ)] = dtL(x(t),ẋ(t)), which encodes the dynamics we aim to describe (see below).

A particular realization of the process leading to extinction can be understood as a path passing through the state of zero density (the absorbing state), where the dynamics ceases. The leading contribution is given by the most probable path starting in the neighborhood of a deterministic attractor and ending at the absorbing state, i.e., the path (x * ,ẋ * ), which obeying such constraints minimizes the action. Up to leading order (in a weak-noise, i.e., low σ , expansion) the mean time to go to extinction is then inversely proportional to the probability of such a path [40] :

Let us now compute the action along the most probable path, following a standard procedure [40, 49, 64] . Given a stochastic process described by a Fokker-Planck equation [e.g., Eq. (7)], the time evolution can be described in terms of an associated Hamiltonian operator, ∂ t P = H P , which, as a rule of thumb, can be written by simply identifying −∂ x → p [40, 49] . In particular, for Eq. (7), (12) Given that the Hamiltonian does not depend explicitly on time (∂H/∂t = 0), it is a constant of motion. Moreover, as the optimal path we are looking for starts from the deterministic attractor (for which p = 0), such a constant is equal to zero [64] .

Imposing these constraints, one finds two solutions of H (x * ,p * ) = 0: the trivial one, p * = 0, corresponding to the deterministic trajectory toward the stable attractor, and another one for p * = 0, describing the most-probable fluctuation driving the system from an initial state to the absorbing state, which, for the case of Eq. (6), is (13) Given that the Lagrangian is the Legendre transform of the Hamiltonian, L(x,ẋ(x,p)) =ẋ(x,p)p − H (x,p) and that H = 0, the action can be easily evaluated as

without the need to explicitly integrate the equations of motion (Hamilton equations). Plugging Eq. (13) into this, and using Eq. (11) [65] , one readily obtains

where we have defined the potential V (x) = − x dx f (x ). This expression, which is valid for a general form of f (x), can be understood as a generalization of the Arrhenius formula with an effective diffusion term equal to τ σ 2 . Figure 3 (a) illustrates different trajectories in the (x,p) plane for the case of the CP with environmental noise under the UCNA approximation. Deterministic trajectories (p = 0; i.e., horizontal axis) push the system toward the stable   FIG. 3. (a) in a path-integral framework, a stochastic process follows trajectories in the (x,p) space. The deterministic dynamics (p = 0) lead the system to a stable attractor x * , from where the system can escape due to fluctuations (p = 0) and go to extinction, which here we identify with the state of one single particle remaining in the population, ρ = N −1 ⇒ x = − log(N − 1). Up to first order, the mean extinction time, T , is the exponential of the action associated to the most probable path, represented by the shaded region. deterministic attractor x * . On the other hand, a stochastic trajectory, i.e., (p = 0) starting arbitrarily close to the attractor, takes the system from there to the absorbing state that-to make it reachable in finite time-we identify with the state with one particle remaining in the system, ρ = 1/N [⇒ x = − log (N − 1) ]. The shaded area in Fig. 3 corresponds to the action S of the most probable stochastic trajectory.

More quantitatively: it is possible to derive the asymptotic behavior of Eq. (15) (i.e., its large-N behavior) for the specific case of Eq. (4). The initial point can be taken arbitrarily, as it does not depend on the system size, and the ending point scales as x f − log(N ) for large system sizes. With that, O(N −1 ) , which was introduced in Eq. (15) and leads to the final result,

for anyλ > μ. Therefore, our conclusion is that environmental noise induces power-law scaling of the mean-extinction time all along the active phase-typical for systems under environmental stochasticity [28, 43, 44, 46 2 N ) , which-neglecting other parameter constants-applies to our setting as we work in the asymptotic limit of large N).

Equation (16) elucidates the meaning of the phase diagram depicted in Fig. 2 : in the weakly active phase (μ <λ < μ + τ σ 2 ), the system makes excursions very close to the absorbing state, and as a consequence T scales sublinearly with the system size, whereas, in the active phase (λ > μ + τ σ 2 ) T scales superlinearly, and extinction becomes more unlikely for large system sizes. The linear case, in between, signals a change in the convexity of the extinction-time versus system-size curves.

We have checked the validity of Eq. (16) with an implementation of the individual-based model using Anderson's next step algorithm [50] (see Appendix B). To this end, we compute T through independent realizations setting as initial condition ρ = 1/2 and = 0, as a function of the system size N, for different values ofλ. As illustrated by Fig. 3(b) , Eq. (16) perfectly captures the asymptotic scaling behavior of T , whereas it fails for small values of N , where demographic stochastic effects (not included in the above calculation) significantly affect the dynamics.

We have presented a mathematical and computational study of a simple model for a well-mixed population where the dynamics is subjected to environmental variability, consisting of a contact process with birth rates modeled as an Ornstein-Uhlenbeck process. Our goal was to explore how the standard phase diagram of the contact process is affected by the introduction of environmental noise and in particular by its temporal autocorrelations. We explored whether its critical point is shifted or not, and what the nature of the emerging phases is. For this, we choose to work in the large system size limit, so that demographic fluctuations are negligible with respect to environmental ones, and the focus was put on phases and phase transitions.

The approach presented here is simple and easy to extend to other models and consists in the successive use of two analytical techniques: an approximation to deal with colored noise that reduces the number of variables, and a way to compute extinction times from the resulting equation. In particular, in the mean field-limit (describing a well-mixed scenario), we employ the unified colored noise approximation [48] (UCNA) to replace the correlated (colored) noise by a δ-correlated Gaussian (white) noise, at the price of introducing, an effective force and an effective diffusion term in the Langevin equation describing the system.

We verified computationally that the UCNA works remarkably well all across the phase diagram, generating steady-state density-distributions, hardly indistinguishable from the exact ones, as obtained from Monte Carlo simulations of the under-lying microscopic model. Numerical analyses were performed using a variation of the (exact) Gillespie integration method, adapted to deal with time-dependent (stochastic) rates, i.e., the so-called Anderson's method [50] . These analyses revealed that the probability distribution becomes a δ-Dirac at ρ = 0 at a critical value of the averaged birth rate coinciding with the critical value for the pure contact process. In other words, the introduction of colored environmental noise does not shift the location of the critical point, separating the absorbing from the active phase.

To proceed further we applied a weak-noise approximation within a path-integral formulation of the effective white-noise problem [40, 64] . Using this standard approach, a second important result is that in the active phase, the mean time required for the system to reach extinction scales as a powerlaw of the system size. This algebraic dependence of extinction times in the presence of uncorrelated environmental noise was first reported by Leigh [67] and was later on scrutinized by Vázquez et al. [43] who introduced the notion of "temporal Griffiths phase" (TGP) to refer to such a sort of active state. Temporal Griffiths phases are the counterpart of standard Griffiths phases [68] , but where the role of spatial quenched disorder is played by temporal one [43] [44] [45] [46] .

A more general study of extinction times in the presence of colored noise was elegantly tackled and solved by Kamenev et al [40] . These authors performed a path-integral formulation to the full problem, including both demographic and environmental stochasticity, and found diverse regimes depending on the ratio between system size and noise correlation time.

In particular, the dependence on system-size of the mean extinction time changes from exponential in the absence of the environmental noise (as corresponds to the Arrhenius law) to a power law for a short-correlated noise (as is the case in our study) and to no dependence whatsoever for noise with very large correlation times. This last regime implies that when there are extremely long periods of adverse external conditions (as compared, using adequate rescaling units, with system size) the system reaches deterministically the absorbing state regardless its size. Let us notice that this situation is not accessible to our approach, as we set the system in the asymptotic limit of large N (where demographic effects can be neglected), and thus τ cannot be much larger than it. For the power-law regime, our simple method gives exactly the same dependence as in Ref. [40] .

An interesting result of our analysis is that the active phase can be divided into two sub-phases. In the "weakly active phase" the probability distribution of ρ has a nonvanishing value around 0 meaning that the system makes excursions to very tiny values, at the edge of the collapse, but then it recovers. For this regime we find that, in finite systems, extinction times scale sublinearly with system size, while in the truly active phase, the scaling is superlinear. Note that the linear case in between signals a change of convexity in the extinction-time versus system-size curves in linear scale.

In summary, the combined use of the UCNA approximation to deal with colored noise and a standard weak-noise approximation for the resulting white noise equation allows us to construct a general approach to analyze particle systems with absorbing states in the presence of colored noise, in a relatively simple and systematic way. We believe that this combined approach should be very useful in applications in theoretical ecology, population dynamics, and epidemic spreading among others. We briefly review the method of the unified colored noise approximation (UCNA) [48] . We start from a general stochastic process with additive colored noise (see below for multiplicative noise):ẋ

where (t) is described by an Ornstein-Uhlenbeck process [Eq. (2)]. To confine the variable around a given bounded interval, we impose that f (x) < 0 for all values of x. (t) can be eliminated from Eq. (A1) by differentiating in time and introducing the expression for˙ , Eq. (2):

Multiplying both sides of this equation by τ and introducing a new time scalet = t/ √ τ , one obtains

where dots now refer to the time derivative in the scalet and where we have introduced the "damping" factor γ τ (x) defined as

Finally, the stochastic term in Eq. (A3) can be replaced by an equivalent one, η( √ τt) → τ −1/4η (t), where η(t)η(t ) = δ(t −t ), as both of them have the same mean value and correlation function, η(

Observe, from the definition of γ τ [Eq. (A4)], that the system becomes over-damped both in the limit τ → 0 and τ → ∞. Therefore, it is possible to perform an adiabatic approximation-valid for either very small or very large values of τ -by neglecting the transient contribution of the term x, and the process becomes approximately equivalent to the following Langevin equation with multiplicative noise:

Let us note that this equation has to be understood in the Stratonovich sense, as the previous derivation has been carried using the rules of the standard calculus [54] . Finally, we remark that the UCNA method becomes exact in the limit τ → 0 and τ → ∞, whereas it is an approximation for intermediate values of τ . Let us notice that, by Eq. (A4), for τ 0 the system is homogeneously over-damped (i.e., independently of x), whereas a dependency of x is preserved for τ 1. Therefore, we should expect that the UCNA provides accurate results independently of x for short-correlated environments, whereas in highly correlated environments we may still find some discrepancies with numerical results for those values of x for which f (x) is large.

The Gillespie algorithm is the most widespread method used to simulate the dynamics of discrete particlelike stochastic processes [69] . Alternatively, when rates explicitly depend on time, one can use the algorithm developed by Anderson [50] , which can be easily adapted to the case in which such a dependency is stochastic. In this appendix, we briefly review Anderson's next reaction method, underlining some technical issues arising from the fact that rates vary in time in a stochastic rather than deterministic way.

In short, Anderson's algorithm keeps track of different "clocks" counting the time remaining for each possible reaction to occur. In the "system of reference" of each reaction, its timer goes at unit speed; an absolute "observer" updates each of these timers according to their corresponding time-dependent rates, accelerating or slowing them, choosing which reaction will take place next and updating the species involved.

At each time, the state of the system is determined by the number of members or particles of each species, X. Each reaction k is characterized by its propensity a k (X,t) and its state-change vector ν k , so that X → X + ν k when reaction k occurs. The algorithm is implemented as follows [50] :

(1) Initialize: set the number of each species, set t = 0 and the internal clocks P k = T k = 0 for each reaction k.

(2) Generate the internal firing times for each k, which are exponentially distributed random variable with unit mean: P k = − log(r k ), where r k is a uniform random number in the interval [0,1].

(3) Calculate propensity functions a k (X,t) for each k.

(4) For each k, find the t k for which t+ t k t ds a k (X,s)ds = P k − T k (see below). (5) Find the minimum time-step = min k t k . We denote α its corresponding reaction.

(6) Update internal times T k = T k + t+ t a k (X,s) (see below), total time t = t + , and finally the number of species according to reaction α, X = X + ν α .

(7) Generate a new internal firing time P α = P α − log(r α ) with r α a uniform random number in (0,1), and go to step 3.

A small complication arises during steps 4 and 6 when rates vary stochastically, as the method has an anticipating nature: when integrating each of the propensity functions one may integrate one or many stochastic functions up to a certain time t k , and then seek for the minimum of such times, . In principle, this step does not require a complete knowledge of the whole stochastic trajectories at intermediate time steps (see below). However, after this, one should reintegrate all the propensity functions up to t + (step 6), but this can be a problem if we did not keep track of the stochastic trajectory with high enough precision.

We believe that this problem may be also solved using the theory of "stochastic bridges" [70] , which in principle would allow us to generate an intermediate point of the stochastic path conditioned to a future value of such a process (information that afterwards could be safely erased). We prefer to not enter into this matter, and we simply keep track of the stochastic trajectories (in our case the Ornstein-Uhlenbeck process and its time integral) taking a moving time window of length 10τ and precision 10 −2 τ , where τ is the correlation of the environment. Intermediate points of such a discretization are calculated using the simple linear interpolation rule.

When implementing Anderson's next reaction method, we use the formulas derived in Ref. [71] to integrate exactly the OU process [Eq. (2)]: 

where N 1 and N 2 are two independent Gaussian random numbers with zero mean and unit variance, and the coefficients θ = exp(− t/τ ), σ 2 1 = σ 2 (1 − θ 2 ), σ 2 2 = 2σ 2 τ 2 ( t/τ − 2(1 − θ ) + (1 − θ 2 )/2) and κ 12 = σ 2 τ (1 − θ 2 ) 2 . τ represents the temporal correlation of the process and σ its standard deviation. Notice that, when implementing Anderson's method, one could simply generate two random Gaussian numbers and look for the t for which t+ t t (s)ds = C, where C is some numerical value given by the algorithm. However, this leads to the already mentioned problem of anticipation, so it is preferable to integrate Eqs. (B1) and (B2) at regular time steps.

Theoretical Ecology: Principles and Applications

The Neutral Theory of Molecular Evolution

Infectious Diseases of Humans: Dynamics and Control

Biophysics: Searching for Principles

Mathematical Biology

Stability and Complexity in Model Ecosystems

Noise-induced Phenomena in the Environmental Sciences

Proc. Natl. Acad. Sci. USA

Proc. Natl. Acad. Sci. USA

Evolution in Changing Environments: Some Theoretical Explorations

Proc. Natl. Acad. Sci. USA

Nonequilibrium Phase Transitions in Lattice Models

Stochastic Interacting Systems: Contact, Voter and Exclusion Processes

Spatial Ecology: The Role of Space in Population Dynamics and Interspecific Interactions

Non-equilibrium Phase Transitions: Absorbing Phase Transitions, Theoretical and Mathematical Physics

As the noise has a finite correlation, the standard rules of calculus apply here

A Guide to First-passage Processes

Field Theory of Non-equilibrium Systems

In the UCNA approximation we have rescaled the time variable, but this change only introduces a √ τ factor in the mean extinction times