key: cord-0456416-l6bl1ui6
authors: Ottino-Loffler, Bertrand; Kardar, Mehran
title: Population Extinction on a Random Fitness Seascape
date: 2020-07-23
journal: nan
DOI: nan
sha: 4b471c1be44e866917a70a3999fed16d2dbaf2af
doc_id: 456416
cord_uid: l6bl1ui6

Models of population growth and extinction are an increasingly popular subject of study. However, consequences of stochasticity and noise in shaping distributions and outcomes are not sufficiently explored. Here we consider a distributed population with logistic growth at each location, subject to"seascape"noise, wherein the population's fitness randomly varies with {it location and time}. Despite its simplicity, the model actually incorporates variants of directed percolation, and directed polymers in random media, within a mean-field perspective. Probability distributions of the population can be computed self-consistently; and the extinction transition is shown to exhibit novel critical behavior with exponents dependent on the ratio of the strengths of migration and noise amplitudes. The results are compared and contrasted with the more conventional choice of demographic noise due to stochastic changes in reproduction.

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 [1] [2] [3] [4] [5] given by dy(x, t) dt = µy − ay 2 + D∇ 2 y .

At each spatial coordinate x, the population size y(x, t) initially increases exponentially in time t, at rate determined by the fitness parameter µ > 0. This (logistic type) growth is eventually slowed down by the nonlinear term, describing resource limitations, with the population asymptoting to µ/a (the capacity). Alternatively, for µ < 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 [6] [7] [8] [9] [10] leads to the stochastic Fisher equation

with σ 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 [11] [12] [13] [14] [15] [16] [17] (with onset of laminar turbulence for µ = a = 0 [15] as a recent application). Indeed, extinction transitions are expected to generically belong to the directed percolation universality class. [18, 19] The main focus of this paper is on a different form of noise. The implicit assumption of the Fisher equation is that the fitness µ 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 µ(x), a time-varying fitness µ(x, t) is referred to as a seascape [20] [21] [22] [23] . Population growth on a random seascape is thus governed by dy = µy − ay 2 + D∇ 2 y dt + σydW,

with σ 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 (σdW (x, t)); with ln y(x, t) satisfying the Kardar-Parisi-Zhang equation, notable in the study of surface growth problems [24] [25] [26] [27] [28] [29] [30] [31] [32] . 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, · · · , N . If population y k of each node can jump to any other node in the network at rate D/N , then (for N 1), evolution in a random seascape is described by dy k = {µy k − ay k 2 + D(ȳ − y k )}dt + σy k dW k , (4) and similarly for demographic noise we have

whereȳ = 1 N 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 → ∞, the probability distribution forȳ converges to a δ-function, reducing the above equations to those of independent variables whose distributions can be obtained by standard methods. The important population averageȳ is then computed self-consistently.

