Dynamic analysis of stochastic delay mutualistic system of leaf-cutter ants with stage structure and their fungus garden

Abstract

In this paper, we propose a stochastic delay mutualistic model of leaf-cutter ants with stage structure and their fungus garden, in which we explore how the discrete delay and white noise affect the dynamic of the population system. The existence and uniqueness of global positive solution are proved, and the asymptotic behaviours of the stochastic model around the positive equilibrium point of the deterministic model are also investigated. Furthermore, the sufficient conditions for the persistence of the population are established. Finally, some numerical simulations are performed to show the effect of random environmental fluctuation on the model.

Keywords:

• Stochastic mutualistic system
• discrete delay
• Lyapunov function
• stability

1. Introduction

Mutualism is a universal phenomenon in the ecosystem, which refers to the cooperative survival relationship among species. For example, euryhaline tilapia co-exist and evolve with chlorella, the similiar relationship also happens to pollinating fine moths and Euphorbiacecae plants, legumes and rhizobium, and so on. Sometimes, predation behaviours also lead to reciprocal effect. Leaf mites can increase cucumber yield by consuming some cucumber seedlings. These interesting phenomena have attracted the attention of scholars. In 1925, Lotka and Volterra [16, 20] have successively established a mathematical model of the interaction between biological populations, which was known as Lotka–Volterra model. After that, many scholars have studied the dynamic behaviour of Lotka–Voterra cooperative model. For example, Yao et al. [24] added impulse and harvest term on the basis of previous studies, and generalized the research results of Lotka–Volterra cooperation model by coincidence degree theory and extension theorem of inequality analysis in 2006. Then in 2017, Feng et al. [8] supplemented delay and impulse to Lotka–Volterra cooperation model, by constructing appropriate Lyapunov functions, sufficient conditions for the existence and uniform asymptotic stability of periodic solutions of the system are obtained. It is worth mentioning that there are also plenty of models based on reality that are different from Lotka–Volterra model, for instance, Kang et al. [12] established a model of leaf-cutter ants and their fungus garden, which introduces in detail about how they cooperate with each other. The model is as follows: (1) $\{\begin{array}{rl}\frac{\mathrm{d}A}{\mathrm{d}t}& =(r_{a}F-d_{a}A)A,\\ \frac{\mathrm{d}F}{\mathrm{d}t}& =(\frac{r_{f}a{A}^2}{b+a{A}^2}-d_{f}F-r_{c}A)F,\end{array}$(1) where $A(t)$ is the total biomass of ants including workers, larvae, pupae and eggs at time t, $F(t)$ is the total biomass of the fungus at time t, $r_c=c{r}_a$ and $a=p^{2}q(1-p)$ can be considered as a parameter measuring the division of labour in the colony of ants. The meaning of other parameters is shown in Table 1, all parameters are positive constants.

Table 1. The meaning of some parameters.

Parameters	Biological meaning	Parameters	Biological meaning
$r_a$	Maximum growth rate of ants	$d_f$	Death rate of fungus
$r_f$	Maximum growth rate of fungus	p	Proportion of ants that are workers
c	Conversion rate between fungus and ants	q	Proportion of workers that take care of fungus
$d_a$	Death rate of ants	b	Half-saturation constant

The mutualism relationship between leaf-cutter ants and their fungus garden that we intend to study in this paper is very fascinating. The ant is very good at cutting leaves with its tail, hence it was named the ‘leaf-cutter ant’. In the Amazon rainforest, they carry leaves back to their nest not for food, but to plant fungus. To supply a place with good environment for fungus' growth, the leaves will be cut into smaller pieces and stack together after they are moved into the nest, then, under the meticulous care of the leaf-cutter ants, the mycelia quickly grows into a fungus garden. In fact, the relationship between the leaf-cutter ants and their fungus garden just like human plant mushroom and eat them, except that they don't sell the fungus to other leaf-cutter ant tribes. After the leaf-cutter ants have scattered their mycelia, if they do not look after the mycelia carefully, such as weeding, fertilizing and spraying, they will not produce a good harvest, exactly as humans grow crops.

The environment inside the nest is very favourable for the growth of the fungus, which also creates hotbeds for the growth of miscellaneous bacteria. Obviously, leaf-cutter ants face two difficulties during the cultivation process. One is how to ensure that the fungus is not contaminated by other bacteria. The other is how do they protect themselves from microbes in an environment surrounded by so many bacteria. In fact, leaf-cutter ant's mouth and forelimb are infested with parasites that can produce streptomycin, they can adjust the active ingredient of streptomycin automatically according to the species of harmful bacteria that have been present. The infected fungus will be picked out and threw away by the ants, if streptomycin does not kill the bacteria. The ants are surrounded by a wide range of pathogenic bacteria all the time, why the ants are not infected by microbes is that they have chosen to live symbiosis with a bacterium of the genus Pseudonocardia. The bacteria attach to the ant's body like armour, they secrete antibiotics once the ant is invaded by the pathogen, which help the ant fight off the intruding pathogen and mould. So far, scientists have isolated a variety of strains from the leaf-cutter ant and extracted the corresponding gene sequence to test whether the corresponding antibodies they produced are active in humans. Scientists hope to use these special bacteria to produce new antibiotics to help humans treat diseases more effectively.

In this paper, the functional response function is defined as the culture rate of fungus by leaf-cutter ants. In 1975, J. R. Beddington and D. L. DeAngelis proposed Beddington–DeAngelis type functional response $(\phi (x,y)=\frac{x}{a+bx+cy})$ [3, 6], referred to as B–D functional response, which is the combination of Holling II functional response and ratio-dependent functional response, has the advantages of prey-dependent functional response and ratio-dependent functional response. Obviously, both of the biomass of leaf-cutter ants and fungus can affect the fungus culture rate. Leaf-cutter ants would reduce the culture rate when there are enough fungus, and on the contrary, they would increase the culture rate when there are less fungus, however, the upper limit of culture rate is related to the number of worker ants among leaf-cutter ants. To have a better explanation of the dynamic behaviour between leaf-cutter ants and their fungus garden, B-D functional response is considered in this paper.

For the sake of making the model more realistic, many scholars introduce stage structure into population model on account of the different abilities and behaviours of species at different ages [1, 13, 15, 18, 21, 23]. Leaf-cutter ants are classified as the maturity and the immaturity, the adult leaf-cutter ants will participate in the cultivation of fungus, while the immaturity, including juvenile ants, eggs and pupae, need to be fed by adult leaf-cutter ants. Then, we get the following model with stage structure and B-D function response (2) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t)y(t)}{m+m_{1}y(t)+m_2{x}_2(t)}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(2) where $x_1(t),x_2(t)$ represent the biomass of immature and mature ants, respectively, $y(t)$ represents the biomass of fungus, $d_i(i=1,2,3)$ is the rate of death, $r_1$, $r_2$ are parameters that measure the maximum growth rate of ants and fungus, respectively, $r_3$ indicates the collection rate of fungus collected by adult ants, a can be considered as a parameter measuring the division of labour in the colony of ants, all parameters are strictly positive.

Time delay is an important factor in the biological mathematical model [2, 11, 19, 22]. In this paper, we introduce gestation time into the model, it is a period that the fungus begins to grow from spores. To cultivate fungus, leaf-cutter ants need to cut leaves first, after the leaves are transported back to the nest, the gardener ants will be responsible for the subsequent cultivation of the fungus garden, the leaves are cleaned and arranged in tight rows and in suitable positions by the gardener ants, and then they are decomposed by the secretions of the gardener, which are rich in enzymes and can form a mucus layer on the leaves' surface. With abundant nutrients, the mucus layer is a natural culture medium in which the gardener ants ‘seed’ a small amount of their own fungus, and then the careful care of ants makes the fungi grow. The period of the fungus from spores to starting growing is a relatively complex gestation time, which can be called ‘gestation time’, then the model (2(2) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t)y(t)}{m+m_{1}y(t)+m_2{x}_2(t)}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(2) ) with time delay can be described as (3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) where τ represents the gestation delay.

