Dynamic environment-induced multistability and critical transition in a metacommunity ecosystem

We study a metacommunity model of consumer-resource populations coupled via dispersal under an environment-dependent framework, and we explore the occurrence of multistability and critical transition. By emphasizing two magnitudes acting on a dynamic environment at temporal and spatial scales, the coupled system with simple diffusive coupling and the nonlinear environmental coupling enables various interesting complex dynamics such as bistability, multistability, and critical transitions. Using the basin stability measure, we find the probability of attaining each alternative state in a multistable region. In addition, critical transitions (one from a high to a low species density and the other from a low to a high species density) are identified at different magnitudes in the presence of stochastic fluctuations. We also explore the robustness of critical slowing-down indicators, e.g., lag-1 autocorrelation and variance, to forewarn the critical transition in the metacommunity model. Further, a network structure also identifies synchronization and multiclustering for a different choice of initial conditions. In contrast with the earlier studies on dynamic environmental coupling, our results based on the defined magnitudes provide important insights into environmental heterogeneity, which determines the set of environmental conditions to predict metacommunity stability and persistence.

Figure 1

Schematic representation of the coupled system (1) with two patches in a dynamic environment. Here ${\alpha }_i$ denotes the consumption (predation) rate of the consumer species in the $i\mathrm{th}$ patch over its local resource $V$, whereas ${\epsilon }_i$ denotes the consumption over the environmental resource $E$.

Figure 2

One- and two-parameter bifurcation diagrams depicting multistability: (a) one-parameter bifurcation diagram of the resource $\left(V_i\right)$ at different fraction $F_{\epsilon \alpha }$ with $\alpha =0.4$ and $d=0.6$, (b) one-parameter bifurcation diagram of the consumer $\left(H_i\right)$ at different fraction $F_{\epsilon \alpha }$ with $\alpha =0.4$ and $d=0.6$, (c) $\epsilon -\alpha$ space as a two-parameter bifurcation diagram, and (d) $\epsilon -d$ space as a two-parameter bifurcation diagram. The dotted lines denote a particular case representing (a) and (b). Here the consumption rates in both of the patches are identical (${\epsilon }_1={\epsilon }_2=\epsilon$). In (a) and (b), green and blue circles denote stable and unstable limit cycles, respectively, whereas solid red and dashed black curves denote stable and unstable steady states, respectively. Here HB, TB, BP, SN, and LP denote Hopf bifurcation, transcritical bifurcation, branch point, saddle-node, and limit point, respectively. BS and MS denote the region of bistability and multistability, respectively. Other parameters are $r=0.55,K=0.55,{\alpha }_i=0.4,B=0.16,\beta =0.75,m=0.2,\gamma =0.6,r_1=0.5,K_1=0.6$, and $C=0.6$.

Figure 3

Time series and the corresponding basin stability measures: In the left panel, each time series is set with an index, whereas the right panel denotes the percentage of corresponding dynamics from ${10}^4$ simulations with random initial conditions. Multistability (a) MS1 for $\epsilon =0.25$ and $\alpha =0.4$, (b) MS2 for $\epsilon =0.45$ and $\alpha =0.4$, (c) MS4 for $\epsilon =0.5$ and $\alpha =0.4$, and (d) MS3 for $\epsilon =0.55$ and $\alpha =0.45$. The stable oscillations shown here are synchronized in both the interacting patches. Other parameters are $d=0.6,r=0.55,K=0.55,B=0.16,\beta =0.75,m=0.2,\gamma =0.6,r_1=0.5,K_1=0.6$, and $C=0.6$.

Figure 4

(a) One-parameter bifurcation diagram for varying the fraction of consumption in the environmental resource ($F_{{\epsilon }_1{\epsilon }_2}$) for $d=0.4$, and (b) one-parameter bifurcation diagram for varying the fraction of consumption in the environmental resource ($F_{{\epsilon }_1{\epsilon }_2}$) for $d=0.8$. Other parameters are $r=0.55,K=0.5,{\alpha }_i=0.4,B=0.16,\beta =0.75,m=0.2,\gamma =0.6,r_1=0.5,K_1=0.6$, and $C=0.6$.

Figure 5

Two-parameter bifurcation diagram with variations in $F_{{\epsilon }_1{\epsilon }_2}$ and $d$: Here Eq and NEq denote equilibrium (steady states) and nonequilibrium dynamics (oscillations), respectively. The bistable and multistable regions are denoted by BS and MS, respectively. Other parameters are $r=0.55,K=0.5,{\alpha }_i=0.4,B=0.16,\beta =0.75,m=0.2,\gamma =0.6,r_1=0.5,K_1=0.6$, and $C=0.6$.

Figure 6

Critical transitions using multiplicative noise for $d=0.4$ when (a) $F_{{\epsilon }_1{\epsilon }_2}$ varies in the forward direction, and (b) $F_{{\epsilon }_1{\epsilon }_2}$ varies in the reverse direction. Critical transition for $d=0.8$ when (c) $F_{{\epsilon }_1{\epsilon }_2}$ varies in the forward direction, and (d) $F_{{\epsilon }_1{\epsilon }_2}$ varies in the reverse direction. Other parameters are $\sigma =0.009,r=0.55,K=0.5,{\alpha }_i=0.4,B=0.16,\beta =0.75,m=0.2,\gamma =0.6,r_1=0.5,K_1=0.6$, and $C=0.6$.

Figure 7

Early warning signals calculated from the simulated time series data of the stochastic model [Eq. (3)] in the case of forward [see Fig. 6] and backward [see Fig. 6] transitions with variations in $F_{{\epsilon }_1{\epsilon }_2}$. The lag-1 autocorrelation [AR(1)] and the variance are calculated using a moving window of half the length of the time series segments [segments are indicated by the boxed regions in Figs. 6 and 6]. (a) In the case of forward transition, neither AR(1) nor variance gives a reliable signal of an upcoming transition. (b) In the case of a backward transition, both indicators show a clear increasing trend and reliably indicate an upcoming critical transition. For more details, see the text.

Figure 8

Synchronization, multiclustering, and species extinction in a network of 16 patches: Spatial dynamics are represented in the top panel, and the respective temporal dynamics are represented in the bottom panel. (a),(d) Synchronized oscillations of 16 patches form 1-cluster; (b),(e) dynamics of stable steady states with three clusters in which one of the clusters reaches the critical zero state and species become extinct; and (c),(f) dynamics of stable steady states with multiple clusters that contain both nonzero and zero species density. The parameter values are $d=0.5,\epsilon =0.06,r=0.55,K=0.55,{\alpha }_i=0.4,B=0.16,\beta =0.6,m=0.15,\gamma =0.7,r_1=0.5,K_1=1.6$, and $C=0.8$.

Figure 9

The stochastic time series of the model (1), generated using the Gillespie algorithm, exhibits flickering between alternative states.