Dynamic behavior on a multi-time scale eco-epidemic model with stochastic disturbances

1.   Introduction

The interaction between prey and predator has become one of the main relationships in ecology due to its universality and importance. Meanwhile, disease is one of the most common influencing factors on prey-predator systems in the natural environment. The study of disease dynamics in an ecological system is an important issue from both mathematical and ecological point of view, as spread of infectious diseases becomes an important factor to regulate animal population sizes. Many researchers have proposed and studied different predator–prey models in the presence of disease. Gupta and Dubey discussed bifurcation and chaos in a delayed eco-epidemic model ^[1]; dynamic behavior of an eco-epidemic model with two delays was analyzed by Jana et al.^[2]; Debasis Mukherjee discussed the Hopf oscillation phenomenon of an eco-epidemic model ^[3].

Meanwhile, as is well-known, the fluctuation in environmental changes is an important component of an ecosystem in natural ecological environments. Due to the stochastic perturbation in the real world, the dynamic behavior described by deterministic models cannot accurately predict the true future behavior of the system. Thus, considering the stochastic models with noise is reasonable. Related work could be found in Mukherjee ^[4]; Khare et al. ^[5]; and Zhang and Wu ^[6,7]. These indicate that many scholars have begun to pay attention to the impact of stochastic factors on ecological populations and have made many academic achievements.

On the other hand, another practical issue in the ecological environment is that the growth rate of prey is generally much faster than that of predator in most systems. Therefore, in most applications, the dynamics of different variables in simple differential equation systems are scaled hierarchically. Based on the stochastic system, many authors have considered system research and management at different timescales. Wang and Roberts discussed the dynamics of a slow-fast stochastic system by using random slow manifold reduction ^[8]. Stochastic P-bifurcation and the asymptotic behaviors of a stochastic model of Alzheimer's disease, which has two timescales, were analyzed by Zhang and Wang ^[9].

Random slow manifolds are geometrical invariant structures of multi-scale stochastic dynamical systems. It has been shown that the deterministic slow manifold has dramatic effect on the overall dynamical behavior. In ^[10], the authors discover that the bifurcation structure of the deterministic slow manifold creates a reaction channel for non-equilibrium transitions, leading to vastly increased transition rates. The random slow manifold also has a profound impact on the stochastic bifurcation for the stochastic dynamical system. It is beneficial to take advantage of random slow manifolds in order to examine dynamical behaviors of multi-scale stochastic systems, either through the slow manifolds themselves or by the reduced systems on these slow manifolds ^[11]. For multidimensional stochastic slow-fast systems, the reduction on random slow manifolds brings great convenience to the analysis of dynamic behavior.

Based on the above discussion, this paper is to discuss a stochastic mathematical model which has two timescales with disease infection in the prey population. As far as we know, the study of the epidemic model considering stochastic noise and multi-time scale is not extensive, especially the further study of the stochastic system. In our work, the further reduction analysis of multi-time scale stochastic system using random slow manifold can analyze its dynamic behavior, such as stochastic bifurcations, while maintaining accuracy. This makes the analysis of the system closer to reality, more biologically significant and more comprehensive. The P-bifurcation analysis after dimension reduction also provides a more specific and further discussion on the dynamic behavior on the dynamical behavior of the ecosystem affected by stochastic noise.

The rest of paper is divided as follows: In Section 2, a stochastic eco-epidemiological system with multi-time scale is presented. The stability and bifurcation of the equilibrium points of the model without noise are discussed in Section 3. In Section 4, the existence, boundedness, and ergodic property of solutions are analyzed on the stochastic model. In addition, the singular perturbation problem with random center-manifold method and the stochastic P-bifurcation are considered in Section 5. Finally, in Section 6, we end the paper with some concluding remarks.

2.   Model formulation

Zhang et al. ^[12] proposed a prey-diseased predator model with stochastic disturbances:

where $S(t), I(t),$ and $P(t)$ represent the number of the susceptible prey, infective prey, and predator at time $t$, respectively. $r$ denotes the increase rate, $K$ is the environmental capacity, $\beta$ denotes the disease transmission coefficient, and $d_{1}$ represents the disease related mortality rate of $I(t)$. $k$ is the conversion efficiency of predator, while $d_{2}$, $c$ represent the natural mortality rate and the attack rate on infected prey of $P(t)$, respectively. $B_{i}(t) > 0$ ($i = 1, 2, 3$) represents the standard Brownian motion with an intensity of $i$, and $\sigma_{i} > 0$ ($i = 1, 2, 3$) represents the intensities of environmental white noise. $B_{i}(t) > 0$ ($i = 1, 2, 3$) are independent of each other.