We carry our this program below initially for seascape noise, in the three cases of: neural evolution (µ = a = 0) in Sec. II, decaying population (µ < 0) in Sec. III, and saturating growth (µ > 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 allimportant 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 σ 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

The above equation can be interpreted as describing a form of synchronization: each y k is being pulled towards the spatial meanȳ with strength D, but they are forced apart by the noise σ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 [33, 34] . While it would be interesting to analyze the case of a finite N , computations simplify in the thermodynamic limit of N → ∞. If starting with a uniform initial condition, say y k (t = 0) ≡ 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 y = y = y 0 . The second equality follows from the conservation of average population in the neutral limit, remaining at N y 0 at all times.

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

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

A stationary distribution (∂ t ρ ≡ 0) is obtained by setting the probability current in the above square

Simulated versus predicted power-law decays of the tail of the probability distribution in the steady-state of Eq. (6) . Circles show numerically obtained decay exponents from simulating 5 × 10 6 replications of Eq. (7); the solid line depicts the predicted decay exponent at large y from Eq. (11) ; and the dashed line shows the maximum error arising from measuring the exponent at finite y, as predicted via Eq. (10). The simulation was run using stochastic Runge-Kutta, with y 0 = 1, D = 0.5, σ = 1.0, and µ = a = 0. The time-step of integration was 0.01, with a total integration time of 25 units. Regression to identify the decay rate of the PDFs began at y * ≈ 13.34, and the final exponent was averaged over the final 5 time units with samples every 0.05 time units.

bracket to zero, leading to

where we define c D = 2D/σ 2 . This is a solvable ordinary differential equation, with a properly normalized solution given by

which is a scaled inverse Chi-squared distribution with ν = 2(c D + 1) and τ = 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

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

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 µ. The nonlinear term is now asymptotically irrelevant as can be seen by the transformation to z = ye −µt , such that the dynamics in this 'moving frame' becomes dz = −az 2 e µt + D (z − z) dt + σzdW .

Since µ < 0, the nonlinear term drops out at large times, leaving us with

which is simply the exact same dynamics as for µ = 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

indicates that the mean ceases to be a good measure of the distribution for 2D/σ 2 < 1.

The most intricate case corresponds to when µ, a, D, and σ are all strictly positive:

The intricacy arises as the noise modifies the mean value of y from the bare value of µ/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 y to be smaller than µ/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 ρ, for all k, i.e.ȳ = yρ(y)dy. The stationary distribution ρ(y) is obtained by examining the Fokker-Plank equation associated with the stochastic differential equation

given by

(18) Setting the probability current to zero leads to the ordinary differential equation

where c D = 2D/σ 2 , c µ = 2µ/σ 2 , and c a = 2a/σ 2 . This admits an unnormalized solution of the form Note that finiteȳ 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

where K γ is a modified Bessel function of the second kind. The mean of the distribution is then given by

leading to the self-consistency equation,

where we define β = c D − c µ and x = √ c D c aȳ 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. 

Going beyond the numerical solution of the selfconsistency condition, we would like to gain analytic insight into the behavior of the population, particularly close to the extinction line. To do so, it is convenient to rescast Eq. (23) as

The solution now explicitly depends on only two parameters, as x = f (β, c D ), and can be plotted as in Fig. 5 . Notably, we focus on the region where x is small, hoping to construct perturbative forms for small x ∝ √ȳ .

As depicted in Fig. 5 , the small x region spans two segments; one corresponding to extinction as µ → 0 along the line c D = β > 0, and another corresponding to c D → 0 for β < 0 corresponding toȳ → 0 as a → ∞. Only the former is of interest, and for which we need the asymptotic behavior of the right hand side of Eq. (24) as x → 0. Not surprisingly, given the power-law in Eq. (20) , the expansion of K β /K β+1 for β > 0 is non-analytic, and given by

The leading term in an expansion of Eq. (24) is thus c D /β, indicating that a non-zero solution exists for c D /β > 1. However, the next order term can be either x 3 or x 1+2β depending on β, setting up two distinct singular behaviors:

• For β > 1 balancing of the cubic and linear terms leads to x 2 β(1 − β)c µ /c D . As µ → 0, the mean population size thus vanishes as

This is the usual "mean-field" singularity close to the (directed percolation) extinction transition for µ → 0, but with amplitude reduced by σ 2 /2D.

• For β < 1, the linear term must be balanced with the non-analytical correction from x 1+2β , resulting in x 2β ∝ c µ . Restoring the prefactors results in a leading singularity of the form

This is very interesting in that seascape noise with σ 2 > 2D leads to a new a universality class with the continuously varying extinction exponent of (σ 2 /2D) > 1.

Thus, excluding the exceptional point at β = 1, the singular behavior of the mean is summarized bȳ

This is confirmed in Fig. 6 , with one case having linear dependence and the other exhibiting a powerlaw singularity. To gain insight into this unusual scaling, note that averaging Eq. (16) over noise, and using y =ȳ in the stationary regime, immediately yields The ratio y 2 / y 2 is a measure of the breadth of the distribution. Keeping track of prefactors, as shown in full in Appendix A, yields

Note that this behavior smoothly connects to that obtained previously for µ ≤ 0: For 2D/σ 2 > 1, we obtain the same finite ratio of y 2 / y 2 on both sides of the transition. For 2D/σ 2 < 1, the moment ratio diverges as µ → 0, matching the infinite ratio from before. The difference in behavior is clearly confirmed in Fig. 7 . This suggests that parameter space can be divided into two regions: one where the variance is large compared to the mean, and one where they are comparable. The line that separates these regions is set by the ratio of the coupling D to the seascape noise variance σ 2 , and is schematically shown in Fig. 8 . (Distinct singular behaviors for moments or cumulants is typically associated with multifractality, and has been observed in critical be- havior of systems with quenched disorder [35, 36] .)

We may well ask if the extinction transition described above persists with the addition of demographic noise, i.e. for the equation For µ = a = 0, because of the conserved form of the deterministic part,ȳ = y = y 0 , leading to the stochastic differential equation 

Then, following the same Fokker-Plank procedure as before, we find the exact distribution for y as

with c D = 2D/σ 2 d . This is a so-called Erlang distribution with a shape parameter k = c D y 0 and rate parameter λ = c D .

Note that the tail of this distribution at large y is exponential, as opposed to the power-law in the case of seascape noise. This leads to well-behaved moments, and the moment ratio

.

As an aside, we note that for demographic noise the probability distribution is singular (albeit normalizable) for y → 0 (as may be seen in Fig. 9 ). Consequently, negative moments y −m are infinite for m > y 0 c D .

B. Demographic noise with µ < 0

Since the population decays to zero, we can ignore the quadratic term as before, and focus on the 

withȳ(t) = y 0 e −|µ|t . For D = 0, random fluctuations at a particular location can set y k = 0 (an absorbing state). At any time t, the global population is then composed of a heterogenous set of extinct nodes, and nodes with non-zero y k . The initial probability distribution thus develops a deltafunction at the origin that gradually grows to encompass the entire probability. A finite D acts as a source that feeds into the population at each location, removing the possibility of local extinction (say of an infection). This is despite the fact that this source (from the global population) is itself decaying exponentially. Unlike other extinction models (e.g., [37] ), the overall population decay does not arise from individual sub-populations going extinct, but rather from all sub-populations decaying at an even pace. Moreover, while the global decay rate is |µ|, Eq. (35) indicates that local nodes relax in response to global and noise fluctuations at a more rapid rate of (D + |µ|).

In the limit of D |µ| we can then take advantage of separation of time scales and assume that local nodes arrive at a quasi-steady state. Ignoring the slow time dependence ofȳ(t), this quasi-steady distribution is given by Eq. constructed by setting the left hand side of Eq. (35) to zero; the right hand side admits a positive solu-tion for y k (t) for any realization of noise, given y > 0.

For µ > 0, we again assume uniformity and use the law of large numbers to assertȳ = y . This reduces the problem to a single stochastic differential equation

The stationary state of the corresponding Fokker-Plank equation satisfies

where c a = 2a/σ 2 d . This admits an unnormalized solution of the form

While the distribution still diverges as a power-law close to the origin, a finiteȳ renders the PDF normalizable. For computing the normalization and moments of the distribution, we note the identity

where r = (c µ − c D )/ 2c a and z = c Dȳ /2, and 1 F 1 is the Kummer confluent hypergeometric function. While more cumbersome than in the seascape case, the mean can still be obtained from the selfconsistency condition

by numerical and perturbative analysis, which is further detailed in Appendix B. However, unlike the case of seascape noise, a finite solution forȳ does not exist for all µ > 0. Instead, a critical value of µ c should be exceeded for a finite mean, as depicted in Fig. 10 . An implicit equation determining µ c is provided in the supplementary material, in Appendix C. In the limit of small noise, it leads to µ c ∝ σ 2 d , while for large noise µ c ∝ σ d √ ln σ d . To characterize critical behavior near the extinction transition, we consider µ = µ c + δµ for δµ µ c . As indicated in the supplements, the self-consistency condition of Eq. (40) yields to an analytic expansion in the vicinity of µ c , which can be rearranged to yieldȳ ∝ δµ. This behavior is confirmed numerically in Fig. 11a .

We may also inquire about the behavior of higher moments. From averaging the equations of motion it is easy to check that y 2 = µ y /a. It can also be checked that (unlike seascape noise) all moments vanish at criticality in the same manner, i.e. y m ∝ y ∝ (δµ). Consequently, the ratio y 2 y 2 diverges as 1/(δµ) on approaching µ c , as depicted in Fig. 12a .

The novel critical behaviors at extinction observed for seascape noise thus give away to standard behavior of directed percolation with demographic noise. an independent physical motivation, a mixed-noise model is quite realistic. The two noise elements can alternatively be represented by a single noise, which after replacingȳ for the population average, leads to the stochastic differential equation

Analysis of the corresponding Fokker-Planck equation leads to a stationary distribution proportional tô

where q = 2aσ d 2 /σ 4 . Note that the (integrable) divergence close to zero (proportional to y −1+ȳc D ) is identical that of pure demographic noise. On the other hand, the exponential decay at large y is the same as that of seascape noise, although the subleading power-law is modified by an additional power of y 1−ȳc D +q . Since the extinction transition is dominated by the behavior for y close to zero, its critical behavior turns out to be the same as that of demographic noise. Indeed,ȳ = y vanishes proportionately to δµ (as seen in Fig 11b) , and the moment ratio, y 2 / y 2 , still scales inversely with δµ (as verified in Fig. 12b ). This is shown explicitly in Appendix D. Such scaling is regardless of the ratio of σ d to σ, indicating that any amount of demographic noise will destroy the unusual extinction transition observed in the seascape case. Thus the novel universality and broad distributions, characteristic of seascape noise close to an extinction transition, disappear in the presence of (expected) demographic noise.

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 y and variance y 2 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 [27, 30, 38] , while in dimensions larger than 2, a critical value of D c separates broad and narrow distributions [25, 32, 39] . 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 σ generates a power-law tail in the distribution for large populations which falls off with exponent 2+2D/σ 2 . Much like the case of DPRM at high dimensions, this signals a transition from broad to a narrow distributions for D > D c = σ 2 /2.

For a = 0, a finite µ simply makes the mean,ȳ, time dependent (proportional to e µt ) without modifying the shape of the distribution. For µ > 0, a finite a is needed to curtail the exponential growth of the population; a self-consistency requirement now enables computing a stationaryȳ. 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 µ → 0: For D > D c , the mean (and all other moments) vanish linearly with µ, as expected for a mean-field extinction transition. However, for D < D c the mean vanishes as µ σ 2 /(2D) , with higher moments vanishing with other exponents (harking back to the broad distribution expected for µ = 0).

While the above unconventional critical behavior is an accurate description of Eq. (4) as µ → 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 |µ| D, the same distribution, albeit with a time dependentȳ(t) can be used to describe growing or shrinking populations.) However, the boundary between growth and extinction is no longer at µ = 0 in the presence of demographic noise, and the self-consistency condition should be used to identify a critical value of µ c . The nature of the extinction transition at µ 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 µ c ) is again conventional. This is not too surprising as extinction characterizes behavior when y → 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 µ = 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 [40] [41] [42] [43] . With notable exceptions (e.g. Ref. [37] ), 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 selfconsistency condition in the mean-field limit, for multiple stochastic equations, generalizing Eq. (41). Even for a single quantity, we may generalize the exponent of the stabilizing nonlinear term from 2 to an arbitrary α. 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 [44, 45] . Indeed, most results presented in this paper have natural generalizations for α > 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 [16] . 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. Wellmixed 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).

