The Role of Stochasticity in Noise-Induced Tipping Point Cascades: A Master Equation Approach

Many systems in nature can transition into a qualitatively distinct dynamic regime when a critical threshold is approached. The associated threshold of such a system, defined in the context of bifurcation theory, is the bifurcation point or tipping point of the system. Tipping points manifest in systems including lake eutrophication (ecology), social contagion (sociology), disease spread (epidemiology), epileptic seizures (physiology), stock market crashes (finance), and even the earth system (climate science) (Scheffer, 2009; Scheffer et al., 2012; O’Regan, 2018; Lenton, 2020; Klose et al., 2020). Tipping points arise in the presence of strongly self-amplifying (mathematically positive) feedbacks (Lenton, 2020). This characterization of tipping points is important in the sense that only sufficiently strong, self-propelling feedbacks are recognized.

A prime example of a tipping point is an ecological system with an Allee effect (Courchamp et al., 2008; Hastings and Gross, 2012; Johnson and Hastings, 2018; Vortkamp et al., 2020). Allee effects occur in populations with low abundances and are believed to be commonplace in ecological systems (Courchamp et al., 2008; Drake and Kramer, 2011). In a population exhibiting an Allee effect, the per capita growth rate is a unimodal function of the population abundance with a global maximum. The sign of the growth rate at low population levels distinguishes a weak Allee effect from a strong one. In particular, a weak Allee effect does not result in a negative growth rate for small population sizes but its strong counterpart presents with a negative rate. Therefore, the existence of a strong Allee effect implies the existence of a critical threshold for survival. The possibility of alternative stable states in systems that are analogous to those with an Allee effect has significant implications at many levels, from the microscopic scale of budding yeast (Dai et al., 2012) to the macroscopic scale of tipping elements for the Earth system (Lenton, 2020; Klose et al., 2020).

In order to investigate tipping points induced by sources of variation in systems with an Allee effect, one needs to formulate a model that can exhibit either a weak or strong Allee effect. As the quality of the environment is degraded, a weak Allee effect becomes a strong one through a transcritical bifurcation. A fundamental question in this setting concerns the propagation of tipping points through an ecosystem consisting of multiple patches. To what extent are interacting populations interdependent and how is this relationship influenced by the parameters governing the model behavior? Stated more simply, how do tipping points cascade through systems?

A simple analysis considers multiple populations in a network that are connected by passive, symmetric diffusion. In this work, we study a model consisting of two populations to analyze how their dynamics are related, conditional on the quality of their internal environments. One case of interest is when both populations are above the Allee threshold to see when the collapse of either population is less likely than for isolated systems. A second case of interest is when one population has fallen below the Allee threshold and we ask whether the next transition puts both populations above, or alternatively, below, the Allee threshold. A similar study was conducted in a deterministic setting (Johnson and Hastings, 2018). That study laid the foundations for an eventual treatment in a stochastic setting. Hence, in our work, we adopt the definition of ecological resilience formulated by Holling (Holling, 1973) as well as the framework proposed by Johnson and Hastings, but we also account for stochasticity in the system.

The literature on single, isolated tipping points is vast. However, studies of cascades of tipping points are less common. We begin with the background needed to formulate a stochastic, individual-based model that accounts for Allee effects of varying strengths. Our model follows a Markovian birth-death process, rooted in the theory of chemical reaction networks. We find that this representation naturally lends itself to a treatment with Markov chains and first passage times. Then, we employ dimensionality reduction techniques and various approximations to derive a reduced model that we study in detail. The section that follows consists of the results obtained using the above methods. Finally, we summarize our findings with conclusions.

Much of the deterministic modeling work done on systems with one population that manifest Allee effects are based on either phenomenology or empirical observations. The crux of these models dates back to 1954 (Odum and Allee, 1954), where the observed per capita growth rate was fit using a suitable function. The general form of a deterministic model in continuous time is the following ordinary differential equation:

$$\dfrac{d\rho(t)}{dt}=\rho f(\rho),$$

(2.1)

where $\rho(t)$ denotes the average number of individuals at time $t,$ and $f(\rho)$ is a function specifying the form of the per-capita population growth rate at size $\rho$. There have been many functional forms proposed for $f(\rho)$ in the literature. A review of various specifications used in models can be found (Boukal and Berec, 2002).

In particular, a simple and important model was proposed early in the twentieth century (Volterra, 1938), where $f(\rho)$ is a quadratic polynomial function of $\rho$. The underlying assumptions of the model are as follows. Given a constant sex ratio, the number of meetings between males and females is proportional to $\rho^{2}$. The ratio of births to meetings can be affected by the population density and is hence assumed to be linearly decreasing in $\rho$. Also, the influx and efflux of individuals are represented by birth and death events that occur at constant per capita rates. Hence, the Volterra model takes the following form:

$$\dfrac{d\rho(t)}{dt}=-a_{1}\rho+(a_{2}-a_{3}\rho)\rho^{2}=-a_{1}\rho+a_{2}\rho^{2}-a_{3}\rho^{3},$$

(2.2)

