How to avoid an extinction time paradox

Abstract

An essential topic in theoretical ecology is the extinction of populations subject to demographic stochasticity. Mechanistic models of demographic stochasticity, such as birth-death processes, can be analytically intractable, so are frequently approximated with stochastic differential equations (SDEs). Here, we consider two pitfalls in this type of approximation. First, familiar deterministic models are not always appropriate for use in an SDE. Second, the common practice of starting directly from an SDE without explicitly constructing a mechanistic model leaves the noise term up to the modeler’s discretion. Since the stability of stochastic models depends on the global properties of both the noise and the deterministic model, overly phenomenological deterministic models, or heuristic choices of noise, can lead to models that are unrealistically stable. The goal of this article is to provide an example of how both of these effects can undermine seemingly reasonable models. Following Dennis et al. (Theor Ecol 9:323–335 2016) and Levine and Meerson (PRE, 87:032127 2013), we compare the persistence of stochastic extensions of standard logistic and Allee models. We show that, for common choices of noise, stochastic logistic models become exponentially less extinction prone when a strong Allee effect is introduced. This apparent paradox can be resolved by recognizing that common models of an Allee effect introduce overcompensation that dominates the extinction dynamics, even when the deterministic model is rescaled to account for overcompensation. These problems can be resolved by mechanistic treatment of the deterministic model and the noise.

Appendix

A.1 Extinction time estimators

In Section “Theory,” we derived an exact expression for the first passage time from carrying capacity, K, to extinction, then introduced a series of estimators based on the exact solution. These estimators were formed by Taylor expanding the potential near populations 0, A, and K. The details are provided here.

First, consider the logistic models. Let ζ = K + h for some small h < δ_ζ and 𝜖 < s < δ_s be small. Following the approach of Bresloff (2014) we can approximate the mean first passage time (15) by:

$$ \tau_{\epsilon}(x) \!\approx\! {\int}_{\epsilon}^{\delta_{s}} {\int}_{K-\delta_{\zeta}}^{K + \delta_{\zeta}} \frac{2}{\sigma^{2}} \frac{1}{v(\zeta)} \exp{\left( \frac{1}{T} \left( S(s) - S(\zeta) \right) \right)} d\zeta ds . $$

(36)

Now, since the intervals do not overlap we can separate the two integrals:

$$ \begin{array}{@{}rcl@{}} &&\tau_{\epsilon}(x) \approx {\int}_{\epsilon}^{\delta_{s}} \exp{\left( \frac{1}{T} S(s) \right)} ds \times {\ldots} \ \\ &&{\ldots} {\int}_{\delta_{\zeta}}^{\delta_{\zeta}} \frac{2}{\sigma^{2}} \frac{1}{v(K + h)} \exp{\left( \frac{1}{T} \left( - S(K + h) \right) \right)} dh . \end{array} $$

(37)

Since S(x) is minimized at x = K and is twice differentiable, it has a Taylor expansion of the form:

$$ S(K + h) = S(K) + \frac{1}{2} S^{\prime\prime}(K) h^{2} + \mathcal{O}(h^{3}). $$

(38)

Therefore, the second integral can be rewritten:

$$ \exp\left( -\frac{1}{T}S(K) \right) {\int}_{-\delta_{\zeta}}^{\delta_{\zeta}} \frac{2}{\sigma^{2}} \frac{1}{v(K + h)} \exp{\left( \frac{1}{T} \left( - \frac{1}{2 T} S^{\prime\prime}(K) h^{2} \right) \right)} dh $$

(39)

Provided S^″(K) is large enough, this integrand decays rapidly for large h. As a result, we can approximate the integral by replacing the bounds with −∞ and ∞ and replacing $\frac {1}{v(K + h)}$ with $\frac {1}{v(K)} - \frac {v^{\prime }(K)}{v(K)} h$:

$$ \frac{2}{\sigma^{2}} \frac{\exp\left( -\frac{1}{T}S(K) \right)}{v(K)} {\int}_{-\infty}^{\infty} \left[ 1 - \frac{v^{\prime}(K)}{v(K)} h \right] \exp{\left( \frac{1}{T} \left( - \frac{1}{2 T} S^{\prime\prime}(K) h^{2} \right) \right)} dh. $$

(40)

Notice that the exponential term is Gaussian. Therefore, the integral can be rewritten:

$$ \frac{2}{\sigma^{2}} \frac{\exp{\left( -\frac{1}{T} S(K)\right)}}{v(K)} \sqrt{\frac{2 \pi T}{|S^{\prime\prime}(K)|}} \mathbb{E}\left[1 - \frac{v^{\prime}(K)}{v(K)} h + \mathcal{O}(h^{2}) \right] $$

(41)

where the expected value is evaluated over a Gaussian distribution. Since the Gaussian distribution is symmetric, all odd order terms in the expansion vanish. Therefore, to second order accuracy in h:

$$ \tau_{\epsilon}(x) \approx \frac{2}{\sigma^{2}} \sqrt{\frac{2 \pi T}{|S^{\prime\prime}(K)|}} \frac{\exp{\left( -\frac{1}{T} S(K)\right)}}{v(K)} {\int}_{\epsilon}^{\delta_{s}} \exp{\left( \frac{1}{T} S(s) \right)} ds. $$

(42)