We will elaborate upon the case of seascape noise with positive fitness, that is,

with µ > 0. In the main body of the paper we find that, in the appropriate limit, the distribution of y is proportional tô

where we used the law of large numbers to assert y = y . We further define c D = 2D/σ 2 , c a = 2a/σ 2 , and β = 2(D − µ)/σ 2 .

In order to get the right normalization Z for the distribution, and to choose a self-consistent value of y, we need to use the unnormalized moments

Therefore, we have that

as the self-consistent restriction on y =ȳ. Because we have assumed thatȳ = y andȳ is constant in time, then we can take the expectation of both sides of the dynamics to get µ y = a y 2 . Therefore, we have that

So if we want to understand the ratio y 2 / y 2 in the extinction limit, that is, as µ → 0 andȳ → 0, then all we need to do is understand the asymptotics of (A4) alone.

For the sake of brevity we will define x = √ c D c aȳ .

The most convenient form of K, the modified Bessel function of the second kind, is

Thus, limiting to the lowest order terms in x, the right hand side of Eq. (A4) becomes

(A7)

Disambiguating which terms are actually most relevant depends on the choice of β, with there being three regimes of note. We shall elaborate on each of these in turn.

In this case, the x −β and x −β+2 terms dominate the numerator of (A7), whereas the denominator is dominated by the x −β−1 and x −β+1 terms. Therefore, we find that