where $a_{2},a_{3}>0$. The sign of $a_{1},$ however, is determined by the magnitude of the difference between the birth and death rates. If we define the two real-valued roots,

	$\displaystyle k_{1}$	$\displaystyle=\dfrac{1}{2a_{3}}[a_{2}-\sqrt{a_{2}^{2}-4a_{1}a_{3}}],$		(2.3a)
	$\displaystyle k_{2}$	$\displaystyle=\dfrac{1}{2a_{3}}[a_{2}+\sqrt{a_{2}^{2}-4a_{1}a_{3}}],$		(2.3b)

with $a_{2}^{2}>4a_{1}a_{3}$, then the model is often shown in the following form (Courchamp et al., 1999)

$$\dfrac{d\rho(t)}{dt}=a_{1}\rho\left(1-\dfrac{\rho}{k_{2}}\right)\left(\dfrac{\rho}{k_{1}}-1\right),$$

(2.4)

resembling the logistic model with the addition of a new unstable steady state, $\rho=k_{1}.$

If $a_{1}>0$, the Volterra model has three steady states: two stable states at $\rho=0$ and $\rho=k_{2}$ and an unstable state at $\rho=k_{1}$. In this case, if the initial population size exceeds $k_{1},$ the population grows over time and converges to the stable steady state $\rho=k_{2}$, the carrying capacity of the system. If the initial population size is less than $k_{1}$, the population decays over time to the stable steady state $\rho=0$ and goes extinct. This scenario describes the strong Allee effect. The dynamics when $a_{1}<0$ are different. In this scenario, $k_{1}$ becomes negative and the steady state $\rho=k_{1}$ lacks biological meaning. Here, the Volterra model has only two steady states, an unstable state at $\rho=0$ and a stable state at $\rho=k_{2}$. This case corresponds to the weak Allee effect.

In order to investigate the effects of internal fluctuations or demographic stochasticity in a two-population system with both weak and strong Allee effects, we consider the temporal evolution of the system as specified by a Markovian birth-death process. Demographic stochasticity is included both in the dynamics of the locations and in the dispersal parameter. We account for the individual reaction kinetics explicitly, in a mechanistic manner, without relying on phenomenological considerations. Here, we follow the approach of Méndez and colleagues (Méndez et al., 2019) by casting the system as a chemical reaction network that results in an individual-based model (IBM). The minimal IBM that displays both the weak and strong Allee effect and also accounts for dispersal can be described as follows. It consists of two birth processes (linear and binary birth), a ternary competition process, a linear death process, and an exchange process. We provide our reaction scheme below.

		$\displaystyle X_{1}\overset{\mu_{1}}{\rightharpoonup}(1+b)X_{1}$		(2.5a)
	$\displaystyle 2$	$\displaystyle X_{1}\overset{\lambda_{1}}{\rightharpoonup}(2+a)X_{1}$		(2.5b)
		$\displaystyle X_{1}\overset{\gamma_{1}}{\rightharpoonup}\emptyset$		(2.5c)
	$\displaystyle 3$	$\displaystyle X_{1}\overset{\tau_{1}}{\rightharpoonup}(3-c)X_{1}$		(2.5d)
		$\displaystyle X_{1}\overset{d}{\rightharpoonup}X_{2}$		(2.5e)
		$\displaystyle X_{2}\overset{\mu_{2}}{\rightharpoonup}(1+b)X_{2}$		(2.5f)
	$\displaystyle 2$	$\displaystyle X_{2}\overset{\lambda_{2}}{\rightharpoonup}(2+a)X_{2}$		(2.5g)
		$\displaystyle X_{2}\overset{\gamma_{2}}{\rightharpoonup}\emptyset$		(2.5h)
	$\displaystyle 3$	$\displaystyle X_{2}\overset{\tau_{2}}{\rightharpoonup}(3-c)X_{2}$		(2.5i)
		$\displaystyle X_{2}\overset{d}{\rightharpoonup}X_{1}$		(2.5j)

The first reaction is a linear birth process, which occurs at a constant rate $\mu_{1}$, and describes the baseline reproductive success of the first population in the absence of cooperative effects. It accounts for the fact that the typical individual produces $b$ offspring that reach reproductive age. The second reaction is a binary process that occurs at a constant rate $\lambda_{1}$. It describes cooperative interactions, such as breeding, antipredator behavior, or foraging, that result in producing $a$ additional offspring which reach reproductive age. The third reaction is a linear death process, occurring at constant rate $\gamma_{1}$, which accounts for mortality due to natural causes. The fourth reaction is a ternary competition process, accounting for the results of overcrowding and resource depletion, where $c$ individuals die at rate $\tau_{1}$. Note that $c=1,2,3$ are the only meaningful values. The next reaction is an exchange process of symmetric dispersal between the two populations. This occurs at a constant rate of $d$. The last five reactions in the scheme describe the dynamics of the second population, respectively.

The reaction scheme (2.5) defines a Markovian process, and the temporal evolution of $P(n_{1},n_{2},t)$, the probability of having $n_{i}$ individuals from the $i^{th}$ population at time $t$ for $i=1,2$, is described by the following master equation, also known as the forward Kolmogorov equation (Gardiner, 2004):

