Stochastic sensitivity analysis of noise-induced transitions in a predator-prey model with environmental toxins

1.   Introduction

With the rapid development of modern industries, the problem of industrial pollution caused by human activities becomes more and more severe. In aquatic ecosystems, the accumulation of toxic compounds in living organisms has caused great harm to aquatic species and human health. To accurately assess the risk of pollutants, over the past several decades, many mathematical toxicity-extrapolation models have been formulated, including population models, ecosystem models and landscape models ^[2,3,4,5,6].

Toxin-dependent IBM (individual-based models) and matrix population models are widely used to evaluate the effects of toxic chemicals on individual organisms and population dynamics ^[4,7]. In Ref. ^[8], Huang et al. proposed a toxin-dependent aquatic population model by incorporating the effects of toxin on growth and recruitment, and the authors further obtained threshold values of toxin concentration in the environment to maintain populations and prevent extirpation. All above-mentioned models are single-species models where populations are assumed to take up toxin only from exposure to the environment. It is well known that the main source of toxins in high-level trophic organisms is through ingestion of food. The higher-level animals in a food chain consume and store more pollutants (normally called biomagnification) ^[9,10]. Taking into account the cumulative effects of toxic pollutants in organisms, Huang et al. ^[1] mechanistically extended the single-species toxin-dependent model ^[8] to a predator-prey model with toxin effect as follows:

where $x$, $y$ are the concentrations of the prey (small fish) and the predator (rainbow trout) with environmental toxins, respectively; $u$, $v$ are the concentrations of toxins in prey and rainbow trout, respectively; $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ are maximum reproduction rate, effect of toxins on the growth and crowding effect of prey, respectively; $k_{1}$, $k_{2}$ are effect of toxins on the mortality of prey and rainbow trout, respectively; $m_{1}$, $m_{2}$ are natural mortality rate of prey and rainbow trout, respectively; $\xi$ is per capita feeding rate; $\eta$ is half-saturation constant; $\beta_{1}$, $\beta_{2}$ are reproduction efficiency and effect of toxins on the reproduction of rainbow trout, respectively; $a_{1}$, $a_{2}$ are uptake coefficient for prey and rainbow trout, respectively; $\sigma_{1}$, $\sigma_{2}$ are depuration rate for prey and rainbow trout, respectively.

For convenience, the nondimensionalized model (1.1) was introduced by letting

Dropping the tildes for notational simplicity, the model (1.1) can be rewritten in its dimensionless form:

In reality, the dynamics of the depuration due to metabolism of toxin have a much faster time scale than the dynamics of population biomass growth. This means the uptake of toxin and depuration may approach a quasi-equilibrium state on a fast time scale ^[1]. Therefore $\epsilon = \alpha_{1}/\sigma_{1}$ can be defined by a small perturbation parameter. Model (1.2) is a high dimensional system, letting $\epsilon\rightarrow 0$ in (1.2), then

In order to ensure the persistence of prey and predator, we assume that $m_{1} < 1, T < \min\{1, \sigma_{2}/c\}, m_{2} < \beta_{1},$ substituting (1.3) into the first and second equations of (1.2), we obtain the following two-dimensional predator-prey model with environmental toxins:

with

In ^[1], the authors pointed out that high toxin concentration in the environment is harmful to both species, and may lead to extirpation of both species, while low toxin concentration produces counterintuitive results. What'more, the phenomena of bistability occur when appropriate parameter values are taken.

Undoubtedly, environmental noises can change the qualitative behaviors of a deterministic model. In order to better and more robustly understand the impact of toxins on the persistence of populations, we assume that the noises only depend on the predation term, then we extend the deterministic model (1.4) to incorporate the noises as the following stochastic differential equations:

where $B_{1}(t)$ and $B_{2}(t)$ are standard one-dimensional independent Brownian motions, $\delta_{1}$ and $\delta_{2}$ are the noise intensities. For the simplicity, in the remaining parts of this paper, we assume $\delta_{1} = \delta_{2} = \delta$.