Since x ∝ √ȳ , the √ȳ on the left hand side of Eq. (A4) cancels out the leading x, hence the need to include the second order correction. From here, we can easily solve forȳ to get

Expanding β in the small-µ limit then leads tō

which is linear in µ, as stated in the main text. Moreover, we get

which is consistent with the results from the µ = 0 analysis for 2D > σ 2 .

b. −1 < β < 1

Here, the numerator of (A7) is dominated by the x β and x −β terms, and the only term from the denominator that survives is the x −β−1 term. So

The leading x cancels out theȳ on the left hand side of (A4), leving us with the mean being set bȳ

In the small µ limit, β → c D , so the lowest order behavior in µ is given bȳ

So if |β| < 1, the mean goes as µ σ 2 /(2D) . Moreover, the moment ratio is given by

(A15) This provides all the prefactors negelected in the main body Eq. (30) .

For this region, the numerator of (A7) is dominated by the x β and x β+2 terms, and the denominator is dominated by the x β+1 and x β+3 terms. So

As before, the second order expansion is required. Solving forȳ gives us

and therefore

(A18)

Unlike the previous two cases, there is no need to do a small-µ analysis here. The only way for β = c D − c µ to be negative is for µ to be greater than D. This is not possible in the limit of small µ, so this regime is irrelevant for our purposes.