The remaining integral can be handled in much the same way. If the potential is approaching a local maximum at x = 0 then replace S(s) with a second order Taylor expansion about s = 0. This gives:

$$ \tau_{\epsilon}(x) \!\approx\! \frac{2}{\sigma^{2}} \frac{2 \pi T}{\sqrt{|S^{\prime\prime}(0)||S^{\prime\prime}(K)|}} \frac{\exp{\left( \frac{1}{T}(S(0) \!- S(K))\right)}}{v(K)}. $$

(43)

If S(s) is not approaching a local maximum as s goes to zero then the derivative of S(s) does not vanish at s = 0, so the Taylor expansion is dominated by the first order term. This approximates S(s) with an exponentially decaying function, and gives:

$$ \tau_{\epsilon}(x) \!\approx \frac{2}{\sigma^{2}} \frac{T}{|S^{\prime}(0)|}\sqrt{\frac{2 \pi T}{|S^{\prime\prime}(K)|}} \frac{\exp{\left( \frac{1}{T}(S(0) - S(K))\right)}}{v(K)}. $$

(44)

The same technique can be applied for the Allee models, expanding about s = A instead of s = 0. Since S(x) is maximized at x = A for both Allee models we approximate S(s) near A with a second order Taylor expansion to get (Bresloff 2014; Nolting and Abbott 2016):

$$ \tau_{\epsilon}(x) \!\approx\! \frac{2}{\sigma^{2}} \frac{2 \pi T}{\sqrt{|S^{\prime\prime}(A)||S^{\prime\prime}(K)|}} \frac{\exp{\left( \frac{1}{T}(S(0) - S(K))\right)}}{v(K)}. $$

(45)

These three Eqs. 43, 44, 45 all fit the general form Eq. 19.

A.2 Prefactors

In Section “Theory” we derived first passage time estimators of the general form Eq. 19. The qualitative behavior of the extinction times reflects the exponential scaling in this formula, independent of the prefactors C(K) or C(A, K). Nevertheless, these factors are important for quantitative agreement between theory and simulation (see Fig. 5). Here we compute the prefactors for each of the four models. The prefactors generally depend on the first or second derivatives of the potential at 0, A, or K. For the four models, the prefactors are:

$$ \begin{array}{lll} & \text{1a. }C(K) = \frac{2}{\sigma^{2}}\frac{1}{v(K)} \frac{2 \pi T}{\sqrt{ |S^{\prime\prime}(0)| |S^{\prime\prime}(K)|}} \\ & \text{1b. }C(K) = \frac{2}{\sigma^{2}}\frac{1}{v(K)} \frac{T}{|S^{\prime}(0)|} \sqrt{\frac{2 \pi T}{|S^{\prime\prime}(K)|}} \\ & \text{2a. }C(A,K) = \frac{2}{\sigma^{2}}\frac{1}{v(K)} \frac{2 \pi T}{\sqrt{ |S^{\prime\prime}(A)| |S^{\prime\prime}(K)|}} \\ & \text{2b. }C(A,K) = \frac{2}{\sigma^{2}}\frac{1}{v(K)} \frac{2 \pi T}{\sqrt{ |S^{\prime\prime}(A)| |S^{\prime\prime}(K)|}} . \end{array} $$

(46)

The derivatives evaluated at 0, A, and K are:

$$ \begin{array}{lll} & \text{1a. } S^{\prime\prime}(x) = \frac{1}{K}(2x - K), \quad |S^{\prime\prime}(0)| = 1, \quad |S^{\prime\prime}(K)| = 1 \\ & \text{1b. } S^{\prime}(x) = \frac{1}{K}(x - K), \quad |S^{\prime}(0)| = 1, \quad S^{\prime\prime}(x) = \frac{1}{K} = |S^{\prime\prime}(K)| \\ & \text{2a. } S^{\prime\prime}(x) = \frac{1}{A K}\left[(x - A)(K - x) + 2x - (K + A) \right], \quad |S^{\prime\prime}(A)| = \frac{1}{AK}|K - A| = |S^{\prime\prime}(K)| \\ & \text{2b. } S^{\prime\prime}(x) = \frac{1}{A K}\left[2x - (K + A) \right], \quad |S^{\prime\prime}(A)| = \frac{1}{AK}|K - A| = |S^{\prime\prime}(K)| . \end{array} $$

(47)

Therefore, the prefactors are:

$$ \begin{array}{lll} & \text{1a. }C(K) = \frac{4 \pi T}{\sigma^{2}} \\ & \text{1b. }C(K) = \frac{2 T}{\sigma^{2}} \sqrt{\frac{2 \pi T}{K}} \\ & \text{2a. }C(A,K) = \frac{4 \pi T}{\sigma^{2}} \frac{ A K}{K - A} \\\ & \text{2b. }C(A,K) = \frac{4 \pi T}{\sigma^{2}} \frac{A}{K - A} . \end{array} $$

(48)

Or, substituting $T = \frac {\sigma ^{2}}{2 r}$:

$$ \begin{array}{lll} & \text{1a. }C(K) = \frac{2 \pi }{r} \\ & \text{1b. }C(K) = \frac{1}{r} \sqrt{\frac{2 \pi T}{K}} \\ & \text{2a. }C(A,K) = \frac{2 \pi }{r} \frac{ A K}{K - A} \\ & \text{2b. }C(A,K) = \frac{2 \pi }{r} \frac{A}{K - A} . \end{array} $$