$$\dfrac{dP(\bm{n},t)}{dt}=\sum\limits_{\bm{r}}[W(\bm{n}-\bm{r},\bm{r})P(\bm{n}-\bm{r},t)-W(\bm{n},\bm{r})P(\bm{n},t)],$$

(2.6)

where $P(\bm{n}<\bm{0},t)=P(n_{1}<0,n_{2}<0,t)=0.$ Here $W(\bm{n},\bm{r})$ are the transition rates between the states with $\bm{n}$ and $\bm{n}+\bm{r}$ individuals, where

$$\bm{r}=\{\bm{r_{1}},\bm{r_{2}},\ldots,\bm{r_{10}}\}=\{(b,0),(a,0),(-1,0),(c,0),(-1,0),(0,b),(0,a),(0,-1),(0,c),(0,-1)\}$$

is the vector of transition increments corresponding to the system given by Eq. (2.5). The transition rates corresponding to each reaction, $W(\bm{n},\bm{r})$, are obtained from the reaction kinetics (Van Kampen (1992); Gardiner (2004)):

	$\displaystyle W(\bm{n},\bm{r_{1}})$	$\displaystyle=\mu_{1}n_{1},$		(2.7a)
	$\displaystyle W(\bm{n},\bm{r_{2}})$	$\displaystyle=\dfrac{\lambda_{1}}{2}n_{1}(n_{1}-1),$		(2.7b)
	$\displaystyle W(\bm{n},\bm{r_{3}})$	$\displaystyle=\gamma_{1}n_{1},$		(2.7c)
	$\displaystyle W(\bm{n},\bm{r_{4}})$	$\displaystyle=\dfrac{\tau_{1}}{6}n_{1}(n_{1}-1)(n_{1}-2),$		(2.7d)
	$\displaystyle W(\bm{n},\bm{r_{5}})$	$\displaystyle=dn_{2},$		(2.7e)
	$\displaystyle W(\bm{n},\bm{r_{6}})$	$\displaystyle=\mu_{2}n_{2},$		(2.7f)
	$\displaystyle W(\bm{n},\bm{r_{7}})$	$\displaystyle=\dfrac{\lambda_{2}}{2}n_{2}(n_{2}-1),$		(2.7g)
	$\displaystyle W(\bm{n},\bm{r_{8}})$	$\displaystyle=\gamma_{2}n_{2},$		(2.7h)
	$\displaystyle W(\bm{n},\bm{r_{9}})$	$\displaystyle=\dfrac{\tau_{2}}{6}n_{2}(n_{2}-1)(n_{2}-2),$		(2.7i)
	$\displaystyle W(\bm{n},\bm{r_{10}})$	$\displaystyle=dn_{1}.$		(2.7j)

Deterministic ODEs for the average population size can be obtained from Eq. (2.6). Multiplying Eq. (2.6) by $n_{1}n_{2}$, using transition rates (2.7), and summing over all values of $n_{1}$ and $n_{2}$, we find

$$\dfrac{d\rho_{i}}{dt}=(\mu_{i}b-\gamma_{i}-d)\rho_{i}+\dfrac{a\gamma_{i}}{2}\rho_{i}^{2}-\dfrac{c\tau_{i}}{6}\rho_{i}^{3},$$

(2.8)

for $i=1,2$, where $\rho_{i}=\langle n_{i}\rangle$ is the mean number of individuals in population $i$. We note that this deterministic equation holds strictly when the demographic fluctuations vanish, which occurs in the thermodynamic limit as the population sizes increase to infinity. Writing $\bm{k_{1}}=(k_{11},k_{12})$ and $\bm{k_{2}}=(k_{21},k_{22})$, Eq. (2.8) can be cast in the form of Eq. (2.4) with the definitions

	$\displaystyle\bm{k_{1}}$	$\displaystyle=\dfrac{3}{2c\tau_{i}}[a\lambda_{i}-\sqrt{a^{2}\lambda_{i}^{2}+8c\tau_{i}(b\mu_{i}-\gamma_{i}-d)/3}],$
	$\displaystyle\bm{k_{2}}$	$\displaystyle=\dfrac{3}{2c\tau_{i}}[a\lambda_{i}+\sqrt{a^{2}\lambda_{i}^{2}+8c\tau_{i}(b\mu_{i}-\gamma_{i}-d)/3}]$

At the deterministic level, the enzymatic reaction scheme (2.5) of the IBM gives rise to the Volterra rate equation Eq. (2.4).