In this paper, the different timescales are applied to the model due to multiple timescales revealing that the rate of reproduction and death of prey is often much faster than that of predators in actual situations. We take into account the timescale $\varepsilon$, which is selected as the ratio of the product of the predator's predation rate and predation conversion rate to the intrinsic growth rate of prey. Rescale the variables as

and assume

still using $t$ to represent $\tilde{t}$ and $P$ to represent $\tilde{P}$.

Moreover, considering that the susceptible prey and the infected prey are often subjected to the same interference as the same biological population in the same environment, the following dimensionless system can be obtained:

where $\varepsilon$ describes the timescale separation and $0 < \varepsilon \ll 1$. $a$ represents the interspecific competition coefficient of $I(t)$. $B_{i}(t) > 0$ ($i = 1, 2$) represents the standard Brownian motion with an intensity of $i$ while $B_{1}(t)$ and $B_{2}(t)$ are independent of each other. $\sigma_{i} > 0$ ($i = 1, 2$) represents the intensities of environmental white noise. This model is applicable to infectious diseases that spread only in the prey population and have a high mortality rate while having no effect on predators, for example, mouse hepatitis virus and rabbit hemorrhagic disease virus.

3.   Equilibrium analysis of model without noise

When $\sigma_{1} = \sigma_{2} = 0$, the corresponding deterministic model of system (2.1) is

In this section, the dynamics of the deterministic situation is mainly focused.

By calculating, there are four equilibrium points for system (3.1): $E_{1}(0, 0, 0)$, $E_{2}(K, 0, 0)$, $E_{3}(\frac{K(\delta+q\alpha)}{q^{2}K+\delta}, \frac{qK-\alpha}{q^{2}K+\delta}, 0)$, and $E_{4}(S^{*}, I^{*}, P^{*})$, where $S^{*} = K(1-qh)$, $I^{*} = h$, $P^{*} = qK(1-qh)-\alpha -\delta h$.

It is easy to see that the Jacobian matrix of (3.1) is given as

For $E_{1}(0, 0, 0)$, the corresponding Jacobian matrix is

Obviously, the eigenvalues of $J(E_{1})$ are $\lambda_{1} = 1 > 0$, $\lambda_{2} = -\alpha < 0$, $\lambda_{3} = -\varepsilon h < 0$. Therefore, $E_{1}$ is a saddle and always unstable.

Theorem 3.1. For $E_{2}$,

(a) if $q < \frac{\alpha}{K}$, $E_{2}$ is a stable node.

(b) if $q > \frac{\alpha}{K}$, $E_{2}$ is a saddle.

(c) if $q = q_{TC} = \frac{\alpha}{K}$, the system (3.1) undergoes a transcritical bifurcation.

Proof. For $E_{2}\left(K, 0, 0 \right)$,

the eigenvalues of $J\left(E_{2}\right)$ are $\lambda_{1} = -1 < 0$, $\lambda_{2} = -h\varepsilon < 0$, $\lambda_{3} = qK-\alpha$. If $q < \frac{\alpha}{K}$, $\lambda_{3} < 0$, $E_{2}$ is stable. If $q > \frac{\alpha}{K}$, $\lambda_{3} > 0$, $E_{2}$ is a saddle.

If $q = q_{TC} = \frac{\alpha}{K}$, $\lambda_{3} = 0$. Let $V$ and $W$ be two eigenvectors corresponding to the eigenvalue $\lambda_{3}$ for the matrices $J\left(E_{2}\right)$ and $J\left(E_{2}\right)^{\mathrm{T}}$, respectively. After calculation, we have

Moreover,

Clearly, $V$ and $W$ satisfy

By the Sotomayor theorem ^[13], when $q = q_{TC}$, the transcritical bifurcation occurs at $E_{2}$, and the two equilibrium points $E_{2}$ and $E_{3}$ coincide. The proof of Theorem 3.1 is finished.

Theorem 3.2. For $E_{3}$,

(a) if $h > \frac{qK-\alpha}{q^{2}K+\delta}$, $E_{3}$ is a stable node.

(b) if $h < \frac{qK-\alpha}{q^{2}K+\delta}$, $E_{3}$ is a saddle.

(c) if $h = h_{TC} = \frac{qK-\alpha}{q^{2}K+\delta}$, the system (3.1) undergoes a transcritical bifurcation.