An ecosystem is inevitably disturbed by random fluctuations such as temperature, humidity, rainfall, and sometimes even artificial factors in the wild [5, 7, 10, 11, 14, 17, 25]. The sudden torrential rain will bring devastating disaster to the leaf-cutter ants, and the drought will increase the death rate of the fungus garden, hence, it is reasonable to consider the random fluctuations in the mutualistic system and it will enrich the dynamic behaviour of the population. Therefore, we consider that the mortality of leaf-cutter ants and fungus are affected by the natural environment in our model.

Let $-d_i\to -d_i+{\sigma }_i{\dot{B}}_i(t),i=1,2,3,$where $B_i(t)(i=1,2,3)$ denote the independent standard Brownian motion, and ${\sigma }_i^2(i=1,2,3)$ denote the intensity of white noise [17, 25]. Then, the stochastic system takes the form (4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) and (5) $x_i(\theta )={\varphi }_i(\theta ),y(\theta )=\psi (\theta ),\theta \in [-\tau ,0],i=1,2,$(5) where ${\varphi }_i(\theta ),(i=1,2)$ and $\psi (\theta )$ are nonnegative continuous function on $[-\tau ,0]$, ${\varphi }_i(\theta )>0,\psi (\theta )>0,\theta \in [-\tau ,0],(i=1,2)$.

We have the following problems to be solved about this system: first, investigate whether the model has a globally unique positive solution, and how about the asymptotic behaviour of the model, then we want to know whether the sufficient conditions for the persistence of the model can be established. Naturally, we hope to analyse the influence of the environment disturbance on the model.

2. Existence and uniqueness of global positive solutions

In this section, we study the existence and uniqueness of positive solution of system, which is the basis of studying the long-time behaviour of the model. First of all, we present some known results and lemma which will be used later.

In general, consider the d-dimensional stochastic differential equation (6) $\text{d}X(t)=f(X(t),t)\text{d}t+g(X(t),t)\text{d}B(t),\mathrm{for}t\ge t_0,$(6) with the initial value $X(0)=x_0\in {\mathbb{R}}^d$. $B(t)$ denotes a d-dimensional standard Brownian motion defined on the complete probability space $(\mathrm{\Omega },\mathcal{F},\{{\mathcal{F}}_t{\}}_{t\ge 0},\mathbb{P})$. The family of all nonnegative functions $V(X,t)$ defined on ${\mathbb{R}}^d\times [t_0,\mathrm{\infty })$, such that they are continuously twice differentiable in X and once in t. The differential operator L of Equation (6(6) $\text{d}X(t)=f(X(t),t)\text{d}t+g(X(t),t)\text{d}B(t),\mathrm{for}t\ge t_0,$(6) ) is defined by Mao [17]. $L=\frac{\mathrm{\partial }}{\mathrm{\partial }t}+\sum _{i=1}^d{f}_i(X,t)\frac{\mathrm{\partial }}{\mathrm{\partial }{X}_i}+\frac{1}{2}\sum _{i,j=1}^d[g^T(X,t)g(X,t){]}_{ij}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{X}_i\mathrm{\partial }{X}_j}.$If L acts on a function $V\in {\mathbb{C}}^{2,1}({\mathbb{R}}^d\times [t_0,\mathrm{\infty });{\mathbb{R}}_{+})$, then $LV(X,t)=V_t(X,t)+V_X(X,t)f(X,t)+\frac{1}{2}\mathrm{trace}[g^T(X,t)V_{XX}(X,t)g(X,t)],$where $V_t=\frac{\mathrm{\partial }V}{\mathrm{\partial }t}$, $V_X=(\frac{\mathrm{\partial }V}{\mathrm{\partial }{X}_1},\dots ,\frac{\mathrm{\partial }V}{\mathrm{\partial }{X}_d})$, $V_{XX}=\frac{{\mathrm{\partial }}^{2}V}{\mathrm{\partial }{X}_i\mathrm{\partial }{X}_j}(i,j=1,2,\dots ,d)$. In view of Itô formula, if $X(t)\in {\mathbb{R}}^d$, then $\text{d}V(X,t)=LV(X,t)\mathrm{d}t+V_X(X,t)g(X,t)\text{d}B(t).$

Theorem 2.1

For any given initial condition (5(5) $x_i(\theta )={\varphi }_i(\theta ),y(\theta )=\psi (\theta ),\theta \in [-\tau ,0],i=1,2,$(5) ), model (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) has a unique global positive solution $(x_1(t),x_2(t),y(t))$ on $t\ge -\tau$, and the solution will remain in ${\mathbb{R}}_{+}^3$ with probability 1.

Proof.

Since the coefficients of system satisfy the local Lipschitz condition, then for any initial value (5(5) $x_i(\theta )={\varphi }_i(\theta ),y(\theta )=\psi (\theta ),\theta \in [-\tau ,0],i=1,2,$(5) ), there exists a unique local solution $(x_1(t),x_2(t))$, $y(t))\in {\mathbb{R}}_{+}^3$ on $[-\tau ,{\tau }_e)$, where ${\tau }_e$ is the explosion time. To verify that this solution is global, we only need to prove that ${\tau }_e=\mathrm{\infty }$ a.s. For this purpose, let $n_0>0$ be large enough such that every component of $({\varphi }_1(\theta ),{\varphi }_2(\theta ),\phi (\theta )),\theta \in [-\tau ,0]$ lay within the interval $[\frac{1}{n_0},n_0]$. For each integer $n>n_0$, we define stopping time: $\begin{array}{rl}{\tau }_n& =inf\{t\in [-\tau ,{\tau }_e):min\{x_1(t),x_2(t),y(t)\}\le \frac{1}{n}\mathrm{or}\\ & max\{x_1(t),x_2(t),y(t)\}\ge n\}.\end{array}$We can see that, ${\tau }_n$ is increasing as $n\to \mathrm{\infty }$, define $\mathrm{\varnothing }$ ( denotes the empty set). Let ${\tau }_{\mathrm{\infty }}=\underset{n\to \mathrm{\infty }}{lim}{\tau }_n$ and ${\tau }_{\mathrm{\infty }}\le {\tau }_e$ a.s. Hence to complete the proof, we need to show that ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$.