(49)

Notice that the prefactors for both Allee models diverge if A approaches K. The prefactor diverges because models 2a and 2b have quartic and cubic potentials (respectively). As A approaches K the second derivative of the potential at A and K vanishes, so a second order Taylor expansion about S(A) and S(K) is not accurate. Therefore we restrict to A ≲ K/2. More accurate approximations could be made by carrying the Taylor expansions out to higher order; however, since the potentials are polynomials this is equivalent to solving for the first passage times exactly.

A.3 Rescaling

In Section “Analysis,” we discussed three rescalings of m(x). This Appendix derives the necessary rescalings, and analyzes their effect on the persistence of the Allee model. In general, consider r_A = s(K, A)r_L where s(K, A) is some scaling factor. We will consider s(K, A) such that the two models have the same:

• 1. maximum per capita growth rate (Gruntfest et al. 1997),

• 2. maximum absolute growth rate, or

• 3. linear stability at carrying capacity.

The per capita growth rate for the logistic model is r_L(K − x)/K and is maximized at x_∗ = 0. Therefore, the maximum per capita growth rate for the logistic model is r_L. The per capita growth rate for the Allee model is the quadratic r_A(x − A)(K − x)/(AK) and is maximized halfway between the Allee threshold and the carrying capacity at x_∗ = (K + A)/2. Therefore, the maximum per capita growth rate for the Allee model is r_L(K − A)²/(4AK). It follows that the Allee model has the same maximum per capita growth rate as the logistic model if $s(A,K) = 4 \frac {A K}{(K-A)^{2}}$.

Plugging into Eq. 25, the arguments of the exponential term in the extinction times (work to extinction divided by temperature) become:

$$ \begin{array}{lll} & \text{2a. } \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \frac{2 r_{L}}{\sigma^{2}} \frac{4 AK}{(K-A)^{2}} \frac{1}{6 A K}\left[A^{3} \left( K - \frac{1}{2} A \right) - K^{3} \left( A - \frac{1}{2} K \right) \right] \\ & \text{2b. } \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \frac{2 r_{L}}{\sigma^{2}} \frac{4 AK}{(K-A)^{2}} \frac{1}{2 A K}\left[A^{2} \left( K - \frac{1}{3} A \right) - K^{2} \left( A - \frac{1}{3} K \right) \right] \end{array} $$

(50)

For large K − A, both of these polynomials are dominated by the second term in the brackets, so taking the limit as K goes to infinity:

$$ \begin{array}{lll} & \text{2a. } \lim\limits_{K \rightarrow \infty} \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{2 r_{L}}{\sigma^{2}} \frac{4 AK}{(K-A)^{2}} \frac{1}{6 A K} K^{3} \left( \frac{1}{2} K - A \right) \\ & \text{2b. } \lim\limits_{K \rightarrow \infty} \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{2 r_{L}}{\sigma^{2}} \frac{4 AK}{(K-A)^{2}} \frac{1}{2 A K} K^{2} \left( \frac{1}{3} K - A \right) \end{array} $$

(51)

Replacing K − A with K, $\frac {1}{2} K - A$ with $\frac {1}{2} K$, and canceling constants gives:

$$ \begin{array}{lll} & \text{2a. } \lim\limits_{K \rightarrow \infty} \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{2}{3} K^{2} \\ & \text{2b. } \lim\limits_{K \rightarrow \infty} \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{4}{3} K. \end{array} $$

(52)

Notice that both of these arguments are one order lower when compared to the same terms in the unscaled models. On the other hand, the arguments of the exponential terms in the logistic models are:

$$ \begin{array}{lll} & \text{1a. } \lim\limits_{K \rightarrow \infty} \frac{2 r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{1}{3} K^{2} \\ & \text{1b. } \lim\limits_{K \rightarrow \infty} \frac{2 r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} K. \end{array} $$

(53)

Therefore, the argument of the exponential term in the Allee models converges to twice the argument of the exponential term in the logistic model when the noise is constant, and four-thirds the argument when the noise is linear. The prefactors for the rescaled Allee models are $C(A,K) = \frac {2 \pi }{r_{L}} (K-A)$ and $C(A,K) = \frac {2 \pi }{ r_{L}} \frac {K-A}{K}$. As usual these are rational in K and A, so are dominated by the exponential terms when K is large. Therefore, while the scaling successfully reduces the order of the argument of the exponential term in the Allee models, it does not resolve the paradox.

The same analysis can be repeated for the second rescaling. The logistic growth rate is quadratic and is maximized at x_∗ = K/2. Therefore, the maximum absolute growth rate for the logistic model is r_LK/4. The Allee growth rate is cubic, and is maximized at $x_{*} = \frac {1}{3}[(K+A) + \sqrt {K^{2} + KA - A^{2}}]$. For large K − A, this approaches $x_{*} \approx \frac {2}{3} K + \frac {1}{3} A$. Therefore, the maximum absolute growth rate of the Allee model is approached (from below) by $\frac {2 r_{A}}{27 AK} (K-A)^{2} (2K + A)$. Therefore, the scaling s(A, K) approaches $\frac {27}{8} \frac {A K^{2}}{(K - A)^{2}(2K + A)}$ from below when K is large.

Plugging into Eq. 25 and taking the limit of large K gives:

$$ \begin{array}{lll} & \text{2a. } \lim\limits_{K \rightarrow \infty} \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{27}{96} K^{2} = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{9}{32} K^{2} \\ & \text{2b. } \lim_{K \rightarrow \infty} \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{27}{48} K = \lim_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{9}{16} K. \end{array} $$

(54)

Once again, rescaling lowers the order of the argument of exponential term in the Allee models. Compared to the equivalent logistic models, the argument of the exponential term in the Allee models converge to 28/32 = 7/8 times the argument in the logistic model with constant noise, and 9/16 times the argument in the logistic model with linear noise. Since both of these factors are less than one, the exponential term in the Allee model is smaller than in the logistic model, resolving the paradox when K is large enough that C(A, K) can be ignored.

Equating the linear stability of the two models at carrying capacity requires rescaling the Allee model so that the derivative of the absolute growth rate at x = K is equal to the derivative of the absolute growth rate of the logistic model at x = K. At x = K the logistic growth rate has slope − r_L. At x = K the growth rate of the Allee model has slope $-r_{A} \frac {K-A}{A}$. Therefore $s(A,K) = \frac {A}{K-A}$.

Plugging into Eq. 25 and taking the limit of large K gives:

$$ \begin{array}{lll} & \text{2a. } \lim\limits_{K \rightarrow \infty} \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{1}{6} K^{2} \\ & \text{2b. } \lim\limits_{K \rightarrow \infty} \frac{2 s(A,K) r_{L}}{\sigma^{2}} W(K \rightarrow 0) = \lim\limits_{K \rightarrow \infty} \frac{ r_{L}}{\sigma^{2}} \frac{1}{3} K . \end{array} $$

(55)

The comparable terms in the logistic model are twice and three times as large; therefore, the logistic model is more stable when K is large enough that C(A, K) can be ignored.

A.4 Diffusion approximation

In Section “How to avoid the paradox” we showed that the extinction time paradox is resolved by deriving the SDE models from birth-death processes, without discussing the details of the approximations involved. This is a well studied area, and there are multiple ways of constructing such an approximation. Accordingly, this appendix will distinguish which of these methods is most appropriate in the context of this paper.

The most familiar method for approximating a birth-death process with an SDE is a system size expansion (Van Kampen 1992). The system size expansion assumes that both the birth and death parameter scale naturally in some system size, Ω taken to represent a characteristic population. In our case Ω could be set to either the carrying capacity or the Allee threshold. By rescaling the population in the system size it is often possible to perform an asymptotic expansion of the birth-death process. As Ω becomes large, the master equation converges to a Fokker-Planck equation (Van Kampen 1992). This Fokker-Planck equation governs the diffusion of probability for an SDE with instantaneous drift m(x) = λ(x) − μ(x) and diffusion coefficient $v(x) = \frac {1}{\Omega }(\lambda (x) + \mu (x))$ (Bresloff 2014).

Notice that this method assumes that the system size is large. As a result, it is only appropriate if the dynamics of interest occur at large populations. This is not the case for an extinction process since a population must necessarily be small before going extinct. Albeit, if the Allee threshold, A, is large, then the extinction process can be separated into two phases. First, the population escapes from the carrying capacity to the Allee threshold; then, after crossing the Allee threshold, it goes extinct with high probability. In general, the escape process is much slower than the process of descending from the Allee threshold to extinction. Therefore, the mean extinction time of an Allee model is dominated by the model’s behavior at populations near, or greater than, the Allee threshold. In that case, if A is large, it is possible to approximate the mean extinction time based solely on the behavior at large populations. Then, a system size expansion is appropriate. However, our goal is to compare extinction times between Allee models and logistic models. The logistic model does not have a threshold past which extinction becomes highly likely. As a result, mean extinction times for the logistic model depend critically on the behavior of the model at small populations, even if the carrying capacity is large. Therefore, a system size expansion is not appropriate. In some cases, it is possible to use a system size expansion coupled to a small system correction (Ovaskainen and Meerson 2010); however, this is unnecessarily complicated in the context of this paper.

Also notice that the mean drift and diffusion coefficient scale differently in Ω. Unlike the mean, the diffusion coefficient is vanishing as Ω becomes large (Bresloff 2014).⁷ As a result, large system size implies low temperature, and long extinction times. This limit is counterproductive since our goal was to study the roles of K, A, and temperature separately. Moreover, if the system size is greater than one it leads to an SDE that is artificially less noisy than the original birth-death process. One solution is to take the system size limit to derive expressions for m(x) and v(x) in terms of λ(x), μ(x) and Ω, then set Ω = 1. This, however, defeats the purpose of the system size expansion since the system size expansion is only accurate⁸ if it is assumed that Ω is large.

Lastly, there exists a discrete potential for any birth-death process (see Eq. 64) analogous to the continuous potentials discussed in Section “Analysis.” This discrete potential does not match the continuous potential at small populations when the continuous potential is derived from a system size expansion (Bresloff 2014). As a result, both the potentials and the passage time statistics of the SDE differ from the potential and passage time statistics of the original birth-death process (Doering et al. 2005).

In summary, the system size approximation constructs a family of discrete birth-death processes indexed by a parameter that converge to an SDE in a limit. The discrete process is closely approximated by an SDE when the parameter is near its limiting value. This ensures that the SDE is a good approximation for the system in a particular limit; however, it does not guarantee that the birth-death process in that limit closely approximates the original birth-death process. In particular, there is no guarantee that a birth-death process with a large system size gives a good approximation for a birth-death process at small or finite system sizes.

An alternative approach is proposed in Allen and Allen (2002) and Cresson and Sonner (2018). Instead of approximating the discrete process in a particular limit where it converges to an SDE, we could attempt to approximate the original process with an SDE without taking a limit. There is no guarantee that this approximation will be accurate, since the discrete process is not an SDE; however, this approach avoids introducing errors due to an artificial system size. In fact, since the discrete process is not an SDE there will necessarily be error in any approximation, so the choice of SDE is not unique. In Allen and Allen (2002) Allen proposes a tried-and-true heuristic for picking a SDE that is conceptually consistent with the original birth-death process, which we now describe.

Any SDE is uniquely specified by its mean drift m(x) and diffusion coefficient v(x). Let $\bar {x} = \mathbb {E}[X]$. If the initial state of the system is known then $m(x) = \frac {d}{dt} \mathbb {E}[X]$ and $v(x) = \frac {d}{dt} \mathbb {E}[(X - \bar {x})^{2}]$. Given the initial state it is possible to compute both the rate of change in the expected population, and the rate of change of the diffusion coefficient in population for a birth-death process. Therefore, it is natural to approximate the birth-death process with an SDE that has the same mean drift and diffusion coefficient.

Suppose the discrete process is in state x at time 0. Then, from the master equation, it is easy to check that:

$$ \begin{array}{lll} & m(x) = \frac{d}{dt} \mathbb{E}[X|X(0) = x] = \lambda(x) - \mu(x) \\ & v(x) = \frac{d}{dt} \mathbb{E}[(X - x)^{2}|X(0) = x] = \lambda(x) + \mu(x). \end{array} $$

(56)

This motivates (29). This is philosophically different than the system size approximation since it uses an SDE to approximate the discrete birth-death process as it is, not as it is in a particular limit. When formulated in this way the passage time statistics for the SDE and the birth-death processes are generally similar (Allen and Allen 2002).⁹

Like the system size expansion, it is possible to construct a family of discrete birth-death processes that converge to this SDE in a particular limit. Suppose that we refine the discrete process by introducing fractional populations. By modifying the birth and death rates at the same time, it is possible to construct a sequence of refined birth-death processes whose mean drift and diffusion coefficient do not depend on the size of the refinement.

Consider a discrete birth-death process where each event produces 0 ≤Δx ≤ 1 individuals. Assume that the birth and death rates are smooth continuous functions of x, and let p(x, t) denote the probability that the population has x individuals at time t. Define the modified rates $\tilde {\lambda }(x|{\Delta } x)$ and $\tilde {\mu }(x|{\Delta } x)$ such that $\tilde {\lambda }(x|1) = \lambda (x)$ and $\tilde {\mu }(x|1) = \mu (x)$. For concision, we will repress the dependence on Δx unless necessary. Then, the refined master equation at x reads:

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} p(x,t) &=& \tilde{\lambda}(x - {\Delta} x) p(x-{\Delta} x, t) - (\tilde{\lambda}(x)) + \tilde{\mu}(x)) p(x,t)\\ &&+ \tilde{\mu}(x + {\Delta} x) p(x + {\Delta} x, t). \end{array} $$