To ensure the existence of $E_{3}$, $qK > \alpha$ must be satisfied, and $h$ is chosen to be the bifurcation parameter for ease of calculation. Then, the situation is similar to $E_{2}$. By calculating, the following results are obtained:

By the Sotomayor theorem ^[13], when $h = h_{TC}$, the transcritical bifurcation occurs at $E_{3}$, and the two equilibrium points $E_{3}$ and $E_{4}$ coincide.

For $E_{4}$, its existence conditions are $q > \frac{\alpha}{K}$, $h < \frac{1}{q}$, $h < \frac{qK-\alpha}{q^{2}K+\delta}$. The corresponding Jacobian matrix is

and the characteristic equation is

where

Due to the existence conditions, it is easy to know $a_{1}, a_{2}, a_{3} > 0$ and $a_{1}a_{2}-a_{3} > 0$. According to the Routh-Hurwitz criterion, $E_{4}$ is asymptotically stable.

4.   Analysis of model with noise

4.1.   The existence and uniqueness of positive solutions

Considering the nonnegative nature of the ecological population, it is necessary to prove the existence of unique positive global solution for any positive initial condition of the system (2.1).

Theorem 4.1. For any initial value $(S(0), I(0), P(0)) \in R^{3}_{+}$, system (2.1) admits a unique global positive solution $(S(t), I(t), P(t))$ for all $t \geq 0$, and the solution remains in $R^{3}_{+}$ with probability one, namely,

Proof. Due to the coefficients of the system (2.1) satisfying the local Lipschitz condition, it always has solution $(S(t), I(t), P(t)) \in R^{3}_{+}$ at $t \in [0, \tau_{e})$ for any initial value $(S(0), I(0), P(0)) \in R^{3}_{+}$, where $\tau_{e}$ is the moment of explosion. To prove that it is almost a global solution, it's simply needed to show that $\tau_{e} = \infty$. Let $k_{0} > 0$ be sufficiently large such that each component of $(S(0), I(0), P(0))$ lies within the interval $[\frac{1}{k_{0}}, k_{0}]$. For any integer $k > k_{0}$, the stopping time is defined:

where $inf \emptyset = \infty$, and $\emptyset$ generally represents the empty set.

Obviously, $\tau_{k}$ monotonically increases as $k$ increases. Set $\tau_{\infty} = \lim\limits_{k\to \infty}\tau_{k}$, then $\tau_{\infty} \leq \tau_{e}$ is valid almost surely. It is known that if $\tau_{\infty} = \infty$ holds almost surely, then $\tau_{e} = \infty$ will hold almost surely and there is $(S(t), I(t), P(t)) \in R^{3}_{+}$ for any $t \geq 0$. The following proof will use the method of proof by contradiction. If the statement is not true, there exist a pair of constants $T > 0$ and $\varepsilon_{1} \in (0, 1)$ such that $P\{\tau_{k} \leq T\} \geq \varepsilon_{1}$ for any integer $k \geq k_{0}$.

Define a $C^{2}$ function $V: R^{3}_{+} \rightarrow R_{+}$

from inequality $u-1-lnu \geq 0$($u > 0$), it can be inferred that $V(S(t), I(t), P(t)) \geq 0$. Applying the generalized $It\hat{o}$ formula to $V$, we have

By using inequality $x-a-aln\frac{x}{a} \geq 0$ for any positive constant $a$, and combining formula (4.1) with positive constant $\xi = 2\varepsilon K+1 > 0$, we get

where $\overline{K} = \xi |ln\frac{\xi}{\varepsilon K}|+K+1+\alpha+\delta K+\varepsilon h+\sigma_{1}^{2}+\frac{\varepsilon \sigma_{2}^{2}}{2} > 0$, and further obtain

Set

and obtain that

For any $k\geq k_{0}$ and $t\in [0, T]$, integrating the two sides of the inequality (4.4) from 0 to $\tau_{k} \wedge T = min\{\tau_{k}, T\}$, we obtain

where

are three local martingales. Then, take expectation of the inequality (4.5) due to the fact that the solution of the system (2.1) is $\mathcal{F}_{t}$ adaptive, and we have

Based on the Gronwall inequality, it yields that

Therefore, for any $k \geq k_{0}$, we define $\Omega_{k} = \{\omega \in \Omega_{k} :\tau_{k} = \tau_{k}(\omega) \leq T\}$, then we have $P(\Omega_{k}) \geq \varepsilon_{1}$. Note that for each $\omega \in \Omega_{k}$, there is at least one in $S(\tau_{k}, \omega)$, $I(\tau_{k}, \omega)$, $P(\tau_{k}, \omega)$ that equals $k$ or $\frac{1}{k}$. Consequently, we have