For the sake of simplicity, in what follows, we focus on the simplest version of this IBM. Namely, we treat the Markovian process as a single-step process with $a=b=c=1.$ The set of reactions (2.5) then becomes

		$\displaystyle X_{1}\overset{\mu_{1}}{\rightharpoonup}2X_{1}$		(2.9a)
	$\displaystyle 2$	$\displaystyle X_{1}\overset{\lambda_{1}}{\rightharpoonup}3X_{1}$		(2.9b)
		$\displaystyle X_{1}\overset{\gamma_{1}}{\rightharpoonup}\emptyset$		(2.9c)
	$\displaystyle 3$	$\displaystyle X_{1}\overset{\tau_{1}}{\rightharpoonup}2X_{1}$		(2.9d)
		$\displaystyle X_{1}\overset{d}{\rightharpoonup}X_{2}$		(2.9e)
		$\displaystyle X_{2}\overset{\mu_{2}}{\rightharpoonup}2X_{2}$		(2.9f)
	$\displaystyle 2$	$\displaystyle X_{2}\overset{\lambda_{2}}{\rightharpoonup}3X_{2}$		(2.9g)
		$\displaystyle X_{2}\overset{\gamma_{2}}{\rightharpoonup}\emptyset$		(2.9h)
	$\displaystyle 3$	$\displaystyle X_{2}\overset{\tau_{2}}{\rightharpoonup}2X_{2}$		(2.9i)
		$\displaystyle X_{2}\overset{d}{\rightharpoonup}X_{1}$		(2.9j)

The mean-field rate equation corresponding to (2.9) is

$$\dfrac{d\rho_{i}}{dt}=(\mu_{i}-\gamma_{i}-d)\rho_{i}+\dfrac{\gamma_{i}}{2}\rho_{i}^{2}-\dfrac{\tau_{i}}{6}\rho_{i}^{3},$$

(2.10)

for $i=1,2$.

For this set of reactions, the master equation can be explicitly obtained as

	$\displaystyle\dfrac{dP(n_{1},n_{2},t)}{dt}$	$\displaystyle=(n_{1}-1)\left[\dfrac{\lambda_{1}}{2}(n_{1}-2)+\mu_{1}\right]P(n_{1}-1,n_{2},t)$		(2.11)
		$\displaystyle+(n_{1}+1)\left[\dfrac{\tau_{1}}{6}n_{1}(n_{1}-1)+\gamma_{1}+d\right]P(n_{1}+1,n_{2},t)$
		$\displaystyle+(n_{2}-1)\left[\dfrac{\lambda_{2}}{2}(n_{2}-2)+\mu_{2}\right]P(n_{1},n_{2}-1,t)$
		$\displaystyle+(n_{2}+1)\left[\dfrac{\tau_{2}}{6}n_{2}(n_{2}-1)+\gamma_{2}+d\right]P(n_{1},n_{2}+1,t)$
		$\displaystyle-n_{1}\left[\dfrac{\lambda_{1}}{2}(n_{1}-1)+\mu_{1}+\dfrac{\tau_{1}}{6}(n_{1}-1)(n_{1}-2)+\gamma_{1}+d\right]P(n_{1},n_{2},t)$
		$\displaystyle-n_{2}\left[\dfrac{\lambda_{2}}{2}(n_{2}-1)+\mu_{2}+\dfrac{\tau_{2}}{6}(n_{2}-1)(n_{2}-2)+\gamma_{2}+d\right]P(n_{1},n_{2},t)$

We note that the master equation (2.11) includes only single-step processes where the transitions take place between the states $\bm{n}$ and $\bm{n}\pm(1,0)$ or $\bm{n}$ and $\bm{n}\pm(0,1)$. We can define new dimensionless quantities in terms of the reaction rates as follows:

$$N_{i}=\dfrac{3\lambda_{i}}{2\tau_{i}},\ \ \ \delta_{i}^{2}=1+\dfrac{8\tau_{i}(\mu_{i}-\gamma_{i}-d)}{3\lambda_{i}^{2}},\ \ \ \ R_{0}^{(i)}=\dfrac{\mu_{i}}{\gamma_{i}+d},$$

(2.12)

Note that $N_{i}$ defines the scale of the typical size of population $i$ prior to extinction. The identities (2.12) establish a relation between the microscopic $(\lambda_{i},\mu_{i},\gamma_{i},\tau_{i},d)$ and macroscopic $(N_{i},\delta_{i},R_{0}^{(i)})$ parameters, which are obtainable through field observations. We now realize that the IBM displays both types of Allee effects:

	$\displaystyle\text{Weak Allee}:\ \ \ \mu_{i}>\gamma_{i}+d\ \ \text{or}\ \ R_{0}^{(i)}>1,$		(2.13a)
	$\displaystyle\text{Strong Allee}:\ \ \ \mu_{i}<\gamma_{i}+d\ \ \text{or}\ \ R_{0}^{(i)}<1.$		(2.13b)

Note that for $R_{0}^{(i)}>1\ (R_{0}^{(i)}<1)$ we have $\delta_{i}>1\ (\delta_{i}<1)$, where for the strong Allee effect, we must also demand that $\delta_{i}>0.$

As mentioned in the Introduction, our starting point for the model used in this work is the deterministic skeleton for the non-dimensionalized model described by Johnson and Hastings, reproduced here for convenience.

	$\displaystyle\dot{X_{1}}=X_{1}(\beta_{1}-(X_{1}-1)^{2})+D(X_{2}-X_{1})$		(2.14a)
	$\displaystyle\dot{X_{2}}=X_{2}(\beta_{2}-(X_{2}-1)^{2})+D(X_{1}-X_{2})$		(2.14b)

In the model above, the parameter $\beta_{i}$ represents a measure for the quality of the environment by the population denoted by $X_{i}$. The parameter $D$ denotes passive diffusion in the system and is an indicator of network connectivity. See (Johnson and Hastings, 2018) for a detailed exposition of the model.