Our main goal of this paper is to study the phenomena of noise-induced transitions from coexistence to prey-only extirpation. In fact, many deterministic models with noise-induced transitions have been studied. The presence of coexisting attractors under random perturbations can generate new dynamic regimes, which have no analogues in the deterministic case ^[11,12]. In ^[13,14], noise-induced perturbations change the coexistence state to extinction for a chemostat model. Perturbations effects of limit cycles caused by noises were studied in ^[15,16]. Nonlinear dynamical models also show new phenomena, such as stochastic resonance ^[17,18], noise-induced transitions ^[19], noise-induced order ^[20,21], noise-induced chaos ^[22], and noise-induced complexity ^{[23,24,25,26,27,28,29]}. The analysis of the noise effects on dynamical systems with multiple stable states attracts attention of many researchers ^[22,30]. Multistable systems exhibit complex dynamics with noise-induced hopping between coexisting attractors and their basins of attraction ^{[31,32,33,34]}. The sensitivity analysis of random forced oscillations is important for investigating these transitions. For constructive analytical description of randomly forced equilibria and cycles of discrete-time systems, a new numerical method based on the technique of stochastic sensitivity functions (SSF) has been proposed in ^[35].

The rest of this paper is organized as follows. In Section 2, we analyze qualitative changes in stochastic dynamics using phase portraits. In Section 3, based on the technique of stochastic sensitivity functions and the method of confidence domains, we explore the phenomenon of noise-induced transitions from coexistence to prey-only extirpation, and construct the confidence ellipse and confidence band for model (1.5). Finally, we summarize the paper by a brief conclusion in Section 4.

2.   Attractors and stochastic phenomena in predator-prey model

In this paper, we take parameter values from ^[1]:

for model (1.5) with $\delta = 0$ and a varying parameter $\eta$. For any $\eta$, the model has a extinction equilibrium $E_{0, 1} = (0, 0)$, a prey-only equilibrium $E_{0, 2} = (1.6667, 0)$, and two coexistence equilibria $E_{1} = (0.2814\eta, 0.3244\eta+0.0913)$, $E_{2} = (0.8363\eta, 0.2491\eta+0.1135)$ when we restrict toxin level and the half-saturation constant such that both populations can coexist.

Fixing parameter $\eta = 1$, simple calculations show that the equilibrium $E_{0, 1}$ is unstable, $E_{0, 2}$ and $E_{1}$ are locally asymptotically stable, and $E_{2}$ is a saddle point. Vector field of the deterministic model is presented in Figure 1(a), here the red dashed-dotted line is the separatrix of two attraction basins. For deterministic model, trajectories starting from the inside of the separatrix tend to the coexistence equilibrium $E_{1}$, and trajectories starting from the outside of the separatrix tend to the prey-only extirpation equilibrium $E_{0, 2}$, that is, depending on the initial population levels, either the prey excludes the predator $E_{0, 2}$ or both species coexist at the interior equilibrium $E_{1}$.

Consider dynamics of model (1.4) under random disturbances. In Figure 2 and 3, time series of $x(t)$ and $y(t)$ (a) and random trajectories (b) are plotted. For small noise, the dynamics of the stochastic model is quite regular: a stochastic trajectory with an initial point near the deterministic coexistence equilibrium $E_{1}$ will leave this attractor, but stay in its small neighborhood, see Figure 2 for $\delta = 0.1$. However, as the noise intensity increases, the stochastic dynamics of model (1.5) becomes irregular: a stochastic trajectory exits the neighborhood of the coexistence equilibrium $E_{1}$ and crosses the separatrix of two attraction basins and finally tends to the prey-only equilibrium $E_{0, 2}$, see Figure 3 for $\delta = 0.2$. This phenomenon of noise-induced transition implies that the coexistence of two species can be destroyed by noise disturbances and the predator can be driven to extirpation.