It can be concluded from formula (4.7):

where $1_{\Omega_{k}(\omega)}$ is the indicator function of $\Omega_{k}$. It's clear to have the expression $\infty = \tilde{V}(S(0), I(0), P(0))e^{\xi T} < \infty$ as $k$ tends to infinity, which shows contradiction. Accordingly, there exists a unique positive solution for stochastic system (2.1), and the proof is complete.

4.2.   Boundedness

Theorem 4.2. For any initial value $(S(0), I(0), P(0)) \in R^{3}_{+}$, the solutions $U(t) = (S(t), I(t), P(t))$ of model (2.1) are stochastically ultimately bounded.

Proof. To prove the validity of the property, it is to prove that for any $0 < \varepsilon_{1} < 1$, there is a positive constant $\tilde{\delta} = \tilde{\delta}(\varepsilon_{1})$ such that the solution $U(t) = (S(t), I(t), P(t))$ satisfies

for any initial value $(S(0), I(0), P(0)) \in R_{+}^{3}$.

We define a function $V:$

where $(S, I, P) \in R_{+}^{3}$ and $\theta > 1$.

Applying the generalized $It\hat{o}$ formula to $V$:

let $K-h+\frac{1}{\theta} = \frac{(\theta-1)}{2}\sigma_{2}^{2}$; further, there is

where $\tilde{K}$ is a constant.

Integrating the two sides of the inequality (4.8) from 0 to $\tau_{k} \wedge t$ and taking expectation, we have

This shows that

Meanwhile, it could obtain that

Thus, we could have

which implies that

Obviously, a positive constant $\pi_{1}$ could be found such that

Applying Markov's inequality and taking $\delta = \frac{\pi_{1}}{\varepsilon_{1}}$, for any $0 < \varepsilon_{1} < 1$, we have

Hence, the following inequality holds:

and the proof is complete.

4.3.   Ergodic property of positive recurrence

When considering eco-epidemic models, when the disease will persist is always of interest. In this section, based on the theory of Has'minskki ^[14], an ergodic stationary distribution which reveals that the disease will persist is proved. Here is some theory about the stationary distribution.

Lemma 4.1. (^[14]) The Markov process $X(t)$ has a unique ergodic stationary distribution $\mu(\cdot)$ if there exists a bounded domain $D \subset E_{l}$ with regular boundary $\Gamma$ and

$A_{1}$: there is a positive number $M$ such that $\sum_{i, j = 1}^{l}a_{ij}(x)\xi_{i}\xi_{j} \geq M|\xi|^{2}$, $x \in D$, $\xi \in R^{l}$.

$A_{2}$: there exists a nonnegative $C^{2}$-function $V$ such that $LV$ is negative for any $E_{l}\backslash D$. Then,

for all $x \in E_{l}$, where $f(\cdot)$ is a function integrable with respect to the measure $\mu$.

Theorem 4.3. Assume that $R_{0}^{s} = \frac{1}{\alpha +\sigma_{1}^{2}+\varepsilon h+\frac{\varepsilon}{2}\sigma_{2}^{2}} > 1$, then system (2.1) has a unique stationary distribution $\mu(\cdot)$ and it has the ergodic property.

Proof. In view of Theorem 4.1, it has been known that for any initial value $(S(0), I(0), P(0)) \in R_{+}^{3}$, there exists a unique global solution $(S(t), I(t), P(t)) \in R_{+}^{3}$. In what follows, for the simplification, we denote $S(t), I(t)$, and $P(t)$ as $S, I, P$, respectively.

The diffusion matrix of system (2.1) is given by

Choose $M = \min \limits_{(S, I, P) \in D_{k} \subset R_{+}^{3}} \{\sigma_{1}^{2}S^{2}, \sigma_{1}^{2}I^{2}, \varepsilon \sigma_{2}^{2}P^{2}\}$, and we have

where $D_{k} = (\frac{1}{k}, k)\times(\frac{1}{k}, k)\times(\frac{1}{k}, k)$, then the condition $A_{1}$ in Lemma 4.1 holds.

Next, we focus on the condition $A_{2}$, and define

where $m$ is a sufficiently small constant satisfying $0 < m < min\{\frac{\alpha}{\sigma_{1}^{2}}-1, \frac{h}{\sigma_{2}^{2}\varepsilon}-1\}$.