Here, we will show the extinction behavior of our mean-field model under demographic noise. In particular, we will look at

(B1) As before, we use the law of large numbers to assert that the ensemble average, spatial average, and initial condition are identical. If this holds true, then there is some stationary density ρ that is identical for every y k , and we can transform the above equa- (B4) To get the moments of y, it is expedient to define z = Dȳ/σ d 2 and r = (−D + µ)/ √ aσ d 2 . Using these definitions, we find that

where 1 F 1 is the Kummer confluent hypergeometric function.

We want our analysis to take place near µ = µ c , sō y should be close to zero. The most practical way to analyze the extinction behavior is to look at both the numerator and denominator of the expression (B5) forȳ in the limit of small z, since ifȳ is small, then so is z ∝ȳ. Specficially, we havē

where A, B, E, and F are constants obtained by expanding in z, with A and B given by

and E and F given by

with Γ P being the polygamma, Erf the error function, Erfi the imaginary error function, and γ E ≈ 0.577 is the Euler-Mascheroni constant. Using these approximate forms, it is easy to see thatȳ

Let's assume that µ c is the critical value of µ wherē y goes from zero to nonzero. The above form allows us to obtain the critical µ, since it implies that µ c must satisfy c D A = 2c a , or a π

with r c = (−D + µ c )/ √ aσ d 2 . If we are instead interested in the behavior of the mean near zero, then we can substitute µ = µ c + δµ for infinitesimal δµ, and find

where

This means that our predicted mean for demographic noise scales asȳ ∝ δµ = µ − µ c near extinction, with a slope that is exactly predictable from the above equation.

By taking the expectation of both sides of the dynamics, we attain y 2 = (µ/a) y . So if we plug this into y 2 / y 2 and only keep the lowest order terms in δµ, we get

Therefore, unlike in the case of seascape noise, the moment ratio always diverges as 1/δµ. b. 0 < σ d 1

An efficient way to rewrite Eq. (C1) is

in terms of the scaled complementary error function, Erfcx. Since σ d is small, therefore r D is necessarily large. If we assume r µ is small, then we can take the first order Taylor series of Erfcx to get

Since r D is large, therefore r µ is small and selfconsistent, leading to