In order to make our subsequent analyses feasible, we aim to reduce the dimensionality of the parameter space for our stochastic model. So, for the sake of illustration, we treat the stochastic rate parameters in (2.9) as identical for both populations. This assumption appears to be reasonable since a compelling justification for heterogeneity in the population parameters is unlikely. Moreover, this assumption is consistent with the model parameterization in the work by Johnson and Hastings.

Thus, our matching scheme can be written as follows:

	$\displaystyle\mu_{1}$	$\displaystyle:=\beta_{1},$		(2.15a)
	$\displaystyle\mu_{2}$	$\displaystyle:=\beta_{2},$		(2.15b)
	$\displaystyle\lambda$	$\displaystyle:=\lambda_{1}=\lambda_{2}\equiv 4,$		(2.15c)
	$\displaystyle\gamma$	$\displaystyle:=\gamma_{1}=\gamma_{2}\equiv 1,$		(2.15d)
	$\displaystyle\tau$	$\displaystyle:=\tau_{1}=\tau_{2}\equiv 6,$		(2.15e)
	$\displaystyle d$	$\displaystyle:=D,$		(2.15f)
	$\displaystyle\tilde{N}$	$\displaystyle:=\tilde{N}_{1}=\tilde{N}_{2}$		(2.15g)

So for $i=1$ and $2$,

$$\tilde{N}=1,\ \ \ \delta_{i}^{2}=\beta_{i}-D,\ \ \ \ R_{0}^{(i)}=\dfrac{\beta_{i}}{D+1}.$$

(2.16)

In what follows, we are interested in the case of bistability, which is manifest in the case of the strong Allee effect. Note that this necessitates the following condition:

$$D<\beta_{i}<D+1$$

(2.17)

We begin by noting that our stochastic model operates over a two-dimensional state space. We argue that the dimension of the state space is as low as possible but nevertheless captures the desired phenomena. Due to the presence of Allee effects, there are cubic nonlinearities in the system. So, the existence of an analytic solution to the system is highly unlikely. To circumvent this issue, Johnson and Hastings conducted numerical simulations using the deterministic skeleton in order to produce bifurcation diagrams. They focused on how the transition rate between the two patches, or $D$, determine the bifurcation structure for this system. In this study, however, we restrict our attention to the case of bistability, which is properly addressed with a stochastic model.

As discussed in the previous section, we chose to adapt a master equation approach for our model (Méndez et al., 2019). Instead of specifying a model with carrying capacities, we used an individual-based modeling approach using a chemical reaction network. This allows for a fine-grained representation of the underlying discrete, stochastic process. Using this approach, we wrote down the two-dimensional chemical master equation for the process.

Given that this stochastic process is a continuous-time Markov chain (CTMC), it can be explicitly described by a generator $\mathcal{Q}$-matrix with a countable state space. A nice feature of most ecological models is that they are built around processes that will approach a compact set exponentially fast. Any reasonable ecological model should not have unbounded population growth. Density dependence in ecology models is typically what gives this. This implies that the model effectively operates over a finite state space as the probability of arbitrarily large populations is negligibly small. Using this insight, we were able to obtain a finite-state CTMC in two dimensions.

Since the multi-dimensional master equation was relatively unwiedly to work with, we sought to reduce the two-dimensional state space to one dimension (Allen, 2010; Allen and Allen, 2003). Denoting $N:=N_{1}=N_{2}$ as the maximum number of individuals in either population, the specific mapping function used was

$$f\left((x,y)\right)=(N+1)x+y+1,\ \ \ \ \ 0\leq x,y\leq N.$$

(3.1)

We could also exploit the sparsity of the banded $\mathcal{Q}$-matrix to drastically simplify the effort needed in computations (van Doorn and Pollett, 2013). Specifically, the $\mathcal{Q}$-matrix has size $(N+1)^{2}\times(N+1)^{2}$ with $7N^{2}+4N-2$ nonzero entries, yielding a matrix density of $\mathcal{O}(N^{-2}).$ Thus, the sparsity of the matrix increases quadratically with $N.$

For the ensuing analyses, we obtained population thresholds as follows. The high thresholds were varied from $H:=H_{1}=H_{2}=2$ to $H=N$. Here, $N$ is defined as before, and can be understood as a system size parameter. The low thresholds were varied from $L:=L_{1}=L_{2}=\tilde{N}=1$ to $L=H-1$. This ensures that $L$ is always less than $H.$ Note that the smallest $L$ threshold, $L=1,$ is representative of quasi-extinction, or the typical size of either population before extinction.