By calculating, we have

where

Similarly, it can be obtained that

where

Define a Lyapunov function $\overline{V}:R_{+}^{3} \rightarrow R$ as

where $M_{0}$ is a positive constant satisfying

It is easy to check that

Furthermore, $\overline{V}(S, I, P)$ is a continuous function. Hence, $\overline{V}(S, I, P)$ has a minimum point $(S_{0}, I_{0}, P_{0})$ in the interior of $R_{+}^{3}$. Then, a nonnegative $C^{2}$-function $\hat{V}(S, I, P):R_{+}^{3} \rightarrow R_{+}$ is constructed as follows:

From (4.10) to (4.15), it can be calculated that

Now, we select a bounded closed set $D_{\varepsilon_{1}}$ as

where $\varepsilon_{1} > 0$ is a sufficiently small constant satisfying the following conditions in the set $R_{+}^{3} \backslash D_{\varepsilon_{1}}$:

where

We can get that (4.17) holds from (4.14). For the sake of convenience, $R_{+}^{3}\backslash D_{\varepsilon_{1}}$ is divided into the following six domains:

Obviously, $D_{\varepsilon_{1}}^{C} = D_{1}\cup\cdots \cup D_{6}$. In order to verify $L\hat{V}(S, I, P) \leq -1$ for any $(S, I, P)$ in $D_{\varepsilon_{1}}^{C}$, we will clarify it on the above six domains.

Case 1. For $(S, I, P) \in D_{1}$, it follows from (4.16) and (4.17) that

Case 2. For $(S, I, P) \in D_{2}$, it follows from (4.16) and (4.18) that

Notice that for any $p > 1$, $x \geq 0$, the following inequality holds ^[15]:

thus,

Case 3. For $(S, I, P) \in D_{3}$, it follows from (4.16), (4.19), and (4.23) that

Case 4. For $(S, I, P) \in D_{4}$, it follows from (4.16) and (4.20) that

Case 5. For $(S, I, P) \in D_{5}$, it follows from (4.16) and (4.21) that

Case 6. For $(S, I, P) \in D_{6}$, it follows from (4.16) and (4.22) that

Clearly, from above all it can be obtained that for a sufficiently small $\varepsilon_{1}$,

Therefore, $A_{2}$ in Lemma 4.1 is satisfied. The system (2.1) has a stable stationary distribution and the solution is ergodic by the Lemma 4.1 ^[14]. This completes the proof.

Remark 4.1. The existence of the stationary distribution describes the probability equilibrium state of variables such as biomass and species abundance in the long-term evolution of ecosystems. It provides quantitative tools for key issues such as extinction risk and coexistence probability. The ultimate boundedness characterizes the dynamic resilience of an ecosystem, which refers to its ability to resist collapse or explosive growth under stochastic disturbances (environmental noise and resource fluctuations). Both of them provide theoretical basis for the analysis of ecosystems, and ultimately boundedness often lays the foundation for the existence of stationary distributions, which is a necessary but not sufficient condition for the existence of stationary distributions.

5.   Stochastic bifurcation

In this section, first the dimensionality of deterministic model (3.1) is reduced. Then, the corresponding reduction stochastic model could be obtained by applying the lemma from ^[11]. Finally, the stochastic bifurcation is analyzed.

In terms of the slow time $t = \varepsilon \tau$, the model (3.1) becomes the following slow system:

Note that systems (3.1) and (5.1) are equivalent as long as $\varepsilon > 0$.

When $\varepsilon = 0$ in (3.1), the following equation, which is also called the fast subsystem, could be obtained:

The singular perturbation theory defines the critical manifold of (3.1) as the equilibrium points of the fast subsystem (5.2), which means

By simply calculating, we get the explicit expression of critical manifold:

For further discussion, the critical manifold is decomposed into the following branches:

Then, we have the following result:

Theorem 5.1. Consider $0 < \varepsilon \ll 1$, where the branches $M_{10}^{-}$ and $M_{20}^{-}$ are normally hyperbolic attracting, and $M_{10}^{+}$ and $M_{20}^{+}$ are normally hyperbolic saddle. $M_{10}$ and $M_{20}$ lose their normal hyperbolicity at $Q(K, 0, qK-\alpha)$ and $T(0, \frac{1}{q}, -\frac{\delta}{q}-\alpha)$.

Proof. Along the manifold $M_{10}$, the Jacobian matrix associated with (5.3) is given by