Another way to prove this case is by using the fact that c D is large, so we can approximate the unnormalized distribution (B4) bŷ ρ(y) y −1+c Dȳ e −c D y .

(C7)

The second moment is exactly solvable here, with y 2 =ȳ 2 +ȳ/c D . If we then take the expectation of both sides of the dynamics (B1) and plug in the second moment we find that

This suggests thatȳ converges to zero whenever µ < a/c D , which retrieves the equation for µ c (C6). So for small noise, the critical fitness scales as σ 2 d , but for large noise it scales as σ d ln(σ d ), as described in the main body of the paper.

Here, we look at the case with mixed seascape and demographic noise, that is, dy = µy − ay 2 + D (ȳ − y) dt + σ 2 y 2 + σ 2 d ydW, (D1) where σ and σ d are the respective amplitudes of seascape and demographic noises. As always, we assume uniform initial conditions and the existence of a stationary distribution ρ with a self-consistent mean y =ȳ. By transforming the above into a Fokker-Plank equation, we find the unnormalized stationary distribution ρ(y) = e −cay y −1+z (y + s) −1−t−z ,

where c a = 2a/σ 2 , s = σ 2 D /σ 2 , z = 2Dȳ/σ d 2 , q = 2aσ 2 d /σ 4 , and t = 2(D − µ)/σ 2 − q. In order to get the self-consistent mean, we need to use the following equation y = s sin(πt) π Γ(t)Γ(1 + z ) 1 F 1 (1 + z , 1 − t, q) + q t Γ(−t)Γ(1 + t + z ) 1 F 1 (1 + t + z , 1 + t, q) −Γ(z ) 1 F R 1 (z , −t, q) + q 1+t Γ(1 + t + z ) 1 F R 1 (1 + t + z , 2 + t, q)

,

with 1 F 1 being the Kummer hypergeometric function and 1 F R 1 being its regularized version. In the extinction limit,ȳ must be small, and so z should be too. In such a limit, we find y s sin(πt) π

which simplifies intō

where A = − te q q t Γ(−t)Γ(t) + e q q t Γ(−t)Γ(1 + t) 

where γ E ≈ 0.577 is the Euler-Mascheroni constant, Γ P is the polygamma, and Γ(x, y) is the incomplete gamma function ∞ y u x−1 e −u du. We can then assert that there is some µ c (with corresponding t c ) where the mean transitions from zero to nonzero. Given this, and assuming µ = µ c + δµ, with c δµ = 2δµ/σ 2 , we find for small δµ that In particular, this means that near extinction, the mean scales asȳ ∝ µ − µ c , with the scaling slope known from the above. As in the previous seascape and demographic noise cases, we can attain the second moment simply by taking the expectation of both sides of the dynamics to get y 2 = (µ/a) y . From here it is easy to see that

So despite the presence of seascape noise, the scaling of the moment ratio near extinction is identical to the purely demographic noise case.

Theoretical population biology

Absorbing Phase Transitions (Spring Science + Business Media B

Critical Dynamics: A Field Theory Approach to Equilibrium and Nonequilibrium Scaling Behavior

Proceedings of the National Academy of Sciences

Sync (Hyperion

Synchronization: a universal concept in nonlinear sciences

Containment strategy for an epidemic based on fluctuations in the sir model

The lancet infectious diseases

Appendix A: Moments for seascape noise Appendix C: µc for demographic noise Here we will find approximate forms for the critical fitness µ c , at which the mean population goes from zero to nonzero. In the prior Appendix, we demonstrated that Eq. (B12) sets the condition for µ c . If we letthen this condition can be written aswhere Erf is the error function. In the following we examine separately the large and small noise limits of the above equation.a. σ d 1Since r D ∝ 1/σ d , we might expect this quantity to be small. Let us assume that r µ is large, and futher assume that r µ r D is small. These conditions can be checked for self-consistency at the end.Taking these assumptions in hand, we haveBy assumption, we were allowed to neglect the e −2rµr D term, as well as approximate the error function by a constant. We now obtain r µ − ln(2 √ πr D ), and therefore getNotice that r µ ∝ ln(σ d ) and r µ r D ∝ ln(σ d )/σ d , making our initial assumptions self-consistent.