In order to probe the system under consideration, we computed the mean first passage times (MFPTs) for the model with the state-space parameterized by $N$ for all combinations of the parameters $\beta_{i},D,L,$ and $H$ (Chou and D’Orsogna, 2014; Polizzi et al., 2016). This was done as follows. First, we formed the $\mathcal{Q}$-matrix for each point in parameter space. Then, the rows and columns corresponding to the trap states (e.g. extinction at $(x=0,y=0)$) were removed. We also formed a vector of initial state probabilities $p_{0}$ governing the subsequent evolution of the CTMC. This vector has $(N+1)^{2}$ entries. The corresponding entry was removed from this vector. Next, using the resulting truncated $\mathcal{Q}$-matrix, $\tilde{\mathcal{Q}}$, we computed the matrix of mean residence times, or $-\tilde{\mathcal{Q}}^{-1}$. Finally, we computed the sum of the entries in $-\tilde{\mathcal{Q}}^{-1}p_{0}$ to yield the MFPT from the initial state to the desired trap state.

We could then construct a compartmental system with a reduced state space consisting of just four states: $HH$, $HL$, $LL$, and $LH$. Each of these states corresponds to a combination of population thresholds. For instance, $HL$ means that the first population is at a high abundance and the second population is at a low level. The MFPTs from the original model were used as input rates for the transition rate $\mathcal{S}$-matrix of the reduced model. Thus we used an emulator, or meta-model, as a proxy to analyze the original system.

Throughout this work, we analyze the simplest possible system capable of exhibiting noise-induced tipping point cascades in the vicinity of saddle-node bifurcations. We can use the emulator described previously to construct the schematic diagram in Fig. 1. In the diagram, each $r_{i}$ for $i=1,\ldots,8$ represents the rate of the transition between the relevant compartments. Each rate can be computed as the inverse of its corresponding mean first passage time. Given that the process is represented as a CTMC, we can write down the transition rate matrix $\mathcal{S}$ for the emulator. The $n$-step transition probability matrix $(\mathcal{P}^{n})_{ij}$ gives the $n$-step transition probabilities of the system, from the state corresponding to the $i^{th}$ row to the state corresponding to the $j^{th}$ column.

The $\mathcal{S}$-matrix is given as:

$$\mathcal{S}=\begin{pmatrix}-r_{1}-r_{5}&r_{1}&r_{5}&0\\ r_{8}&-r_{2}-r_{8}&0&r_{2}\\ r_{4}&0&-r_{4}-r_{6}&r_{6}\\ 0&r_{7}&r_{3}&-r_{3}-r_{7}\end{pmatrix}$$

where the ordering of the states is $(HH,HL,LH,LL)$. Note that $\mathcal{S}$ satisfies the properties of a generator matrix, namely:

• •

  All off-diagonal elements are non-negative.

• •

  All diagonal elements are negative.

• •

  Each row sums to zero.

Since the two populations $X_{1}$ and $X_{2}$ interact with each other via the symmetric dispersal parameter $D$ and the thresholds for both populations are equivalent, we can treat the states $LH$ and $HL$ as equivalent. Hence, we can simplify our formalism by omitting the $HL$ state (Fig. 2):

$$\mathcal{\tilde{S}}=\begin{pmatrix}-r_{5}&r_{5}&0\\ r_{4}&-r_{4}-r_{6}&r_{6}\\ 0&r_{3}&-r_{3}\end{pmatrix}$$

where the ordering of the states is now $(HH,LH,LL)$.

Now, we obtain the following results:

	$\displaystyle P(HH\to LH)=1\ \text{with}\ MFPT=\frac{1}{r_{5}}$		(4.1)
	$\displaystyle P(LH\to HH)=\frac{r_{4}}{r_{4}+r_{6}}\ \text{with}\ MFPT=\frac{1}{r_{4}}$		(4.2)
	$\displaystyle P(LH\to LL)=\frac{r_{6}}{r_{4}+r_{6}}\ \text{with}\ MFPT=\frac{1}{r_{6}}$		(4.3)

By the multiplicative rule for probabilities, we obtain:

	$\displaystyle P(HH\to LL)=\frac{r_{6}}{r_{4}+r_{6}}\ \text{with}\ MFPT=\frac{1}{r_{5}}+\frac{1}{r_{6}}$		(4.4)
	$\displaystyle P(HH\to HH)=\frac{r_{4}}{r_{4}+r_{6}}\ \text{with}\ MFPT=\frac{1}{r_{5}}+\frac{1}{r_{4}}$		(4.5)

Hence, the probability of a system collapse differs from the probability of system resiliency in a fundamental manner that is quantified by the parameters

$\displaystyle r_{6}=\frac{1}{MFPT(LH\to LL)},\ \ r_{4}=\frac{1}{MFPT(LH\to HH)}$

(4.6)

If $r_{4}=r_{6}$, then $P(HH\to HH)=r_{4}/(r_{4}+r_{6})=r_{6}/(r_{4}+r_{6})=P(HH\to LL)=1/2$ so the system is equally likely to collapse or recover. We observe the following chain of implications in the case that $r_{4}>r_{6}$. The likelihood of recovery is larger than the threshold value of $\frac{1}{2}$, which implies that the MFPT decreases in the vicinity of the tipping point of the system. This in turn implies that $r_{4}>r_{6}$.