(57)

Now let Δx be small. Taylor expanding in Δx gives:

$$ \begin{array}{@{}rcl@{}} \begin{array}{lll} && \tilde{\lambda}(x - {\Delta} x) p(x-{\Delta} x, t) = \tilde{\lambda}(x) p(x,t) - \frac{\partial}{\partial x} \tilde{\lambda}(x) p(x,t) {\Delta} x \\&&\quad+ \frac{1}{2} \frac{\partial^{2}}{\partial x^{2}} \tilde{\lambda}(x) p(x,t) {\Delta} x^{2} + \mathcal{O}({\Delta} x^{3}) \\ && \tilde{\mu}(x + {\Delta} x) p(x+{\Delta} x, t) = \tilde{\mu}(x) p(x,t) + \frac{\partial}{\partial x} \tilde{\mu}(x) p(x,t) {\Delta} x \\&&\quad+ \frac{1}{2} \frac{\partial^{2}}{\partial x^{2}} \tilde{\mu}(x) p(x,t) {\Delta} x^{2} + \mathcal{O}({\Delta} x^{3}). \end{array} \end{array} $$

(58)

The zeroth order terms cancel, leaving:

$$ \begin{array}{@{}rcl@{}} &&\frac{d}{dt} p(x,t) = - \frac{\partial}{\partial x} \left( (\tilde{\lambda}(x) - \tilde{\mu}(x)) p(x,t) \right) {\Delta} x \\ &&\quad+ \frac{1}{2} \frac{\partial^{2}}{\partial x^{2}} \left( (\tilde{\lambda}(x) + \tilde{\mu}(x)) p(x,t) \right) {\Delta} x^{2} + \mathcal{O}({\Delta} x^{3}). \end{array} $$

(59)

Finally, to rewrite in terms of density, ρ(x, t), divide across by Δx:

$$ \begin{array}{@{}rcl@{}} &&\frac{\partial}{\partial t} \rho(x,t) = - \frac{\partial}{\partial x} (\tilde{\lambda}(x) - \tilde{\mu}(x)) \rho(x,t) + \frac{1}{2} \frac{\partial^{2}}{\partial x^{2}} (\tilde{\lambda}(x) \\&&\quad+ \tilde{\mu}(x)) \rho(x,t) {\Delta} x + \mathcal{O}({\Delta} x^{2}). \end{array} $$

(60)

In the limit of vanishing Δx, Eq. 60 approaches the Ito form of the Fokker-Planck equation with infinitesimal mean and diffusion coefficient:

$$ \begin{array}{lll} & m(x|{\Delta} x) = \tilde{\lambda}(x|{\Delta} x) - \tilde{\mu}(x|{\Delta} x) \\ & v(x|{\Delta} x) = {\Delta} x \left( \tilde{\lambda}(x|{\Delta} x) + \tilde{\mu}(x|{\Delta} x)\right). \end{array} $$

(61)

Now, to ensure the noise intensity, σ² does not vanish as Δx goes to zero set:

$$ \begin{array}{lll} & \tilde{\lambda}(x|{\Delta} x) = \frac{1}{2 {\Delta} x} \left( \lambda(x) + \mu(x) \right) + \frac{1}{2} \left( \lambda(x) - \mu(x) \right) \\ &\tilde{\mu}(x|{\Delta} x) = \frac{1}{2 {\Delta} x}\left( \lambda(x) + \mu(x) \right) - \frac{1}{2} \left( \lambda(x) \!- \mu(x) \right). \end{array} $$

(62)

Then, $\tilde {\lambda }(x|1) = \lambda (x)$, $\tilde {\mu }(x|1) = \mu (x)$ and the infinitesimal mean and diffusion coefficient of the SDE match the infinitesimal mean and diffusion coefficient of the original birth-death process. Defining the modified birth and death rates in this way ensures that the infinitesimal mean and diffusion coefficient of the refined process scale indentically in Δx. This is accomplished whenever:¹⁰

$$ \frac{\tilde{\lambda}(x|{\Delta} x) - \tilde{\mu}(x|{\Delta} x)}{\tilde{\lambda}(x|{\Delta} x) + \tilde{\mu}(x|{\Delta} x)} = \mathcal{O}({\Delta} x). $$

(63)

Condition (63) guarantees that the potential of the birth-death process converges to the potential of the limiting SDE (Allen and Allen 2002). This is not true for all diffusion approximations (Doering et al. 2005), including the systems size expansion (Bresloff 2014). See Appendix A.5 for an example of a discrete birth-death process that, when refined, obeys condition (63) automatically.

A discrete birth-death process has potential:

$$ \phi(x) = -\frac{1}{2}\sum\limits_{y = 0}^{x} \log{\left( \frac{\tilde{\lambda}(y)}{\tilde{\mu}(y + {\Delta} x)} \right)} $$

(64)

and mean times to extinction (Allen and Allen 2002; Leigh 1981):

$$ \tau(K) = \sum\limits_{x = 0}^{K} \sum\limits_{y = x}^{\infty} \exp{\left( 2(\phi(x) - \phi(y)) \right)}. $$

(65)

Notice the similarity between (65) and (15). Taking Δx to zero:

$$ \frac{\tilde{\lambda}(x)}{\tilde{\mu}(x + {\Delta} x)} = \frac{v(x) + m(x) {\Delta} x}{v(x) - m(x) {\Delta} x} = 1 + 2 \frac{m(x)}{v(x)} {\Delta} x + \mathcal{O}({\Delta} x^{2}). $$

(66)

Therefore:

$$ \frac{1}{2}\log{\left( \frac{\tilde{\lambda}(x)}{\tilde{\mu}(x + {\Delta} x)} \right)} = \frac{m(x)}{v(x)} {\Delta} x + \mathcal{O}({\Delta} x^{2}). $$

(67)

Then:

$$ \begin{array}{@{}rcl@{}} \lim\limits_{{\Delta} x \rightarrow 0} \sum\limits_{i = 0}^{j} \frac{1}{2}\log{\left( \frac{\tilde{\lambda}(x)}{\tilde{\mu}(x + {\Delta} x)} \right)} = \lim\limits_{{\Delta} x \rightarrow 0} \sum\limits_{i = 0}^{j} \frac{m(x)}{v(x)} {\Delta} x \\+ \mathcal{O}({\Delta} x^{2}) = {{\int}_{0}^{x}} \frac{m(s)}{v(s)} ds = {{\int}_{0}^{x}} f(s) ds. \end{array} $$

(68)

Now, from Eq. 11, it follows that the discrete potential converges to the continuous potential:

$$ \lim_{{\Delta} x \rightarrow 0} \phi = S. $$

(69)

This same technique is used in Leigh (1981) to approximate ϕ analytically. As Δx goes to zero, the double sum is replaced with a double integral and the extra factor of 2 is absorbed into the temperature. Then, the mean extinction time for the discrete model converges to the mean extinction time for the continuous model.

A.5 Physical analogy

In Section “Theory,” we asserted that the ratios f(x) = m(x)/v(x) and $T = \frac {\sigma ^{2}}{2 r}$ acted like forces and temperature. Here, we develop a specific physical system which is mathematically identical to the population models we consider. We then show that f(x) and T are force and temperature in the analogous physical system.