The eigenvalues are $\lambda_{11} = -1 < 0$ and $\lambda_{12} = qK-P-\alpha\ $. By simply calculating, $M_{10}$ is normally hyperbolic attracting when $P > qK-\alpha$ and normally hyperbolic saddle when $0 < P < qK-\alpha$. In addition, $Q(K, 0, qK-\alpha)$ is a turning point where $M_{10}$ loses its normal hyperbolicity (as Figure 1).

Similarly, along the manifold $M_{20}$, we have

The eigenvalues are $\lambda_{21, 22} = \frac{qI-\delta I-1 \pm \sqrt{(\delta I-qI+1)^{2} -4(qI^{2}(-q^{2}K-\delta)+I(\delta +q^{2}K))}}{2}$. By calculating, $M_{20}$ is normally hyperbolic attracting when $0 < I < \frac{1}{q}$ and normally hyperbolic saddle when $I > \frac{1}{q}$. In addition, $T(0, \frac{1}{q}, -\frac{\delta}{q}-\alpha)$ is a turning point where $M_{20}$ loses its normal hyperbolicity.

Obviously, $T$ has no biological significance because of $-\frac{\delta}{q}-\alpha < 0$. In the following, we will only discuss about point $Q$. Moreover, from definition, it is easy to know that slow flow doesn't exist at $Q$. For $0 < \varepsilon \ll 1$, Fenichel's theorem indicates that $M_{10}$ and $M_{20}$ can be perturbed to $M_{10}^{\varepsilon}$ and $M_{20}^{\varepsilon}$, which the distance to $M_{10}$ and $M_{20}$ is within $O(\varepsilon)$, respectively. As $\varepsilon$ increases, the solution trajectory fluctuates accordingly, but ultimately stabilizes at the internal equilibrium point $E_{4}$, as Figure 2.

To analyze the system (3.1), the central manifold theorem is applied to reduce the dimensionality. Using the linear transformation $w = S+qKI$ yields the following system:

Let

By transformation, the matrix $\partial_{X}f(K, 0, qK-\alpha)$ is transformed to block-diagonal with eigenvalues $-1$ and $0$. Hence, for sufficiently small $\varepsilon$ and in a neighborhood of $(0, qK-\alpha)$, the locally attracting critical manifold $\hat{M_{0}}$ and slow manifold $\hat{M_{\varepsilon}}$ are as follows:

where $h_{0}(I, P) = \frac{K(1+qI+q^{2}KI+\sqrt{(1+qI+q^{2}KI)^{2}-\frac{4}{K}(qKI+q^{3}K^{2}I^{2}+qKIP+qK\alpha I+qK\delta I+K)})}{2}$ and $h(I, P, \varepsilon)$ is a solution of the partial differential equation

and $N$ is a sufficiently small neighborhood of $(I, P)$.

Explicitly, we have

Then, we obtain the following equation:

Through using Maple to compare coefficients, the center manifold expression to second order is

where $k_{4} = -\frac{qK(qK\delta-\alpha \delta+\alpha q+\delta)}{(qK-\alpha +1)^{2}(2qK-2\alpha +1)}$.

By local attractivity, it is sufficient to obtain the reduced system of system (3.1):

The function $h(I, P, \varepsilon)$ satisfies $h(0, qK-\alpha, 0) = K$. In order to make it easier to distinguish between the stochastic system and deterministic system, we use $w_{t}$, $I_{t}$, and $P_{t}$ to represent stochastic system variables here. Meanwhile, we fix a particular solution $(I_{t}^{det}, P_{t}^{det})$ of the deterministic system (5.5).

By local attractivity, the system (2.1) could be rewritten in $(w_{t}, I_{t}, P_{t})$-coordinates as

where

Consider the deviation $\xi _{t} = w_{t}-h(I_{t}, P_{t}, \varepsilon)$ of sample paths from $\hat{M_{\varepsilon}}$. It satisfies an SDE (Stochastic Differential Equation) of the form

and using It$\hat{o}$'s formula, it shows that

where the last term is the second-order term in It$\hat{o}$'s formula. The properties of invariant manifold indicate that all terms in $\hat f^{-}(\xi_{t}, I_{t}, P_{t})$, except the last one, vanish when $\xi_{t} = 0$. Then, we could approximate the linearization of $\hat{f^{-}}$ at $\xi_{t} = 0$ by the matrix