Our aim henceforth is to explore the multi-dimensional region of parameter space that corresponds to the (desired) probability of system resiliency. In symbols, for a fixed $N$, we want to identify all combinations $(D,\beta_{1},\beta_{2},H,L)$ where $D>0,\ \beta_{1}\in(D,D+1),\ \beta_{2}\in(D,D+1),$ $H\in[2,N],$ and $L\in[1,H-1]$ and the system recovers with the aid of rescue effects. In order to perform these numerical simulations, we began by specifying an upper bound of $1$ for the dispersal parameter $D.$ This value for $D$ could be considered as large, since $D$ should be commensurate with the per-capita rate of the system, which was chosen as $r=1$ in the non-dimensionalization of the model Eqs. (2.14a)-(2.14b) (Johnson and Hastings, 2018). This implies that $D\in(0,1).$ Then, we chose $N=10.$ This was the largest $N$ such that the largest condition number among the space of all matrices $\tilde{\mathcal{Q}}$, taken over the set of simulated parameter values, was less than $10^{7}$. For our chosen value of $N,$ we found this condition number to be $8.56\times 10^{6}.$ We note that numerical linear algebra was used here instead of Monte-Carlo methods, including the Doob-Gillespie algorithm (Gillespie, 1976, 1977)). The primary justification for this choice of method was that some of the simulated mean first passage times were very large, causing prohibitively long runtimes with Monte-Carlo simulations. Thus, $\beta_{1}$ and $\beta_{2}$ were restricted accordingly between $D$ and $D+1$. Finally, the high and low thresholds were varied with unit spacing. All computations were performed in MATLAB (MATLAB, 2020).

Defining $r=P(LH\to HH)=\frac{r_{4}}{r_{4}+r_{6}},$ we can analyze the parameter space. For the sake of illustration, we chose three values for the habitability parameters, for each value of $D$ (i.e. $\beta_{i}=D+0.01,D+0.5,D+0.99$). Upon analysis of the simulated data set, we found that the minimum value of $r$ over the entire parameter space was approximately $0.275$, and the maximum value was $1$. This indicates that it is possible to guarantee system recovery for some set of parameter values, but there is no set of parameter values that ensures system failure. In other words, the system is never more than $72.5\%$ likely to fail. Next, we varied the threshold value $\eta$, such that $r>\eta$, for $\eta=0.9,0.95,0.99$. How weak can the environment be, or equivalently, which resource constraints will not obstruct our goal of achieving $r>\eta$? We can intuit that the percentage of parameter space for which the system recovers will decrease monotonically with the threshold $\eta$. Naturally, we would like to find a region of parameter space such that $\eta$ is large and the proportion of parameter space $\nu$ for which $r>\eta$ is large.

Referring to Fig. 3, several observations can be made. For instance, low dispersal ($D=0.01$) guarantees that $\nu=0$ for $\eta\geq 0.9.$ This means that the desired probability of system resiliency is never $90\%$ or higher if the connectivity between the two patches is low. This makes sense, because the aid of rescue effects in such a system is muted, due to a poor degree of communication between the two populations. Also, $\nu$ decreases as $\eta$ is increased. This is plausible since $\eta$ behaves as a constraint on the parameter space and increasing $\eta$ implies greater specificity. In all cases, if both component Allee effects are weakly strong (i.e. $(\beta_{1},\beta_{2})=(D+0.99,D+0.99)$), then $\nu$ is maximized. This is especially noticeable when $D=0.99$, regardless of the value of the threshold parameter $\eta$. Lastly, if each heatmap is treated as a matrix, we see that the column (row) of a given matrix is increasing in $\nu$, as the corresponding column (row) indices are increased. This is also reflected within the coloring scheme of each matrix. As expected, $\nu$ is maximal with a value of $0.78$, for $D=0.99,\eta=0.9,\beta_{1}=D+0.99,$ and $\beta_{2}=D+0.99.$

In Fig. 4, we have made explicit the correspondence between Fig. 3 and the combination of high/low thresholds. In concert with the analysis in Fig. 3, we focused here on the case of high dispersal (i.e. $D=0.99$). We see that the combination $(H,L)=(2,1)$ is systematically undesirable for the system, as it always results in the lowest value of $r$ for a given pair $(\beta_{1},\beta_{2})$. However, as the high threshold is increased while the low threshold is kept at the same level, we see that $r$ increases consistently. This is true for any value of the low threshold (i.e. with only $H$ increased). In general, a moderate level of separation between the high and low thresholds is needed to guarantee a desirable outcome for the system. In other words, in order to tip positive change with full confidence ($r\approx 1$), the difference $H-L$ should be sufficiently large. This is sensible in the context of rescue effects, since the high and low thresholds are effectively separated, resulting in a resilience-averaging effect for the system. If the two thresholds were commensurate, the benefit to either population would be reduced due to a lower margin of error between a potential catastrophe or an eventual success. This can be seen clearly in the top-most diagonal of the matrix for each case, as the values of $r$ are low compared to their neighbors, as expected.