For fixed $\eta = 0.5$, the coexistence equilibrium $E_{1}$ of the deterministic model (1.4) becomes unstable. Meanwhile, the model has a unique stable limit cycle surrounding the equilibrium $E_{1}$, where both population levels fluctuate around a coexistence equilibrium, see Figure 1(b), in which the red dash-dotted homoclinic orbit is the separatrix of two attraction basins. For the deterministic model, trajectories starting from the inside of the homoclinic orbit tend to the limit cycle, and trajectories starting from the outside of the homoclinic orbit tend to the prey-only equilibrium $E_{0, 2}$. Similar to the previous situation, in the stochastic case, the limit cycle generally disappears, but the trajectory stays in a small neighborhood of the deterministic limit cycle for small noise, see Figure 4 for $\delta = 0.007$. As the noise intensity increases, the stochastic trajectory exits the neighborhood of the limit cycle and tends to the prey-only equilibrium $E_{0, 2}$ (extirpation of the predator), see Figure 5 for $\delta = 0.012$.

3.   Analysis of noise-induced transitions via confidence domains

Noise-induced transitions between attraction basins of two equilibria or between attraction basins of a limit cycle and an equilibrium occur for $\delta > 0$. In this section, we want to use the knowledge of SSF technique to theoretically analyze these noise-induced transitions with the aid of confidence domains. We first construct the confidence ellipse for the stochastic model (1.5) with $\eta = 1$, and then estimate the threshold value of noise intensity for noise-induced transition from coexistence to predator extinction. Denote the deterministic coexistence equilibrium $E_{1} = (x_{1}, y_{1})$ and define

where

From $(A.2)$ in ^[36], we know that the stochastic sensitivity matrix

satisfies the following equations:

Moreover, we know that the confidence ellipse equation from $(A.3)$ in ^[36] is

where $\delta$ is the noise intensity and $P$ is a fiducial probability. From the above, we obtain

respectively. The confidence ellipse equation of $E_{1}$ is

For fixed fiducial probability $P = 0.95$, we change the noise intensity and observe how noise affects the confidence ellipse. For $\delta = 0.12$, random states of the stochastic model (1.5) and the corresponding confidence ellipse are plotted in Figure 6(a). We can observe that random states are distributed around the corresponding deterministic coexistence equilibrium, and they belong to the interior of the confidence ellipse with probability $0.95$.

Next, we take different noise intensities $\delta = 0.1$, $0.18$, and $0.3$, resulting in the corresponding confidence ellipses shown in Figure 6(b). Obviously, as the noise intensity increases, the confidence ellipse begins to expand and pass through the separatrix, and finally it enters the attraction basin of the predator extirpation equilibrium, that is, trajectories of the noised model with a high noise intensity can leave the attraction basin of the coexistence equilibrium $E_{1}$ and localize near the predator extirpation equilibrium. The noise intensity corresponding to the intersection of the confidence ellipse and the separatrix can be used as an estimation of the threshold value $\tilde{\delta}$. Here $\tilde{\delta}\approx0.18$. The transition escaping from the deterministic coexistence equilibrium $E_{1}$ to the predator extirpation equilibrium $E_{0, 2}$ has a typical example shown in Figure 3 for $\delta = 0.2 > 0.18$. In this case, the extirpation of the predator is mainly caused by large noise.

We now take $\eta = 0.5$. To find the configurational arrangement of the neighborhood of the deterministic limit cycle for small noise, we construct the confidence band for the stochastic model (1.5), and denote the deterministic limit cycle by $\Gamma(x(t), y(t)), t\in[0, T]$, where $T$ is the period of the periodic orbit. Let

and define matrices $F(t)$, $G(t)$ and $S(t)$ as follows:

where

From (A.4) in ^[36], we know that the stochastic sensitivity function $\mu(t)$ satisfies the following boundary problem:

where

Here,

are elements of a vector function $p(t) = (p_{1}(t), p_{2}(t))^{T}$ orthogonal to vector $(F_{1}(x, y), F_{2}(x, y))^{T}\mid_{\Gamma}$. It follows from (A.5) in ^[36] that the boundaries $\Gamma_{1, 2}(t)$ of the confidence band have the following explicit parametrical form:

where the parameter $k$ is connected with the fiducial probability $P$ by the formula $k = erf^{-1}(P)$, in which $erf(x) = \frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt$ is the error function.

In Figure 7(a), the confidence band around a stable cycle $\Gamma$ with $P = 0.95$ is plotted for $\delta = 0.01$. the rose line is the deterministic limit cycle, the blue points are the random states at different times and the two green dash-dotted lines are the boundaries of the confidence band. Obviously, the random states leave the deterministic limit cycle and form some distribution around it, and they belong to the interior of the confidence band with probability $0.95$.

Similar to confidence ellipse, we take noise intensities $\delta = 0.005, \delta = 0.011, \delta = 0.018$, resulting in the corresponding outer boundary (green lines) of confidence band shown in Figure 7(b). As the noise intensity increases, the outer boundary begins to expand and pass through the separatrix (red dash-dotted line), and finally it enters the attraction basin of the predator extirpation equilibrium, that is, trajectories of the noised model with a high noise intensity can leave the attraction basin of the limit cycle $\Gamma$ and localize near the predator extirpation equilibrium. The noise intensity corresponding to the intersection of the confidence band and the separatrix can be used as an estimation of the threshold value $\tilde{\delta}$. Here $\tilde{\delta}\approx0.011$. The transition escaping from the deterministic limit cycle $\Gamma$ to the predator extirpation equilibrium $E_{0, 2}$ has a typical example shown in Figure 7 for $\delta = 0.012 > 0.011$.

4.   Conclusion

In this paper, we study noise-induced transitions from the zone of coexistence to the zone of predator extirpation for the stochastically forced predator-prey model with environmental toxins. The corresponding deterministic model has two stable attractors ^[1]: the boundary prey-only attractor $E_{0, 2}$ corresponds to the extirpation of predator, and the internal attractor (either the stable equilibrium $E_{1}$ or the limit cycle $\Gamma$) corresponds to the coexistence of prey and predator. Two attraction basins are separated by the separatrix, which is the stable manifold of the saddle point $E_{2}$, that is, solutions starting from one attraction basin eventually approach the corresponding attractor. However, when there exists noise disturbance, the separatrix can be crossed by solutions due to large noise, and the solutions in the attraction basin of the coexistence equilibrium/limit cycle can eventually approach the predator extirpation equilibrium, see Figures 3 and 5.

In order to analyze the phenomenon of noise-induced transitions, the confidence ellipse of our stochastic model is constructed by using the technique of SSF. We provide the general location of the equilibrium and analyze the random states distribution of the stochastic model. For small noises, the confidence ellipse completely belongs to the attraction basin of the coexistence equilibrium, and the random trajectories will not leave the confidence ellipse with high probability. As the noise intensity increases, the confidence ellipse starts to expand. After it crosses the separatrix of two attraction domains, it will leave the attraction basin of the coexistence equilibrium. Consequently, random trajectories entering the confidence ellipse can leave the attraction basin of the coexistence equilibrium and approach the predator extirpation equilibrium with high probability. In addition, we establish the confidence band for the stochastic model (1.5), from which we provide the general location of the stochastic limit cycle and analyze the distribution of random states around the deterministic limit cycle.

We have expounded how the confidence domain method is used to understand the qualitative changes in stochastic dynamics and estimate a threshold value of noise intensity at which noise-induced extirpation can occur with high probability. This method is applicable to nonsmooth predator-prey models and more complex higher dimensional models in terrestrial ecosystems. Our results enrich the study of asymptotic behaviors in predator-prey models and help better understand the predator-prey dynamics in the stochastic perspective.

Acknowledgements

S. Yuan and D. Wu are partially supported by the National Natural Science Foundation of China (No. 11671260), Shanghai Leading Academic Discipline Project (No. XTKX2012), and Hujiang Foundation of China (B14005). H. Wang is partially supported by an NSERC discovery grant.

Conflict of interest

All authors declare no conflicts of interest in this paper.