Since $A^{-}(0, qK-\alpha, 0) = \partial_{w_{t}}f^{-}(K, 0, qK-\alpha) = A^{-} = -1$, the eigenvalues of $A^{-}(I_{t}, P_{t}, \varepsilon)$ have negative real parts, bounded away from zero when we take $N$ and $\varepsilon$ small enough. In the following, the dependence of $A^{-}$ on $\varepsilon$ will be ignored.

Now approximate the dynamics of $(\xi_{t}, I_{t}, P_{t})$ by the following system:

where $F_{0}^{-}(I, P) = F^{-}(0, I, P)|_{Q} = K$ is the value of the diffusion coefficient on the invariant manifold. Moreover, the process $\xi_{t}^{0}$ is Gaussian with zero mean and covariance matrix $X$, where $X$ obeys the equation:

through calculating, we have

As we know, all eigenvalues of $A^{-}(I_{t}, P_{t})$ have negative real parts, and there exists a locally attracting invariant manifold ${X} = {H}(I, P, \varepsilon)$ for $(I, P)\in N$ and $\varepsilon$ small enough. Based on this invariant manifold, define the domain of concentration of paths

and stopping times

where $h_{*}$ is defined in Theorem 5.3.2 ^[11].

Then, a conclusion through the Theorem 5.3.2 could be drawn that sample paths of the stochastic system (2.1) are concentrated in $B(h)$ as long as $(I_{t}, P_{t})$ remains in $N$. To put it another way, it means the sample paths under stochastic condition tend to concentrate in a neighborhood of order $\sigma_{1}$ of the invariant manifold $\hat{M_{\varepsilon}}$. Based on this result, rescale time by $\tau = \frac{1}{\varepsilon}t$, the system could be approximated by its projection

which is also called the reduced stochastic system.

Remark 5.1. This method of reducing the dimensionality of the stochastic system by mapping it to the slow manifold makes the dynamical behavior of multidimensional stochastic slow-fast systems easier to analyze, while also ensuring the accuracy of the analysis. For ecosystems affected by natural stochastic noise and characterized by multiple timescales, the discussion process will not be complicated and the conclusions are more in line with reality.

In the following, the bifurcation of the reduced system (5.6) will be discussed. First of all, there are three equilibrium points containing two boundary equilibrium points $E_{11}(0, 0)$, $E_{12}(I_{12}, 0)$ where $I_{12} = \frac{q^{2}K+\delta -qk_{1}+\sqrt{(qk_{1}-q^{2}K-\delta)^{2}-4qk_{4}(qK-\alpha)}}{2qk_{4}}$, and an internal equilibrium point $E_{13}(I_{r}, P_{r})$ where $I_{r} = h, P_{r} = \frac{q(k_{0}+k_{1}h+k_{4}h^{2})-q^{2}Kh-\alpha -\delta h}{1-k_{5}qh}$. Note that $0 < h < \frac{q^{2}K+\delta -qk_{1}+\sqrt{(qk_{1}-q^{2}K-\delta)^{2}-4qk_{4}(qK-\alpha)}}{2qk_{4}}$ is an existence condition for $E_{13}$.

This subsection mainly focuses on the analysis of the internal equilibrium point $E_{13}$, as it has more biological significance. Let $u = I_{t}-I_{r}$, $v = P_{t}-P_{r}$, and we have

where $\xi_{i}(\tau)d\tau = dB_{i}(\tau)$,

Performing coordinate transformation $u = rcos\theta$ and $v = rsin\theta$ on system(5.7), we have

Based on the Khasminskii limit theorem ^[14], if the noise intensity of white noise $(\sigma_{1}, \sigma_{2})$ is small enough, then the response process ${r(\tau), \theta(\tau)}$ weakly converges with a two-dimensional Markov diffusion process. By using the random averaging method, the following It$\hat{o}$ stochastic differential equation is obtained:

where $B_{r}(\tau)$ and $B_{\theta}(\tau)$ are independent standard Wiener processes, $\left(\begin{matrix} m_{r}\\ m_{\theta}\\ \end{matrix} \right)$ is the drift vector, and $\left(\begin{matrix} \varepsilon_{11} & \varepsilon_{12}\\ \varepsilon_{21} & \varepsilon_{22}\\ \end{matrix} \right)$ is the diffusion coefficient matrix.

By calculating various coefficients, the parameters in (5.8) satisfy:

where the average amplitude $r(\tau)$ is a one-dimensional Markov diffusion process when $\varepsilon_{12}^{2} = \varepsilon_{21}^{2} = 0$. It means we could obtain the following equation:

where