Consider an ion with charge q diffusing down a narrow channel. The channel could be a pore in a cell membrane. The ion is subject to constant thermal noise and, in the absence of a driving potential, would follow a Wiener process (Brownian motion). Assume that the channel is kept at fixed temperature T, and there is only one charged particle in the channel. The channel is not kept at constant voltage, so for any x, there is a voltage V (x). Therefore, if the ion is at x it has electrostatic potential energy U(x) = qV (x).

Coarse grain the channel into a sequence of discrete intervals with width Δx. These could correspond to a sequence of binding sites that are evenly spaced along the channel. Then, the motion of the ion is closely approximated by a discrete random walk with exponentially distributed waiting times. For a given x let λ(x,Δx), denote the rate at which the ion hops to x + Δx and μ(x,Δx) denote the rate at which the ion hops to x −Δx. For concision, we will suppress λ and μ’s dependence on Δx except where necessary. Following (Schnakenberg 1976):

$$ \frac{\lambda(x)}{\mu(x + {\Delta} x)} = \exp{\left( \frac{1}{k_{B} T} (U(x) - U(x + {\Delta} x)) \right)} $$

(70)

where k_B is the Boltzmann constant. This is the essential physics that will link forces and temperature to ratios of m(x), v(x), σ², and r. In order to study the relationship between physical quantities and the SDE, we take a diffusion limit one side at a time.

The derivative of the potential energy U(x) is a force f(x), so, by the intermediate value theorem:

$$ U(x) - U(x + {\Delta} x) = - {\int}_{x}^{x + {\Delta} x} f(s) ds = f(\zeta) {\Delta} x $$

(71)

for some ζ ∈ [x, x + Δx]. In the limit as Δx goes to zero f(ζ) converges to f(x). Now, Taylor expanding the exponential in small Δx:

$$ \exp{\left( -\frac{1}{k_{B} T} f(\zeta) {\Delta} x \right)} = 1 + \frac{1}{k_{B} T} f(\zeta) {\Delta} x + \mathcal{O}({\Delta} x^{2}). $$

(72)

This implies that, in the limit as Δx goes to zero the ratio λ(x,Δx)/μ(x + Δx,Δx) converges to 1. This, in turn, implies that (λ(x,Δx) − μ(x + Δx,Δx))/(λ(x,Δx) + μ(x + Δx,Δx)) is $\mathcal {O}({\Delta } x)$ (converges to zero proportional to Δx). Therefore, this discrete process is mathematically identical to the birth-death process considered in Appendix A.3. It follows that:

$$ \frac{\lambda(x)}{\mu(x + {\Delta} x)} = \frac{\sigma^{2} v(x) + r m(x) {\Delta} x}{\sigma^{2} v(x) - r m(x) {\Delta} x} = 1 + 2 \frac{r m(x)}{\sigma^{2} v(x)} {\Delta} x + \mathcal{O}({\Delta} x^{2}). $$

(73)

Equating the left and right hand sides:

$$ \frac{2 r}{\sigma^{2}} \frac{m(x)}{v(x)} = \frac{1}{k_{B} T} f(x) $$

(74)

This leads immediately to the association:

$$ \begin{array}{lll} & \frac{m(x)}{v(x)} = f(x) \\ & \frac{\sigma^{2}}{2 r} = k_{B} T. \end{array} $$

(75)

A.6 Analysis for a saturating Allee effect

In Section “How to avoid the paradox,” we proposed two different modifications to our Allee model that resolve the passage time paradox: either modify $m(x) = \frac {r}{K}x(K - x)$ with a saturating function or derive v(x) from a birth-death model. Here, we show that these two modifications lead to a potential with the same functional form.

First, suppose we set:

$$ \begin{array}{lll} & m(x) = \frac{r}{K}x(K - x) - \beta (1 - P(x)) x, \quad P(x) = \frac{x}{\theta + x} \end{array} $$

(76)

as in Eq. 27. This model can be derived by assuming fecundity β, probability of finding a mate P(x), background per capita death rate δ₀ and additional death rate due to competition δ₁ (Boukal and Berec 2002; Dennis 2002; 1989; McCarthy 1997). To convert into r and K set $K = \frac {\beta - \delta _{0}}{\delta _{1}}$ and r = β − δ₀ (see Eq. 27).

The corresponding deterministic model has root where m(x) = 0. Therefore, the roots occur at 0, some A, and some K^′ near K. Since the subsequent model will be analyzed in terms of birth and death rates, and since the Allee threshold A and carrying capacity K^′ are no longer explicit parameters of the model it will be convenient to work with the parameters β, δ₀,δ₁,𝜃 from now on. The corresponding birth-death model is:

$$ \begin{array}{lll} & \lambda(x) = \beta \frac{x^{2}}{\theta + x} \\ & \mu(x) = \delta_{0} x + \delta_{1} x^{2}. \end{array} $$

(77)

Therefore, the forces are:

$$ f(x) = \frac{\beta \frac{x}{\theta + x} - \delta_{0} - \delta_{1} x}{\beta \frac{x}{\theta + x} + \delta_{0} + \delta_{1} x} $$

(78)

which can be simplified by multiplying through by 𝜃 + x:

$$ \begin{array}{@{}rcl@{}} f(x) &&= \frac{\beta x - \delta_{0} (\theta + x) - \delta_{1} x (\theta + x) }{\beta x + \delta_{0} (\theta + x) + \delta_{1} (\theta + x)} \\&&= - \frac{(-\beta + \delta_{0} + \delta_{1} \theta) x + \delta_{0} \theta + \delta_{1} x^{2} }{(\beta + \delta_{0} + \delta_{1} \theta) x + \delta_{0} \theta + \delta_{1} x^{2})}. \end{array} $$

(79)

Let $b = \frac {\beta }{\delta _{1}}$, $c = \frac {\delta _{0}}{\delta _{1}} \theta$, and $d = \frac {\delta _{0} + \delta _{1} \theta }{\delta _{1}}$. Then, after dividing through by δ₁ and rearranging:

$$ f(x) = - \frac{ x^{2} + (d - b) x + c }{x^{2} + (b + d) x + c}. $$

(80)

The corresponding potential S(x) is:

$$ S(x) = {{\int}_{0}^{x}} \frac{ s^{2} + (d - b) s + c }{s^{2} + (b + d) s + c} ds. $$

(81)

To integrate, replace the numerator with s² + (d + b)s + c − 2bs. Then:

$$ \begin{array}{@{}rcl@{}} S(x) &&= {{\int}_{0}^{x}} 1 - \frac{2 b s}{s^{2} + (b + d) s + c} ds \\&&= x - 2b {{\int}_{0}^{x}} \frac{s}{s^{2} + (b + d) s + c} ds. \end{array} $$

(82)

Notice that the integrand is vanishing for large s, so the potential is close to linear in x for large x. This is a natural feature of f(x) derived from a birth-death model, and reflects the saturating behavior of f(x). To perform the integral factor the denominator:

$$ s^{2} + (b + d) s + c = (s - r_{1})(s - r_{2}) $$

(83)

where r₁,r₂ are the roots of the polynomial given by:

$$ r_{1,2} = -\frac{(b+d) \pm \sqrt{(b+d)^{2} - 4 c}}{2}. $$

(84)

Notice that since b and d are both positiv, the real part of the roots is always negative.

Next, we can expand the integrand using partial fractions:

$$ \frac{s}{s^{2} + (b + d) s + c} = \frac{r_{1}}{r_{1} - r_{2}}\frac{1}{s - r_{1}} + \frac{r_{2}}{r_{2} - r_{1}}\frac{1}{s - r_{2}}. $$

(85)

Now, the integration is easy:

$$ \frac{r_{1}}{r_{1} - r_{2}} {{\int}_{0}^{x}}\frac{1}{s - r_{1}} ds + \frac{r_{2}}{r_{2} - r_{1}}{{\int}_{0}^{x}} \frac{1}{s - r_{2}} ds = \frac{r_{1}}{r_{1} - r_{2}}\log{\left( - \frac{(x - r_{1})}{r_{1}}\right)} + \frac{r_{2}}{r_{2} - r_{1}}\log{\left( - \frac{(x - r_{2})}{r_{2}} \right)}. $$

(86)

Therefore:

$$ S(x) = x - \left[ \frac{r_{1}}{r_{1} - r_{2}}\log{\left( - \frac{(x - r_{1})}{r_{1}}\right)} + \frac{r_{2}}{r_{2} - r_{1}}\log{\left( - \frac{(x - r_{2})}{r_{2}} \right)} \right]. $$

(87)

Notice the immediate similarity with Eq. 35. This is actually a general form, since any forces that are rational functions whose numerator and denominator are of the same order can be analyzed in the same way. Note that if v(x) is derived from an underlying birth-death model, the birth and death rates are polynomials or rational functions, and the higher order terms in the birth and death rates do not cancel when added or subtracted, then the forces are rational with numerator and denominator of the same order. For example, suppose we had started with the birth-death model (Brassil 2001):

$$ \begin{array}{lll} & \lambda(x) = \beta_{1} x^{2} \\ & \mu(x) = \delta_{0} x + \delta_{2} x^{3} \end{array} $$

(88)

Then:

$$ \begin{array}{lll} & m(x) = \beta_{1} x^{2} - \delta_{0} x - \delta_{2} x^{3} \\ & v(x) = \beta_{1} x^{2} + \delta_{0} x + \delta_{2} x^{3}. \end{array} $$

(89)

Notice that m(x) is cubic, so could be rewritten $m(x) = \frac {r}{AK}(x - A)(K-x) x$ as in models 2a and 2b. The corresponding forces are:

$$ f(x) = -\frac{ \delta_{2} x^{2} - \beta_{1} x + \delta_{0}}{\delta_{2} x^{2} + \beta_{1} x + \delta_{0}} $$

(90)

Let $b = \frac {\beta _{1}}{\delta _{2}}$, $c = \frac {\delta _{0}}{\delta _{2}}$. Then:

$$ f(x) = -\frac{x^{2} - b x + c}{x^{2} + b x + c} = -1 + 2 b \frac{x}{x^{2} + b x + c} $$

(91)

This is essentially the same as the forces we analyzed in the previous case, only with different coefficients b and c. The potential can be derived following the same steps as before, only with roots:

$$ r_{1,2} = -\frac{b \pm \sqrt{b^{2} - 4 c}}{2}. $$

(92)

In either case, the potential is linear in x with a logarithmic correction. It follows that the passage time from carrying capacity to extinction is, at worst, exponential in K.