A key issue which relates to the idea of a rescue effect is the relative likelihood of the transition from $HH$ to $LH$, as compared with the transitions $LH$ to $LL$ and $LH$ to $HH$. Referring to Eqs. (4.1) - (4.3), we can reason as follows. The likelihood of system resiliency in the form of a rescue effect is given by $P(LH\to HH)=r$, whereas the probability of a total failure or catastrophic collapse is given by $P(LH\to LL)=1-r$. Note that (4.1) indicates that a partial failure of the system is inevitable (i.e. with probability $1=r+(1-r)$). This result follows from the nature of the one-step transition probabilities in the emulator framework. Thus, we can reason that the odds of system recovery is equal to $\frac{r}{1-r}$. As illustrated in Fig. 5, we see that all of the conclusions observed in Fig. 4 are consistent here as well. For instance, the system is extremely likely to tip favorably in the presence of high dispersal, if both Allee effects are weakly strong. In symbols, $D=0.99,H=10,L=1$, and $\beta_{1}=\beta_{2}=D+0.99$. The odds of tipping in this case are nearly $50,000.$

Given this series of observations, there is an interesting and intuitive explanation that generalizes the given setting. Noting that $0<r<1$ since $r$ is a probability, we have that

$$\dfrac{r}{1-r}=r+r^{2}+r^{3}+\ldots$$

is a convergent geometric series. Thus, we can write the left-hand side of the equation above as the odds of one network tipping favorably, where the network consists of a system with two patches. The right-hand side, however, can be interpreted as the cumulative probability of all such networks tipping favorably, viz.

$$\dfrac{r}{1-r}=\lim\limits_{N\to\infty}\sum\limits_{n=1}^{N}r^{n},$$

where $n$ is the number of two-patch systems in a network. With this analysis, we thus see an instance of a tipping cascade of networks in the form of a domino effect. Note that the propagation of domino dynamics through a network requires strong connectivity, or high dispersal, in our model (Lenton, 2020).

Our analysis of a minimal, stochastic system that could have cascading tipping points yields some interesting ecological and mathematical insights. A novel result of this work is that noise-induced tipping can distinguish a catastrophic collapse and a successful recovery from the brink of extinction. This idea is expressed succinctly as the distinction between the tipping of positive change and the failure to adequately address an impending critical transition (Lenton, 2020). In our two-population model, we showed that if the system is in the low/high state, a stochastic perturbation due to demographic noise results in either a collapse of both populations (low/low state) or a full recovery to the high/high stable state. This may be a key feature of spatially connected populations with strong Allee effects, and analyzing these dynamics should inform management strategies for similar ecological systems. We find that the population with higher resilience, in the form of a larger Allee threshold, has the capacity to save its counterpart through a rescue effect. This occurs in particular due to the presence of dispersal between the populations. In the presence of noise, patch homogeneity, and strong network connectivity, the system is more likely to exhibit a tipping point cascade (Lenton, 2020).

Some discussion of systems with noise-induced tipping and their features is needed. In situations where tipping one system increases the probability of another system tipping - for example, melting of the Greenland ice sheet increases the likelihood of failure of the Atlantic Meridional Overturning Circulation (AMOC) - the first system should act as a proxy for the entire network, in terms of a call to action. This is an important point, because by definition, our system can tip without warning, thus precluding detection through generic early-warning signals. However, we have been able to characterize the likelihood of its propensity to tip, either favorably or catastrophically (Lenton, 2020). The results from our analysis support the hypothesis that both stochastic and transient dynamics might play a key role in regime shifts and their associated tipping (Rocha et al., 2018; Hastings et al., 2018). In particular, domino effects have relatively slow temporal dynamics and larger spatial scales; in our case, the dispersal parameter dominates the time scale in magnitude. In addition, the one-step transition probabilities in our model support the notion that structural dependencies manifest as one-way interactions for the domino effect. This is a salient feature of regime shift couplings. In our system, we note that a partial collapse is inevitable. Most examples of regime shifts exhibiting domino effects concern the Earth system, including monsoon weakening and thermohaline circulation collapse, as well as nutrient transport mechanisms (Scheffer, 2009; Rocha et al., 2018).

Our analysis is based on an individual-based model that describes the dynamics of the system under consideration and this approach should have general applicability. The mean-field description of this model necessarily includes a quadratic term. In other work (Abraham, 1991; Klose et al., 2020), the proposed tipping element omits the quadratic term, in an attempt to describe a “dangerous” bifurcation in the form of a cusp catastrophe. This would correspond to zero mortality due to natural causes in our model (i.e. $\gamma=0$). Similarly, the absence of the quadratic term would indicate that both the carrying capacity and Allee threshold for either population are identically zero (Johnson and Hastings, 2018). Thus, we see that important details of the underlying mechanisms of a tipping element can be lost in a phenomenological description.

Ways to generalize our work include looking at more than two populations which would allow the study of how specific forms of connectivity could play a role in tipping cascades, or introducing the possibility of more complex dynamics (Strogatz, 2001). Allowing for patch heterogeneity in the two populations by distinguishing their stochastic reaction rates could also produce interesting dynamics. We leave the exploration of these aspects to future work. Given the ubiquity of Allee effects (Courchamp et al., 2008) and the generality of our model, we believe that we have uncovered important ecological conclusions that are robust and should apply in a variety of settings. Within the field of landscape ecology, rescue effects can be crucial in both metapopulations and species augmentation efforts. Other potential areas of application include communication systems in network theory and regulatory networks in systems biology.

The authors declare that they have no conflicts of interest.