From $P_{r} > 0$, it can be seen that $\mu_{3} \neq 0$. According to the It$\hat{o}$'s differential equation, the Fokker-Planck equation of system (5.9) is:

The initial condition is $P(r, \tau | r_{0}, \tau_{0}) \rightarrow \delta(r-r_{0})$, $\tau \rightarrow \tau_{0}$, where $P(r, \tau | r_{0}, \tau_{0})$ is the transition probability density of the diffusion process $r(\tau)$. The steady-state density $P_{ST}(r)$ is the invariant measure of $r(\tau)$, and it is the solution of the following degenerate system:

By calculating, we could obtain

where $\Delta = (8\mu_{1}+\mu_{2})\mu_{4}^{-1}$, $\Gamma(x) = \int_{0}^{\infty} \tau^{x-1} e^{-\tau} d\tau$.

Based on Namachivaya's theory ^[16], when the noise intensity approaches zero, the extreme value of $P_{st}(r)$ approaches the behavior of the deterministic system. We calculate the most likely amplitude $r^{*}$ of (5.9) so that $P_{st}(r)$ has its maximum value at $r^{*}$, i.e.,

where $r^{*} = \sqrt{\frac{-8\mu_{3}}{8\mu_{1}+\mu_{2}-\mu_{4}}}$, ($\frac{-8\mu_{3}}{8\mu_{1}+\mu_{2}-\mu_{4}} < \frac{1}{2}$).

In addition, $P_{st}(r)$ achieves the minimum value at $r = 0$, which means the system (2.1) is almost unstable at the equilibrium point $E_{13}$ when subjected to random excitation. According to the singular boundary theory, the system (2.1) undergoes P-bifurcation at $r = r^{*}$.

Remark 5.2. Through numerical simulation, as shown in Figure 3, the position of P-bifurcation increases as $\mu_{3}$ increases. Figure 4(a)–(d) respectively corresponds to different values of parameter $\mu_{3}$ in Figure 3(a). The influence of $\varepsilon$ is also shown in Figure 3. With the increase of $\varepsilon$, the oscillation amplitude of $x(\tau)$ and $y(\tau)$ decreases, where the impact on $x(\tau)$ is relatively greater.

Remark 5.3. The existence of P-bifurcation indicates that the solution crossing the singularity $Q$ will follow $M_{10}^{\varepsilon +}$ for a distance before being repelled, then the solution will be attracted by $M_{20}^{\varepsilon -}$. By calculating, the equilibrium point $E_{4}$ on $M_{20}^{\varepsilon -}$ is locally asymptotically stable. Hence, the solution converges to the stable equilibrium point $E_{4}$.

6.   Conclusions

In this paper, an stochastic eco-epidemiological model with multi-time scale is mainly studied through the theory analysis and numerical simulation. First, we take a dimensionless transformation to system for simplicity. Then based on the stochastic noise, the stochastic items are considered in order to be more in line with reality. Meanwhile, considering the difference in iteration speed between prey and predator, paramater $\varepsilon$ is added to the model. The stability and bifurcation of equilibrium points on the deterministic system are discussed.

For the system with noise, the existence and uniqueness of solutions is proven by constructing a function $V$, and boundedness is also obtained. Then, the $C^{2}$-function $\overline{V}$ and $\hat{V}$ are constructed to prove the ergodic property of the solutions by the lemma from ^[14]. As for the slow-fast dynamics, the deterministic condition is analyzed. It is obtained that the solution crossing the singularity converges to the stable equilibrium point $E_{4}$ by geometric singular perturbation theory and center-manifold reduction considering $0 < \varepsilon \leq 1$. Applying the theorem from ^[11], it could be known that sample paths of the stochastic system are concentrated in the neighborhood of the deterministic one. Finally, through calculating, P-bifurcation occurs at the singularity, and that's why there is a delay in the trajectory at the singularity of the solution. Numerical simulation verified our theoretical results.

Compared with other eco-epidemic systems, it should be noted that the models proposed in this paper investigate dynamical behavior with multi-time scale and stochastic noises, which makes the work studied in this paper have some new and positive features. It is also interesting to consider the effects of the other types of functional response functions, which will be the subject of our further research. In practical application, the control of disease is also very significant and important. The discussion of biological control is very extensive. One direction that can be explored in the future is model predictive control. It can better tackle the scenario involving the uncertainties or disturbances, and it is worth exploring whether the model predictive control can be used to the proposed model, such as ^[17].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This study was funded by National Natural Science Foundation of China (61703083) and Regional Joint Key Project of National Natural Science Foundation (U23A20324).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.