Define a $C^2$-function $V(x_1,x_2,y)=x_1(t)+1-\ln x_1(t)+x_2(t)+1-\ln x_2(t)+y(t)+1-\ln y(t)$, applying the $\mathrm{It}\hat{\mathrm{o}}$ formula [17], one can compute that $\begin{array}{rl}\mathrm{L}V(x_1,x_2,y)& =(1-\frac{1}{x_1(t)})(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\\ & +(1-\frac{1}{y}(t))(\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t))\\ & +(1-\frac{1}{x_2(t)})(c{x}_1(t)-d_2{x}_2(t))+\frac{{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2}\\ & =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t)-\frac{r_1{x}_2(t)y(t)}{x_1(t)}+d_1+c+c{x}_1(t)-d_2{x}_2(t)\\ & -\frac{c{x}_1(t)}{x_2(t)}+d_2+[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)\\ & -\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{y(t)[m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )]}+d_3+r_3{x}_2(t)]\\ & +\frac{{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2}\\ & \le d_1+c+d_2+d_3+(r_1-r_3)x_2(t)y(t)+(r_3-d_2)x_2(t)\\ & +\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}+\frac{{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2}.\end{array}$Then we have $K_1,K_2$ and $K_3$, where $K_1=max\{d_1+c+d_2+d_3+\frac{{\sigma }_1^2}{2}+\frac{{\sigma }_2^2}{2}+\frac{{\sigma }_3^2}{2},2\mid r_3-d_2\mid \}$, $K_2=\mid r_1-r_3\mid$, $K_3=\frac{r_{2}a}{m}$. (7) $\mathrm{L}V\le K_1(1+\frac{x_2(t)}{2})+K_2{x}_2(t)y(t)+K_3{x}_2(t-\tau )y(t-\tau ).$(7) Note that $x_2(t)\le 2V(x_1(t),x_2(t),y(t))$. We can get $\begin{array}{rl}\mathrm{d}V& \le [K_1(1+V(x_1(t),x_2(t),y(t)))+K_2{x}_2(t)y(t)+K_3{x}_2(t-\tau )y(t-\tau )]\mathrm{d}t\\ & +{\sigma }_1(x_1(t)-1)\mathrm{d}{B}_1(t)+{\sigma }_2(x_2(t)-1)\mathrm{d}{B}_2(t)+{\sigma }_3(x_3(t)-1)\mathrm{d}{B}_3(t).\end{array}$For any $n\ge n_0$ and $t_1\in [0,\tau ]$, we obtain (8) $\mathbb{E}V(x_1(t_1\wedge {\tau }_n),x_2(t_1\wedge {\tau }_n),y(t_1\wedge {\tau }_n))\le K_4+K_1\mathbb{E}{\int }_0^{t_1\wedge {\tau }_n}V(x_1(t),x_2(t),y(t))\mathrm{d}t,$(8) where $\begin{array}{rl}{K}_4& =V({\varphi }_1(\theta ),{\varphi }_2(\theta ),\psi (\theta ))+K_1\tau +K_2{\int }_0^{\tau }{x}_2(t)y(t)\mathrm{d}t\\ & +K_3{\int }_0^{\tau }{x}_2(t-\tau )y(t-\tau )\mathrm{d}t<\mathrm{\infty }.\end{array}$Equation (8(8) $\mathbb{E}V(x_1(t_1\wedge {\tau }_n),x_2(t_1\wedge {\tau }_n),y(t_1\wedge {\tau }_n))\le K_4+K_1\mathbb{E}{\int }_0^{t_1\wedge {\tau }_n}V(x_1(t),x_2(t),y(t))\mathrm{d}t,$(8) ) and the Gronwall inequality [4, 17] imply that $\begin{array}{rl}& \mathbb{E}V(x_1(t_1\wedge {\tau }_n),x_2(t_1\wedge {\tau }_n),y(t_1\wedge {\tau }_n))\le K_4{e}^{\tau K_1},\\ & \mathrm{for}0\le t_1\le \tau ,n\ge n_0.\end{array}$It then follows that (9) $\mathbb{E}V(x_1({\tau }_n\wedge \tau ),x_2({\tau }_n\wedge \tau ),y({\tau }_n\wedge \tau ))\le K_4{e}^{\tau K_1}\mathrm{for}n\ge n_0.$(9) We can hence show that ${\tau }_{\mathrm{\infty }}\ge \tau a.s.$ Additionally, letting $n\to \mathrm{\infty }$ in Equation (9(9) $\mathbb{E}V(x_1({\tau }_n\wedge \tau ),x_2({\tau }_n\wedge \tau ),y({\tau }_n\wedge \tau ))\le K_4{e}^{\tau K_1}\mathrm{for}n\ge n_0.$(9) ) gives (10) $\mathbb{E}V(x_1(t_1),x_2(t_1),y(t_1))\le K_4{e}^{\tau K_1}\mathrm{for}0\le t_1\le \tau .$(10) For $t_2\in (\tau ,2\tau ]$, (11) $\mathbb{E}V(x_1(t_2\wedge {\tau }_n),x_2(t_2\wedge {\tau }_n),y(t_2\wedge {\tau }_n))\le K_5+K_1\mathbb{E}{\int }_0^{t_2\wedge {\tau }_n}V(x_1(t),x_2(t),y(t))\mathrm{d}t,$(11) where $K_5=V({\varphi }_1(\theta ),{\varphi }_2(\theta ),\psi (\theta ))+2K_1\tau +K_2{\int }_{\tau }^{2\tau }{x}_2(t)y(t)\mathrm{d}t+K_3{\int }_{\tau }^{2\tau }{x}_2(t-\tau )y(t-\tau )\mathrm{d}t<\mathrm{\infty }$. Similarly, we obtain that ${\tau }_{\mathrm{\infty }}\ge 2\tau a.s.$ and $\mathbb{E}V(x_1(t_2\wedge {\tau }_n),x_2(t_2\wedge {\tau }_n),y(t_2\wedge {\tau }_n))\le K_5{e}^{2\tau K_1}.$Repeating this procedure, one can show ${\tau }_{\mathrm{\infty }}\ge m\tau$ with probability one for any inter $m\ge 1$. Therefore, ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ a.s.

3. Asymptotic property

After introducing white noise into the deterministic model, the solution of the system will be disturbed under the influence of random fluctuations. Therefore, we will discuss the asymptotic property of stochastic system (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) around positive equilibrium point of the corresponding deterministic model (3(3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) ).

System (3(3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) ) has the same equilibrium point as system (2(2) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t)y(t)}{m+m_{1}y(t)+m_2{x}_2(t)}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(2) ). Hence, we just need to consider the following equations: (12) $\{\begin{array}{rl}& r_1{x}_2^{\ast }{y}^{\ast }-d_1{x}_1^{\ast }-c{x}_1^{\ast }=0,\\ & c{x}_1^{\ast }-d_2{x}_2^{\ast }=0,\\ & \frac{r_{2}a{x}_2^{\ast }{y}^{\ast }}{m+m_1{y}^{\ast }+m_2{x}_2^{\ast }}-d_3{y}^{\ast }-r_3{x}_2^{\ast }{y}^{\ast }=0.\end{array}$(12) We can deduce $x_1^{\ast }=\frac{d_2{x}_2^{\ast }}{c},$ from the second equation of (12(12) $\{\begin{array}{rl}& r_1{x}_2^{\ast }{y}^{\ast }-d_1{x}_1^{\ast }-c{x}_1^{\ast }=0,\\ & c{x}_1^{\ast }-d_2{x}_2^{\ast }=0,\\ & \frac{r_{2}a{x}_2^{\ast }{y}^{\ast }}{m+m_1{y}^{\ast }+m_2{x}_2^{\ast }}-d_3{y}^{\ast }-r_3{x}_2^{\ast }{y}^{\ast }=0.\end{array}$(12) ), we get $y^{\ast }=\frac{d_2(d_1+c)}{r_{1}c}$.

From the third equation of (12(12) $\{\begin{array}{rl}& r_1{x}_2^{\ast }{y}^{\ast }-d_1{x}_1^{\ast }-c{x}_1^{\ast }=0,\\ & c{x}_1^{\ast }-d_2{x}_2^{\ast }=0,\\ & \frac{r_{2}a{x}_2^{\ast }{y}^{\ast }}{m+m_1{y}^{\ast }+m_2{x}_2^{\ast }}-d_3{y}^{\ast }-r_3{x}_2^{\ast }{y}^{\ast }=0.\end{array}$(12) ), we have $\begin{array}{rl}{r}_{2}a{y}^{\ast }{x}_2^{\ast }& =(d_3{y}^{\ast }+r_3{x}_2^{\ast }{y}^{\ast })(m+m_1{y}^{\ast }+m_2{x}_2^{\ast })\\ & =r_3{m}_2{y}^{\ast }{x}_2^{\ast 2}+r_3{y}^{\ast }(m+m_1{y}^{\ast })x_2^{\ast }+d_3{m}_2{y}^{\ast }{x}_2^{\ast }+d_3{y}^{\ast }(m+m_1{y}^{\ast }),\end{array}$that is $r_3{m}_2{x}_2^{\ast 2}+[r_3(m+m_1{y}^{\ast })+d_3{m}_2-r_{2}a]x_2^{\ast }+d_3(m+m_1{y}^{\ast })=0,$where $y^{\ast }=\frac{d_2(d_1+c)}{r_{1}c}$, let $b_1=r_3{m}_2$, $b_2=d_3{m}_2+r_{3}m+\frac{r_3{m}_1{d}_2(d_1+c)}{r_{1}c}-r_{2}a$, $b_3=d_3(m+\frac{m_1{d}_2(d_1+c)}{r_{1}c})$.

Then (13) $b_1{x}_2^{\ast 2}+b_2{x}_2^{\ast }+b_3=0.$(13) Then, if the condition $H_0:\mathrm{△}=b_2^2-4b_1{b}_3\ge 0,b_2<0$ is satisfied, we can get one or two positive solutions of (13(13) $b_1{x}_2^{\ast 2}+b_2{x}_2^{\ast }+b_3=0.$(13) ), simplicity, let's write them all as $x_2^{\ast }$. We denote all the positive equilibria as $E^{\ast }(x_1^{\ast },x_2^{\ast },y^{\ast })$.

To facilitate the study, we define the following assumptions: (14) $H_1:\{\begin{array}{rl}& d_1+\frac{c}{2}-{\sigma }_1^2>0,\\ & d_2-\frac{c}{2}-{\sigma }_2^2>0,\\ & d_3-(r_{2}a{)}^2-{\sigma }_3^2>0,\\ & r_1{r}_{2}ac-r_3{d}_2{m}_1(d_1+c)\le 0.\end{array}$(14)

Theorem 3.1

For any given initial condition (5(5) $x_i(\theta )={\varphi }_i(\theta ),y(\theta )=\psi (\theta ),\theta \in [-\tau ,0],i=1,2,$(5) ), if hypothesis $(H_0)$ and $(H_1)$ are established, the solution $(x_1(t),x_2(t),y(t))$ of system (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) has the property (15) $\underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t[(x_1(\theta )-x_1^{\ast }{)}^2+(x_2(\theta )-x_2^{\ast }{)}^2+(y(\theta )-y^{\ast }{)}^2]\mathrm{d}\theta \le \frac{q}{p},$(15) where $p=min\{a_1(d_1+\frac{c}{2}-{\sigma }_1^2),(d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2),a_3(d_3-(r_{2}a{)}^2-{\sigma }_3^2)\},$and $q=a_1{\sigma }_1^2{x}_1^{\ast }+{\sigma }_2^2{x}_2^{\ast }+a_3{y}^{\ast }{\sigma }_3^2+\frac{a_3{x}_2^{{}^{\ast }2}}{m^2}+\frac{a_3(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}$, $(x_1^{\ast },x_2^{\ast },y^{\ast })$ is positive equilibrium point of model (3(3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) ), $a_1,a_3$ are positive constants to be determined later.

Proof.

Define the function (16) $V=a_1{V}_1+V_2+a_3{V}_3+a_3{V}_4+a_5{V}_5+a_5{V}_6+V_7.$(16) where $\begin{array}{rl}& V_1=\frac{(x_1(t)-x_1^{\ast }{)}^2}{2},V_4={\int }_{t-\tau }^t\frac{(x_2(r)-x_2^{\ast }{)}^2}{2m_1^2}\mathrm{d}r,\\ & V_2=\frac{(x_2(t)-x_2^{\ast }{)}^2}{2},V_5=y,\\ & V_3=\frac{(y(t)-y^{\ast }{)}^2}{2},V_6=r_{2}a{\int }_{t-\tau }^t\frac{x_2(t)y(t)}{m+m_{1}y(t)+m_2{x}_2(t)}\mathrm{d}r,\\ & V_7=a_2{x}_1(t)+a_4{x}_2(t).\end{array}$And $a_i(i=1,2,3,4,5)$ are positive constants to be determined later.

The positive equilibrium satisfies the following properties: (17) $\{\begin{array}{rl}& r_1{x}_2^{\ast }{y}^{\ast }-d_1{x}_1^{\ast }-c{x}_1^{\ast }=0,\\ & c{x}_1^{\ast }-d_2{x}_2^{\ast }=0,\\ & \frac{r_{2}a{x}_2^{\ast }{y}^{\ast }}{m+m_1{y}^{\ast }+m_2{x}_2^{\ast }}-d_3{y}^{\ast }-r_2{x}_2^{\ast }{y}^{\ast }=0.\end{array}$(17) By $\mathrm{It}\hat{\mathrm{o}}$ formula, we obtain $\begin{array}{rl}\mathrm{L}{V}_1& =(x_1(t)-x_1^{\ast })(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))+\frac{1}{2}{x}_1^2(t){\sigma }_1^2\\ & \le -(d_1+c-{\sigma }_1^2)(x_1(t)-x_1^{\ast }{)}^2+r_1{x}_2^{\ast }{y}^{\ast }y(t)+r_1{y}^{\ast }{x}_2(t)y(t)+{\sigma }_1^2{x}_1^{\ast },\\ \mathrm{L}{V}_2& =(x_2(t)-x_2^{\ast })(c{x}_1(t)-d_2{x}_2(t))+\frac{1}{2}{x}_2^2(t){\sigma }_2^2\\ & \le -(d_2-\frac{c}{2}-{\sigma }_2^2)(x_2(t)-x_2^{\ast }{)}^2+\frac{c}{2}(x_1(t)-x_1^{\ast }{)}^2+{\sigma }_2^2{x}_2^{\ast },\\ \mathrm{L}{V}_3& =(y(t)-y^{\ast })(\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t))\\ & +\frac{1}{2}{y}^2(t){\sigma }_3^2\\ & =\left[\frac{\begin{array}{c}(m+m_1{y}^{\ast }+m_2{x}_2^{\ast })x_2(t-\tau )y(t-\tau )\\ -(m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau ))x_2^{\ast }{y}^{\ast }\end{array}}{(m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau ))(m+m_1{y}^{\ast }+m_2{x}_2^{\ast })}\right]\\ & \cdot r_{2}a(y(t)-y^{\ast })-d_3(y(t)-y^{\ast }{)}^2-r_3(y(t)-y^{\ast })(x_2(t)y(t)-x_2^{\ast }{y}^{\ast })+\frac{1}{2}{y}^2(t){\sigma }_3^2\\ & \le -(d_3-(r_{2}a{)}^2-{\sigma }_3^2)(y(t)-y^{\ast }{)}^2+\frac{y^2(t-\tau )(x_2(t-\tau )-x_2^{\ast }{)}^2}{2(m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau ){)}^2}\\ & +\frac{[(m+m_1{y}^{\ast }+m_2{x}_2^{\ast })x_2^{\ast }y(t-\tau )-(m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau ))x_2^{\ast }{y}^{\ast }{]}^2}{2(m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau ){)}^2(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}\\ & +r_3{x}_2^{\ast }{y}^{\ast }y(t)+y^{\ast }{\sigma }_3^2+r_3{y}^{\ast }{x}_2(t)y(t)\\ & \le -(d_3-(r_{2}a{)}^2-{\sigma }_3^2)(y(t)-y^{\ast }{)}^2+\frac{x_2^{{}^{\ast }2}}{m_1^2}+\frac{(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}\\ & +\frac{(x_2(t-\tau )-x_2^{\ast }{)}^2}{2m_1^2}+r_3{x}_2^{\ast }{y}^{\ast }y(t)+r_3{y}^{\ast }{x}_2(t)y(t)+y^{\ast }{\sigma }_3^2,\\ \mathrm{L}{V}_4& =\frac{(x_2(t)-x_2^{\ast }{)}^2}{2m_1^2}-\frac{(x_2(t-\tau )-x_2^{\ast }{)}^2}{2m_1^2},\\ \mathrm{L}{V}_5& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\\ \mathrm{L}{V}_6& =\frac{r_{2}a{x}_2(t)y(t)}{m+m_{1}y(t)+m_2{x}_2(t)}-\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )},\\ \mathrm{L}{V}_7& =a_2(r_1{x}_1(t)y(t)-d_1{x}_1(t)-c{x}_1(t))+a_4(c{x}_1(t)-d_2{x}_2(t)).\end{array}$Therefore, we have (18) $\begin{array}{rl}\mathrm{d}V(x_1,x_2,y)& =\mathrm{L}V\mathrm{d}t+a_1{\sigma }_1(x_1(t)-x_1^{\ast })x_1(t)\mathrm{d}{B}_1(t)+{\sigma }_2(x_2(t)-x_2^{\ast })x_2(t)\mathrm{d}{B}_2(t)\\ & +a_3{\sigma }_3(y(t)-x_1^{\ast })y(t)\mathrm{d}{B}_3(t)+a_3{\sigma }_{3}y(t)\mathrm{d}{B}_3(t).\end{array}$(18) where $\begin{array}{rl}\mathrm{L}V& \le -a_1(d_1+c-{\sigma }_1^2)(x_1(t)-x_1^{\ast }{)}^2+a_1{r}_1{x}_2^{\ast }{y}^{\ast }y(t)+a_1{r}_1{y}^{\ast }{x}_2(t)y(t)+a_1{\sigma }_1^2{x}_1^{\ast }\\ & -(d_2-\frac{c}{2}-{\sigma }_2^2)(x_2(t)-x_2^{\ast }{)}^2+\frac{c}{2}(x_1(t)-x_1^{\ast }{)}^2+{\sigma }_2^2{x}_2^{\ast }-a_3(d_3-(r_{2}a{)}^2-{\sigma }_3^2)\\ & \cdot (y(t)-y^{\ast }{)}^2+\frac{a_3{x}_2^{{}^{\ast }2}}{m_1^2}+\frac{a_3(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}+\frac{a_3(x_2(t-\tau )-x_2^{\ast }{)}^2}{2m_1^2}\\ & +\frac{a_3(x_2(t)-x_2^{\ast }{)}^2}{2m_1^2}-\frac{a_3(x_2(t-\tau )-x_2^{\ast }{)}^2}{2m_1^2}+a_3{r}_3{x}_2^{\ast }{y}^{\ast }y(t)+a_3{r}_3{y}^{\ast }{x}_2(t)y(t)\\ & +a_3{y}^{\ast }{\sigma }_3^2+\frac{a_5{r}_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}\\ & -a_5{d}_{3}y(t)-a_5{r}_3{x}_2(t)y(t)+a_4(c{x}_1(t)-d_2{x}_2(t))\\ & +\frac{a_5{r}_{2}a{x}_2(t)y(t)}{m+m_{1}y(t)+m_2{x}_2(t)}-\frac{a_5{r}_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}\\ & +a_2(r_1{x}_1(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\\ & \le -a_1(d_1+\frac{c}{2}-{\sigma }_1^2)(x_1(t)-x_1^{\ast }{)}^2-(d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2)(x_2(t)-x_2^{\ast }{)}^2\\ & -a_3(d_3-(r_{2}a{)}^2-{\sigma }_3^2)(y(t)-y^{\ast }{)}^2+(a_1{r}_1{x}_2^{\ast }{y}^{\ast }+a_3{r}_3{x}_2^{\ast }{y}^{\ast }-a_5{d}_3)y(t)\\ & +(a_1{r}_1{y}^{\ast }+a_3{r}_3{y}^{\ast }+a_2{r}_1-a_5{r}_3)x_2(t)y(t)+a_1{\sigma }_1^2{x}_1^{\ast }+{\sigma }_2^2{x}_2^{\ast }+a_3{y}^{\ast }{\sigma }_3^2\\ & +\frac{a_3{x}_2^{{}^{\ast }2}}{m_1^2}+\frac{a_3(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}+(\frac{a_5{r}_{2}a}{m}-a_4{d}_2)x_2(t)\\ & +(-a_2(d_1+c)+a_{4}c)x_1(t).\end{array}$Choosing $\begin{array}{rl}{a}_1& \le \frac{-a_2{r}_1+a_5{r}_3}{2r_1{y}^{\ast }},\\ & max\{\frac{a_5{r}_{2}ac}{d_2{m}_1(d_1+c)},\frac{a_5{r}_3{x}_2^{\ast }-a_5{d}_3}{r_1{x}_2^{\ast }}\}\le a_2\le \frac{a_5{r}_3}{r_1},\\ a_3& \le min\{2(d_2-\frac{c}{2}-{\sigma }_2^2)m_1^2,\frac{-a_2{r}_1+a_5{r}_3}{2r_3{y}^{\ast }}\},\\ & \frac{a_5{r}_{2}a}{d_2{m}_1}\le a_4\le \frac{a_2(d_1+c)}{c},\\ a_5& =1.\end{array}$then details calculation process of $a_1-a_5$ is given in Appendix, then we have $\begin{array}{rl}& d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2\ge 0,\\ & a_1{r}_1{x}_2^{\ast }{y}^{\ast }+a_3{r}_3{x}_2^{\ast }{y}^{\ast }-a_5{d}_3\le 0,\\ & \frac{a_5{r}_{2}a}{m_1}-a_4{d}_2\le 0,\\ & -a_2(d_1+c)+a_{4}c\le 0,\\ & a_1{r}_1{y}^{\ast }+a_3{r}_3{y}^{\ast }+a_2{r}_1-a_5{r}_3\le 0.\end{array}$We can obtain (19) $\begin{array}{rl}\mathrm{L}V& \le -a_1(d_1+\frac{c}{2}-{\sigma }_1^2)(x_1(t)-x_1^{\ast }{)}^2-(d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2)(x_2(t)-x_2^{\ast }{)}^2\\ & -a_3(d_3-(r_{2}a{)}^2-{\sigma }_3^2)(y(t)-y^{\ast }{)}^2+a_1{\sigma }_1^2{x}_1^{\ast }+{\sigma }_2^2{x}_2^{\ast }+a_3{y}^{\ast }{\sigma }_3^2+\frac{a_3{x}_2^{{}^{\ast }2}}{m_1^2}\\ & +\frac{a_3(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}.\end{array}$(19) Integrating both sides of (18(18) $\begin{array}{rl}\mathrm{d}V(x_1,x_2,y)& =\mathrm{L}V\mathrm{d}t+a_1{\sigma }_1(x_1(t)-x_1^{\ast })x_1(t)\mathrm{d}{B}_1(t)+{\sigma }_2(x_2(t)-x_2^{\ast })x_2(t)\mathrm{d}{B}_2(t)\\ & +a_3{\sigma }_3(y(t)-x_1^{\ast })y(t)\mathrm{d}{B}_3(t)+a_3{\sigma }_{3}y(t)\mathrm{d}{B}_3(t).\end{array}$(18) ) from 0 to t and taking the expectation, we get (20) $\begin{array}{rl}EV(t)-EV(0)& \le -E{\int }_0^t{a}_1(d_1+\frac{c}{2}-{\sigma }_1^2)(x_1(t)-x_1^{\ast }{)}^2\mathrm{d}\theta \\ & -E{\int }_0^t(d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2\phantom{\frac{c}{2}})(x_2(t)-x_2^{\ast }{)}^2\mathrm{d}\theta \\ & -E{\int }_0^t{a}_3(d_3-(r_{2}a{)}^2-{\sigma }_3^2)(y(t)-y^{\ast }{)}^2\mathrm{d}\theta \\ & +[a_1{\sigma }_1^2{x}_1^{\ast }+{\sigma }_2^2{x}_2^{\ast }+a_3{y}^{\ast }{\sigma }_3^2+\frac{a_3{r}_{2}a{x}_2^{{}^{\ast }2}}{m^2}+\frac{a_3{r}_{2}a(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}]t.\end{array}$(20) Dividing both sides by t and taking the limit superior, we have (21) $\begin{array}{rl}& \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t{a}_1(d_1+\frac{c}{2}-{\sigma }_1^2)(x_1(t)-x_1^{\ast }{)}^2\mathrm{d}\theta \\ & +\underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t(d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2)\\ & \cdot (x_2(t)-x_2^{\ast }{)}^2\mathrm{d}\theta +\underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t{a}_3(d_3-(r_{2}a{)}^2-{\sigma }_3^2)(y(t)-y^{\ast }{)}^2\mathrm{d}\theta \\ & \le a_1{\sigma }_1^2{x}_1^{\ast }+{\sigma }_2^2{x}_2^{\ast }+a_3{y}^{\ast }{\sigma }_3^2+\frac{a_3{x}_2^{{}^{\ast }2}}{m^2}+\frac{a_3(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}.\end{array}$(21) Obviously, (22) $\underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t(x_1(t)-x_1^{\ast }{)}^2+(x_2(t)-x_2^{\ast }{)}^2+(y(t)-y^{\ast }{)}^2\mathrm{d}\theta \le \frac{q}{p},$(22) where (23) $\begin{array}{rl}p& =min\{a_1(d_1+\frac{c}{2}-{\sigma }_1^2),(d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2),a_3(d_3-(r_{2}a{)}^2-{\sigma }_3^2)\},\\ q& =a_1{\sigma }_1^2{x}_1^{\ast }+{\sigma }_2^2{x}_2^{\ast }+a_3{y}^{\ast }{\sigma }_3^2+\frac{a_3{x}_2^{{}^{\ast }2}}{m^2}+\frac{a_3(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}.\end{array}$(23) We have completed the proof.

Remark 3.1

Theorem 3.1 shows that if the initial condition (5(5) $x_i(\theta )={\varphi }_i(\theta ),y(\theta )=\psi (\theta ),\theta \in [-\tau ,0],i=1,2,$(5) ) holds, the solution of system (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) oscillates around the equilibrium point $(x_1^{\ast },x_2^{\ast },y^{\ast })$ of system (3(3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) ) under certain conditions, and the upper bound of the solution of system (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) is positively correlated with the intensity of ${\sigma }_i^2,(i=1,2,3)$.

4. Persistence

In nature, whether the population has the conditions for persistence is an issue of great concern to us. Before discussing the persistence of stochastic system, we give the following assumption: (24) $\begin{array}{rl}{H}_2:\mu & =max\{{\sigma }_1,{\sigma }_2,{\sigma }_3\}<min\{x_1^{\ast }\sqrt{\frac{p}{2q_0}},x_2^{\ast }\sqrt{\frac{p}{2q_0}},y^{\ast }\sqrt{\frac{p}{2q_0}}\},\\ & a_3<min\{x_1^{\ast 2}\frac{p}{2q_1},x_2^{\ast 2}\frac{p}{2q_1},y^{\ast 2}\frac{p}{2q_1}\},\end{array}$(24) where $q_0=a_1{x}_1^{\ast }+x_2^{\ast }+a_3{y}^{\ast }$, $q_1=\frac{x_2^{{}^{\ast }2}}{m^2}+\frac{(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}$.

Theorem 4.1

For any given initial condition (5(5) $x_i(\theta )={\varphi }_i(\theta ),y(\theta )=\psi (\theta ),\theta \in [-\tau ,0],i=1,2,$(5) ), if assumptions $(H_0)$, $(H_1)$ and $(H_2)$ hold, the system (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) is persistent, that is (25) $\begin{array}{rl}& \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t{x}_1(\theta )\mathrm{d}\theta >0,\\ & \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t{x}_2(\theta )\mathrm{d}\theta >0,\\ & \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^{t}y(\theta )\mathrm{d}\theta >0.\end{array}$(25)

Proof.

According to (15(15) $\underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t[(x_1(\theta )-x_1^{\ast }{)}^2+(x_2(\theta )-x_2^{\ast }{)}^2+(y(\theta )-y^{\ast }{)}^2]\mathrm{d}\theta \le \frac{q}{p},$(15) ), we have (26) $\begin{array}{rl}& \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t(x_1(t)-x_1^{\ast }{)}^2\mathrm{d}\theta \le \frac{q}{p},\\ & \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t(x_2(t)-x_2^{\ast }{)}^2\mathrm{d}\theta \le \frac{q}{p},\\ & \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t(y(t)-y^{\ast }{)}^2\mathrm{d}\theta \le \frac{q}{p}.\end{array}$(26) As we know, $x_1(t)\ge 0$ and $x_1^{\ast }>0$, from $2x_1(t)x_1^{\ast }\ge (x_1^{\ast }{)}^2-(x_1(t)-x_1^{\ast }{)}^2$, we yield (27) $x_1(t)\ge \frac{x_1^{\ast }}{2}-\frac{(x_1(t)-x_1^{\ast }{)}^2}{2x_1^{\ast }}.$(27) By the condition of $\mu <x_1^{\ast }\sqrt{\frac{p}{2q_0}},a_3<x_1^{\ast 2}\frac{p}{2q_1}$, we have (28) $\begin{array}{rl}\underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t{x}_1(\theta )\mathrm{d}\theta & \ge \frac{x_1^{\ast }}{2}-\underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t\frac{(x_1(t)-x_1^{\ast }{)}^2}{2x_1^{\ast }}\mathrm{d}\theta \\ & \ge \frac{x_1^{\ast }}{2}-\frac{q}{2p{x}_1^{\ast }}\ge \frac{x_1^{\ast }}{2}-\frac{{\mu }^2{q}_0}{2p{x}_1^{\ast }}-\frac{a_3{q}_1}{2p{x}_1^{\ast }}>0.\end{array}$(28) Similarly, when $\mu <min\{x_2^{\ast }\sqrt{\frac{p}{2q_0}},y^{\ast }\sqrt{\frac{p}{2q_0}}\},a_3<min\{x_2^{\ast 2}\frac{p}{2q_1},y^{\ast 2}\frac{p}{2q_1}\}$, $\begin{array}{rl}& \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^t{x}_2(\theta )\mathrm{d}\theta >0,\\ & \underset{t\to \mathrm{\infty }}{lim sup}\frac{1}{t}E{\int }_0^{t}y(\theta )\mathrm{d}\theta >0.\end{array}$This completes the proof of Theorem 4.1.

5. Numerical simulations

In this section, we will verify theory results and discuss the effect of environmental noise. We use the Milstein's method to simulate the stochastic model (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) [9], the numerical scheme for stochastic model (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) is given by (29) $\{\begin{array}{rl}{x}_1(n+1)& =x_1(n)+\mathrm{\Delta }t(r_1{x}_2(n)y(n)-d_1{x}_1(n)-c{x}_1(n))+{\sigma }_1{x}_1(n)\sqrt{\mathrm{\Delta }t}\iota (k)\\ & +\frac{{\sigma }_1^2}{2}{x}_1(n)({\iota }^2(k)-1)\mathrm{\Delta }t,\\ x_2(n+1)& =x_2(n)+\mathrm{\Delta }t(c{x}_1(n)-d_2{x}_2(n))+{\sigma }_2{x}_2(n)\sqrt{\mathrm{\Delta }t}\eta (k)\\ & +\frac{{\sigma }_2^2}{2}{x}_2(n)({\eta }^2(k)-1)\mathrm{\Delta }t,\\ y(n+1)& =y(n)+\mathrm{\Delta }t(\frac{r_{2}a{x}_2(n+1-\frac{\tau }{\mathrm{\Delta }t})y(n+1-\frac{\tau }{\mathrm{\Delta }t})}{m+m_{1}y(n+1-\frac{\tau }{\mathrm{\Delta }t})+m_2{x}_2(n+1-\frac{\tau }{\mathrm{\Delta }t})}\\ & -d_{3}y(n)-r_3{x}_2(n)y(n))\\ & +{\sigma }_{3}y(n)\sqrt{\mathrm{\Delta }t}\xi (k)+\frac{{\sigma }_3^2}{2}y(n)({\xi }^2(k)-1)\mathrm{\Delta }t,\end{array}$(29) where $\iota (k),\eta (k),\xi (k)(k=1,2,\dots ,n)$ are the Gaussian random variables $N(0,1)$, and time increment $\mathrm{\Delta }t>0$.

We adopt the following parameter values: (30) $\begin{array}{rl}& r_1=0.96,r_2=0.25,r_3=0.3,d_1=0.2,d_2=0.4,d_3=0.1,\\ & m=0.1,m_1=0.2,m_2=0.01,a=0.45,c=0.65,\tau =4.\end{array}$(30) We fix the initial value $(x_1(0),x_2(0),y(0))=(2,1,3)$ except for the other specification, we obtain $E^{\ast }(9.813,15.946,0.545)$. Next we take into account the effect of environmental noise intensities on the system by comparing the following different situations. Use matlab to plot diagrams of those cases.

• In case 1, we take ${\sigma }_1=0.05,{\sigma }_2=0.05,{\sigma }_3=0.01$, we conclude that $\mu =0.01<min\{x_1^{\ast }\sqrt{\frac{p}{2q_0}},x_2^{\ast }\sqrt{\frac{p}{2q_0}},y^{\ast }\sqrt{\frac{p}{2q_0}}\}=min\{7.299,0.738,0.52\}$, and $a_3=0.00029<min\{0.00054,0.0014,0.0041\}$, it is obvious that the conditions of $H_0$, $H_1$, and $H_2$ in Theorem 4.1 are satisfied, the system is persistence. It is easy to see that the solution of stochastic model oscillate around the equilibrium point of the deterministic model and maintain persistent in Figure 1, when the white noise is small. To satisfy the conditions of Theorem 4.1, we should select the value of $a_3$ as smaller as possible, so there is no need to verify $a_3$ in the following cases because we just take it sufficient small.

• Then, to analyse the effect of ${\sigma }_1$, we change ${\sigma }_1=0.4$ and ${\sigma }_1=1.4$ in case 2 and case 3, respectively, the value of $\mu =0.4<min\{15.0358,24.4331,0.8348\}=0.8348$ in case 2, and $\mu =1.4>min\{15.0358,24.4331,0.8348\}=0.8348$ in case 3, we can conclude that case 2 satisfy Theorem 4.1 but case 3 do not. After comparing Figures 1 and 2, we find out that the change of ${\sigma }_1$ will affect the fluctuation of the whole system. More specifically, when ${\sigma }_1$ increases appropriately, the solution of the leaf-cutter ants and the fungus remain persistent below the equilibrium point of the deterministic model, and all three populations of the stochastic model are rapidly extinct when ${\sigma }_1$ is large enough, which indicates that a random variation of the death rate of immature leaf-cutter ants will turn the persistent population into extinction.

• Furthermore, we chose ${\sigma }_2=0.2$ and ${\sigma }_2=0.5$ in case 4 and case 5 to explore the effect of  ${\sigma }_2$ on the system (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ), $\mu =0.2<min\{10.5349,17.1181,0.5849\}=0.5849$ can satisfy Theorem 4.1, but $\mu =0.5<min\{8.4490,13.7297,0.4691\}=0.4691$ cannot satisfy Theorem 4.1. Which can be concluded by comparing Figures 1 and 3 is that the small changes of ${\sigma }_2$ will have a great impact on the whole system, and even lead to the extinction of the fungus (see Figure 3 c,d).

• At last, comparing Figures 1 and  4 that plotted by case 6 which ${\sigma }_3=0.4$ and case 7 which ${\sigma }_3=0.9$, we can make inquiry for the impact of ${\sigma }_3$ on system (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ), $\mu =0.4<min\{15.0358,24.4331,0.8348\}=0.8348$ can satisfy Theorem 4.1, but $\mu =0.9<min\{15.0358,24.4331,0.8348\}=0.8348$ cannot satisfy Theorem 4.1, we can conclude that ${\sigma }_3$ have the same impact as ${\sigma }_1$ on the system, that is to say when ${\sigma }_3$ increases appropriately, the solution of the leaf-cutter ants and the fungus remain persistent below the equilibrium point of the deterministic model, and all populations are rapidly extinct when ${\sigma }_3$ is large enough, which indicates that a random variation of the death rate of immature leaf-cutter ants will turn the persistent population into extinction.

6. Conclusion

In this paper, we propose a stochastic delay mutualistic model of the leaf-cutter ants with stage structure and their fungus garden. First of all, we analyse the existence and uniqueness of the positive solution, and consider the asymptotic property of the positive solution, then we give the sufficient conditions for the persistence of the population. By exploring some numerical simulations, we conclude the impact of the impact of environmental disturbances on the system. In nature, the secretion of leaf-cutter ants contains many unknown antibodies, so it is of great scientific significance to maintain the persistence of leaf-cutter ants and their fungus. The research results of this paper provide a good solution for preventing the extinction of species. That is, to artificially reduce the impact of environmental disturbances on the system, such as watering the fungus nursery during drought.

Figure 1. The evolution of $(x_1(t),x_2(t),y(t))$ for model (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) and its corresponding deterministic model (3(3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) ) with initial value $(x_1(0),x_2(0),y(0))=(2,1,3)$: (a) time series diagram of $x_1(t),x_2(t),y(t)$ in case 1: ${\sigma }_1=0.05,{\sigma }_2=0.05,{\sigma }_3=0.01$ and (b) phase diagram of $x_1(t),x_2(t),y(t)$ in case 1: ${\sigma }_1=0.04,{\sigma }_2=0.03,{\sigma }_3=0.02$.

Figure 2. The evolution of $(x_1(t),x_2(t),y(t))$ for model (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) and its corresponding deterministic model (3(3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) ) with initial value $(x_1(0),x_2(0),y(0))=(2,1,3)$: (a) time series diagram of $x_1(t),x_2(t),y(t)$ in case 2:${\sigma }_1=0.4,{\sigma }_2=0.05,{\sigma }_3=0.01$; (b) phase diagram of $x_1(t),x_2(t),y(t)$ in case 2:${\sigma }_1=0.4,{\sigma }_2=0.05,{\sigma }_3=0.01$; (c) time series diagram of $x_1(t),x_2(t),y(t)$ in case 3:${\sigma }_1=1.4,{\sigma }_2=0.05,{\sigma }_3=0.01$; (d) phase diagram of $x_1(t),x_2(t),y(t)$ in case 3:${\sigma }_1=1.4,{\sigma }_2=0.05,{\sigma }_3=0.01$.

Figure 3. The evolution of $(x_1(t),x_2(t),y(t))$ for model (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) and its corresponding deterministic model (3(3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) ) with initial value $(x_1(0),x_2(0),y(0))=(2,1,3)$: (a) time series diagram of $x_1(t),x_2(t),y(t)$ in case 4: ${\sigma }_1=0.05,{\sigma }_2=0.2,{\sigma }_3=0.01$; (b) phase diagram of $x_1(t),x_2(t),y(t)$ in case 4:${\sigma }_1=0.05,{\sigma }_2=0.2,{\sigma }_3=0.01$; (c) time series diagram of $x_1(t),x_2(t),y(t)$ in case 5: ${\sigma }_1=0.05,{\sigma }_2=0.5,{\sigma }_3=0.01$ and (d) phase diagram of $x_1(t),x_2(t),y(t)$ in case 5: ${\sigma }_1=0.05,{\sigma }_2=0.5,{\sigma }_3=0.01$.

Figure 4. The evolution of $(x_1(t),x_2(t),y(t))$ for model (4(4) $\{\begin{array}{rl}\mathrm{d}{x}_1(t)& =(r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t))\mathrm{d}t+{\sigma }_1{x}_1(t)\mathrm{d}{B}_1(t),\\ \mathrm{d}{x}_2(t)& =(c{x}_1(t)-d_2{x}_2(t))\mathrm{d}t+{\sigma }_2{x}_2(t)\mathrm{d}{B}_2(t),\\ \mathrm{d}y(t)& =[\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t)]\mathrm{d}t\\ & +{\sigma }_{3}y(t)\mathrm{d}{B}_3(t),\end{array}$(4) ) and its corresponding deterministic model (3(3) $\{\begin{array}{rl}\frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}& =r_1{x}_2(t)y(t)-d_1{x}_1(t)-c{x}_1(t),\\ \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}& =c{x}_1(t)-d_2{x}_2(t),\\ \frac{\mathrm{d}y(t)}{\mathrm{d}t}& =\frac{r_{2}a{x}_2(t-\tau )y(t-\tau )}{m+m_{1}y(t-\tau )+m_2{x}_2(t-\tau )}-d_{3}y(t)-r_3{x}_2(t)y(t),\end{array}$(3) ) with initial value $(x_1(0),x_2(0),y(0))=(2,1,3)$: (a) time series diagram of $x_1(t),x_2(t),y(t)$ in case 6:${\sigma }_1=0.05,{\sigma }_2=0.05,{\sigma }_3=0.4$; (b) phase diagram of $x_1(t),x_2(t),y(t)$ in case 6:${\sigma }_1=0.05,{\sigma }_2=0.05,{\sigma }_3=0.4$; (c) time series diagram of $x_1(t),x_2(t),y(t)$ in case 7:${\sigma }_1=0.05,{\sigma }_2=0.05,{\sigma }_3=0.9$; (d) phase diagram of $x_1(t),x_2(t),y(t)$ in case 7:${\sigma }_1=0.05,{\sigma }_2=0.05,{\sigma }_3=0.9$.

Acknowledgments

We are grateful to the editor and reviewers for their valuable comments and suggestions that greatly improved the presentation of this paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported by National Natural Sciences Foundation of China (11971405, 22072057).

Appendix

If Equation (19(19) $\begin{array}{rl}\mathrm{L}V& \le -a_1(d_1+\frac{c}{2}-{\sigma }_1^2)(x_1(t)-x_1^{\ast }{)}^2-(d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2)(x_2(t)-x_2^{\ast }{)}^2\\ & -a_3(d_3-(r_{2}a{)}^2-{\sigma }_3^2)(y(t)-y^{\ast }{)}^2+a_1{\sigma }_1^2{x}_1^{\ast }+{\sigma }_2^2{x}_2^{\ast }+a_3{y}^{\ast }{\sigma }_3^2+\frac{a_3{x}_2^{{}^{\ast }2}}{m_1^2}\\ & +\frac{a_3(x_2^{\ast }{y}^{\ast }{)}^2}{(m+m_1{y}^{\ast }+m_2{x}_2^{\ast }{)}^2}.\end{array}$(19) ) holds, the following equations should be satisfied: (A1) $\begin{array}{rl}& d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2\ge 0,\end{array}$(A1) (A2) $\begin{array}{rl}& a_1{r}_1{x}_2^{\ast }{y}^{\ast }+a_3{r}_3{x}_2^{\ast }{y}^{\ast }-a_5{d}_3\le 0,\end{array}$(A2) (A3) $\begin{array}{rl}& \frac{a_5{r}_{2}a}{m_1}-a_4{d}_2\le 0,\end{array}$(A3) (A4) $\begin{array}{rl}& -a_2(d_1+c)+a_{4}c\le 0,\end{array}$(A4) (A5) $\begin{array}{rl}& a_1{r}_1{y}^{\ast }+a_3{r}_3{y}^{\ast }+a_2{r}_1-a_5{r}_3\le 0.\end{array}$(A5) From Equation (A1(A1) $\begin{array}{rl}& d_2-\frac{c}{2}-\frac{a_3}{2m_1^2}-{\sigma }_2^2\ge 0,\end{array}$(A1) ), we can obtain $a_3\le 2(d_2-\frac{c}{2}-{\sigma }_2^2)m_1^2$.

It is easy to conclude that Equation (A5(A5) $\begin{array}{rl}& a_1{r}_1{y}^{\ast }+a_3{r}_3{y}^{\ast }+a_2{r}_1-a_5{r}_3\le 0.\end{array}$(A5) ) holds, if we have the following equations: $\begin{array}{rl}& a_2\le \frac{a_5{r}_3}{r_1},\\ & a_1\le \frac{-a_2{r}_1+a_5{r}_3}{2r_1{y}^{\ast }},\\ & a_3\le \frac{-a_2{r}_1+a_5{r}_3}{2r_3{y}^{\ast }}.\end{array}$Considering Equations (A2(A2) $\begin{array}{rl}& a_1{r}_1{x}_2^{\ast }{y}^{\ast }+a_3{r}_3{x}_2^{\ast }{y}^{\ast }-a_5{d}_3\le 0,\end{array}$(A2) ) and (A5(A5) $\begin{array}{rl}& a_1{r}_1{y}^{\ast }+a_3{r}_3{y}^{\ast }+a_2{r}_1-a_5{r}_3\le 0.\end{array}$(A5) ), if we have $\begin{array}{rl}{a}_2{r}_1{x}_2^{\ast }-a_5{r}_3{x}_2^{\ast }& \ge -a_5{d}_3,\\ a_2& \ge \frac{a_5{r}_3{x}_2^{\ast }-a_5{d}_3}{r_1{x}_2^{\ast }},\end{array}$then Equation (A2(A2) $\begin{array}{rl}& a_1{r}_1{x}_2^{\ast }{y}^{\ast }+a_3{r}_3{x}_2^{\ast }{y}^{\ast }-a_5{d}_3\le 0,\end{array}$(A2) ) holds.

From Equations (A4(A4) $\begin{array}{rl}& -a_2(d_1+c)+a_{4}c\le 0,\end{array}$(A4) ) and (A5(A5) $\begin{array}{rl}& a_1{r}_1{y}^{\ast }+a_3{r}_3{y}^{\ast }+a_2{r}_1-a_5{r}_3\le 0.\end{array}$(A5) ), $\frac{a_5{r}_{2}a}{d_2{m}_1}\le a_4\le \frac{a_2(d_1+c)}{c}$, at last we need $a_2\ge \frac{a_5{r}_{2}ac}{d_2{m}_1(d_1+c)}$, combine with $a_2\le \frac{a_5{r}_3}{r_1}$, so we need $\frac{r_3}{r_1}\ge \frac{r_{2}ac}{d_2{m}_1(d_1+c)}$, that is $r_1{r}_{2}ac-r_3{d}_2{m}_1(d_1+c)\le 0$.

In summary, we have $\begin{array}{rl}& a_1\le \frac{-a_2{r}_1+a_5{r}_3}{2r_1{y}^{\ast }},\\ & max\{\frac{a_5{r}_{2}ac}{d_2{m}_1(d_1+c)},\frac{a_5{r}_3{x}_2^{\ast }-a_5{d}_3}{r_1{x}_2^{\ast }}\}\le a_2\le \frac{a_5{r}_3}{r_1},\\ & a_3\le min\{2(d_2-\frac{c}{2}-{\sigma }_2^2)m_1^2,\frac{-a_2{r}_1+a_5{r}_3}{2r_3{y}^{\ast }}\},\\ & \frac{a_5{r}_{2}a}{d_2{m}_1}\le a_4\le \frac{a_2(d_1+c)}{c},\\ & a_5=1.\end{array}$