A reliable numerical algorithm for fractional Lienard equation arising in oscillating circuits

1.   Introduction

According to the definition of ^[1], mutualism is the interaction of two/many species that benefits both/each other. As a common occurrence in nature, the mutualism interaction has an important impact and is well documented in many types of communities. Mutualism can be obligate or facultative, more specifically, an obligate mutualist is a species which requires the presence of another species for its survival ^[2] while a facultative mutualist is one which benefits in the same way from the association with another species but will be survive in its absence ^[3]. In recent years many mutualism models have been studied intensively and some good results have been obtained, for details to see stability and bifurcation ^{[4,5,6,7,8,9,10,11,12]}, persistence and extinction ^{[5,7,8,11,13,14,15,16,17,18,19,20]}, periodic solution and almost periodic solution ^{[20,21,22,23,24]}, optimal control ^[25,26], and stationary distribution ^[27,28].

Recalling many of the above studies, we can see that the distributed delay does't been taken into account. In fact the evolution of a species may reply on an average over past population or the cumulative effect of the past history, and distributed delays are often incorporated into populations models, for details to see references ^{[30,31,32,33,34,35,36,37,38,39,40]}. Particularly, the following Gamma distribution initially given by MacDonald ^[32]

is usually used for the delay kernels. It is well known that there exist two types of kernels: weak kernel and strong kernel, which respectively, represented by

Both weak kernel and strong kernel own different biological meanings: the former implies that the maximum weighted response of the growth rate is due to current population density while past densities have exponentially decreasing influence, and the latter indicate that the maximum influence on the growth rate response at some time is due to population density at the previous time (see ^[41]).

On the other hand, the deterministic models may be necessary to incorporate the environmental noises into these models. Nisbet and Gurney ^[42] and May^[43] suggested that the growth rates in population systems should own stochasticity and emerge random fluctuation to a certain degree. Thus, some noise sources were incorporated and then corresponding stochastic models were established. However, in many excellent investigations the authors assumed that one noise source only had an effect on the intrinsic growth rate of one species. Obviously, a reasonable idea is to consider that one noise source has influence not only on the intrinsic growth rate of one species but also on that of other species.

Inspired by the above arguments, in the next section we introduce a stochastic facultative mutualism model with distributed delays and strong kernels (see model (2.3)). Survival analysis and stationary distribution will become two topics of our whole research because the survival analysis reveals the persistence or extinction of one or more species in random environment and the stationary distribution is concerned with the stochastic statistical characteristic of the long-term behaviours of the sample trajectories. To the best of our knowledge, there are few published papers concerning model (2.3). The rest of this work is organized as follows. In Section 3, we present the main results including extinction exponentially, persistence in the mean, permanent in time average. In Section 4, we devote to investigating the existence and uniqueness of stationary distribution. In Section 5, numerical simulations are given to support our findings. A brief discussion on the biological meanings is shown in Section 6.

2.   Model and preliminaries

For the final export of the model we will discuss, let us first introduce the following facultative mutualism model with saturation effect which corresponds to a deterministic competitive model proposed by Gopalsamy ^[29]

where $x_i\ (i = 1, 2)$ are the densities of two species, $a_i > 0$ denote the intrinsic growth rates, $b_i > 0$ are the intraspecific competition rates, $c_i > 0$ are the interspecific mutualism rates and the nonlinear term $c_1x_2/(1+x_2)$ (or $c_2x_1/(1+x_1)$) reflects a saturation effect for large enough $x_2$ (or $x_1$).

With the idea of distributed delays and strong kernels, model (2.1) becomes a delayed version

Similar to ^[44], we use two coupling noise sources to model the random perturbations and derive a new stochastic model

with initial values $x_i(s) = \phi_i(s)\geq0$, $s\in(-\infty, 0]$ and $\phi_i(0) > 0$, where $\phi_i$ are continuous bounded functions on $(-\infty, 0]$. $B_i(t)$ are standard independent Brownian motions defined on a complete probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t\geq0}, P)$ with a filtration $\{\mathcal{F}_t\}_{t\geq0}$ satisfying the usual conditions. And $\alpha_{1i}^2$, $\alpha_{2i}^2$ denote the coupling noise intensities.

Assign

With the help of chain techniques, the delayed stochastic facultative mutualism model (2.3) is transformed into an equivalent undelayed stochastic six-dimensional system

with initial value $(x_1(0), x_2(0), m_1(0), m_2(0), n_1(0), n_2(0))$, where

To show the novelty of our work, we explicate the following two facts:

(I) Zuo et al. ^[36] recently investigated the following stochastic two-species cooperative model with distributed delays and weak kernels

Obviously, there exists a limitation in model (2.6): with the increase of one cooperator's density, its cooperative capacity will increase and tend to infinity (see $c_i\int_{-\infty}^t\sigma_je^{-\sigma_j(t-s)}x_j(s)ds$). But in real life, this interaction between different species should be upper-bounded (a similar argument can be seen in ^[45]). In our model (2.3), the interspecific mutualism terms $c_i\int_{-\infty}^t(t-s)\sigma_i^2e^{-\sigma_i(t-s)}\frac{x_j(s)}{1+x_j(s)}ds$ show saturation effects. Also, by the chain techniques, model (2.6) can be transformed into an equivalent four-dimensional system (see model (2.2) in ^[36]) whose dimension is lower than the above six-dimensional system (2.5). Finally, we must point out that there are two noise sources in model (2.6), but one noise source only has an effect on one species. It is easy to see that one noise source affects two species at the same time in model (2.3).

(II) In a recent investigation, Ning et al. ^[40] discussed a stochastic competitive model with distributed delays and weak kernels

Clearly, there exists a remarkably different mechanism of action between model (2.3) and model (2.7) because the former is interspecific mutualism (see the positive feedback parameters $c_i$) and the latter is interspecific competition (see the negative feedback parameters $-c_i$). The strong kernel functions of our model (2.3) differs distinctly from the weak kernel functions of the above model (2.7) (also see Eq (1.1)). By the chain techniques, an equivalent undelayed six-dimensional system to model (2.3) is also different from that of the undelayed four-dimensional system to model (2.7) (see model (3) in ^[40]).

As a continuation of previous work ^[40], our main purpose of this contribution is to investigate the effects of both coupling noise sources on the long-time behaviors of facultative mutualism model (2.3) with distributed delays and strong kernels by analyzing its equivalent system (2.5). For the convenience of the subsequent analysis, we list the following two definitions.

Definition 2.1. Signals and abbreviations are defined in Table 1.

Definition 2.2. Survival results of species are defined in Table 2.

3.   Survival analysis

This section is determined to analyze the survival of system (2.5). For the convenience of the subsequent discussion, we first estimate $m_i(t)$ and $n_i(t)$.

Lemma 3.1. $m_i(t)$, $n_i(t)\leq 1$ and $\lim_{t\rightarrow +\infty}{m_i(t)}/{t} = \lim_{t\rightarrow +\infty}{n_i(t)}/{t} = 0$, $i = 1, 2$.

Proof. We first consider the upper bound of $m_i(t)$. It follows from Eq (2.4) that

Also, we obtain from Eq (2.4) that

Obviously, Eqs (3.1) and (3.2) imply that $\lim_{t\rightarrow +\infty}{m_i(t)}/{t} = 0$, $\lim_{t\rightarrow +\infty}{n_i(t)}/{t} = 0$, $i = 1, 2$.

We continue to give the following fundamental lemma on the global existence and uniqueness of positive solution to system (2.5).

Lemma 3.2. For any initial value $X(0) = (x_1(0), x_2(0), m_1(0), m_2(0), n_1(0), n_2(0)) > 0$, system $(2.5)$ admits a unique global solution $X(t) = (x_1(t), x_2(t), m_1(t), m_2(t), n_1(t), n_2(t)) > 0$ for $t\geq 0$ a.s.

Proof. Using Itô's formula, we obtain from system (2.5) that

Define a $C^2$-function $U(X(t))$ by

Applying Itô's formula to Eq (3.4) leads to

where

and

It follows from Lemma 3.1 that $m_1\leq 1$ and $m_2\leq 1$, and then $K$ is bounded when $x_1, x_2\in (0, +\infty)$. As a consequence, $LU(X(t))$ is bounded. The rest proof is similar to that of Theorem 2.1 in ^[46], and hence we omit it.

Assign $\xi_1 = 0.5(\alpha_{11}^2+\alpha_{12}^2)$, $\xi_2 = 0.5(\alpha_{21}^2+\alpha_{22}^2)$. The following Theorems 3.1-3.4 focus on the survival results of both species.

Theorem 3.1. Both species are extinct exponentially if $a_1+c_1 < \xi_1$ and $a_2+c_2 < \xi_2$.

Proof. An integration of Eq (3.3) over $[0, t]$ leads to

Dividing by $t$ and using Lemma 3.1, one has

The strong law of local martingales ^[47] states $\lim_{t\rightarrow +\infty}t^{-1}B_i(t) = 0$, and moreover, we derive from Eq (3.6) that

Thus, both species are extinct exponentially by the assumptions of Theorem 3.1.

Theorem 3.2. Assume that $a_1+c_1 < \xi_1$ and $a_2 > \xi_2$, then species $x_1$ goes to exponential extinction while species $x_2$ is persistent in the mean and $\lim_{t\rightarrow +\infty}\langle x_2(t)\rangle = (a_2-\xi_2)/{b_2}\ a.s.$

Proof. If $a_1+c_1 < \xi_1$, then we obtain from Theorem 3.1 that species $x_1$ will be extinct exponentially and $\limsup_{t\rightarrow +\infty}t^{-1}\ln x_1(t) < 0$. Thus, we further derive that

An integration of the last four equations of system (2.5) on both sides results in

Consequently, we have

furthermore, it follows from Lemma 3.1 that $\lim_{t\rightarrow +\infty}{m_i(t)}/{t} = 0$ and $\lim_{t\rightarrow +\infty}{n_i(t)}/{t} = 0$. And $\lim_{t\rightarrow +\infty}{m_i(0)}/{t} = 0$, $\lim_{t\rightarrow +\infty}{n_i(0)}/{t} = 0$. So we have

Next, by Eq (3.7) one gets for arbitrarily small $\varepsilon > 0$, there is $T > 0$ such that for $t\geq T$,

which together with Eq (3.8) leads to

Let $-\varepsilon < t^{-1}\ln x_2(0) < \varepsilon/2$ for $t\geq T$. We obtain from Eqs (3.5) and (3.9) that for $t\geq T$,

An application of Lemma 4 in ^[48] to the above two inequalities gives that

Thus

is acquired by the arbitrariness of $\varepsilon$.

Theorem 3.3. Suppose that $a_1 > \xi_1$ and $a_2+c_2 < \xi_2$, then species $x_2$ goes to exponent extinction while species $x_1$ is persistent in the mean and $\lim_{t\rightarrow +\infty}\langle x_1(t)\rangle = (a_1-\xi_1)/{b_1}$ a.s.

Proof. The proof is similar to that of Theorem 3.2, and hence we omit it.

Theorem 3.4. Assume that $a_1 > \xi_1$ and $a_2 > \xi_2$, then both species will be permanent in time average and $(a_1-\xi_1)/{b_1}\leq\langle x_1\rangle_{*}\leq\langle x_1\rangle^{*}\leq(a_1-\xi_1+c_1)/{b_1}$, $(a_2-\xi_2)/{b_2}\leq\langle x_2\rangle_{*}\leq\langle x_2\rangle^{*}\leq(a_2-\xi_2+c_2)/{b_2}$ a.s.

Proof. Recalling Eq (3.5) and Lemma 3.1, we have

It follows from Lemma 4 in ^[48] that

So the desired conclusion is obtained.

4.   Stationary distribution

In this section, to discuss the stationary distribution of system (2.5) we make some preliminaries. Consider the following integral equation

Lemma 4.1. ^[49]. Assume that the coefficients of Eq $(4.1)$ are independent of $t$ and satisfy the following conditions for a constant $\kappa$:

in $O_\mathbb{R}\subset\mathbb{R}_+^d$ and there exists a nonnegative $C^2$-function $V(x)$ in $\mathbb{R}_+^d$ satisfying $LV(x)\leq-1$ outside some compact set. Then Eq $(4.1)$ exists a solution which has a stationary distribution.

Remark 4.1. ^[36]. The condition in Lemma 4.1 may be replaced by the global existence of the solution to Eq $(4.1)$ according to Remark 5 in Xu et al. ^[50].

We first give Lemma 4.2 which is important for the subsequent discussions.

Lemma 4.2. Assume that $X(t)$ is a solution to system $(2.5)$ with initial value $X(0) > 0$. Then there is a positive constant $Q_q$ such that for $q > 0$,

Proof. Let

Applying Itô's formula to Eq (4.2), one can derive that

in which

Obviously, we get, by Young's inequality, that

It follows from Lemma 3.1 that $m_1\leq1$ and $m_2\leq1$. By Eq (4.3) we have

For a constant $\eta > 0$, we have

Choosing the above constant $\eta$ small enough such that $b_i\eta/(2\sigma_i)-b_i/2 < 0$, and noting that $b_i\eta/(2\sigma_i)n_i^{q+1}\leq b_i\eta/(2\sigma_i)$ (see Lemma 3.1, $n_i\leq 1$), we further obtain

where

Integrating Eq (4.4) from 0 to $t$ and then taking the expectation, one has

which together with the continuity of $V(X(t))$ and the boundedness of $e^{-\eta t}V(X(0))$ and $G_1/\eta$, implies that there exists a constant $G_2 > 0$ such that for all $t\geq0$

We further obtain from Eq (4.2) that $\mathbb{E}[x_i^q/q]\leq \mathbb{E}[V(X(t))]\leq G_2$, and hence

Also, it follows from Eq (4.2) that $\mathbb{E}[m_i^{q+1}]\leq2\sigma_iG_2/b_i$. And by the Young's inequality, there exist $A_{1i} > 0$ such that

Similarly, we can prove there exist $A_{2i}$ such that

Let

then for $q > 0$,

The proof is complete.

Lemma 4.3. Suppose that $X(t)$ is a solution to system $(2.5)$ with $X(0) > 0$, then almost every path of $X(t)$ to system $(2.5)$ will be uniformly continuous.

Proof. First let us consider $x_1(t)$. For any $0\leq t_1\leq t_2$, an integration of the first equation of system (2.5) yields

Let $p > 2$, by the elementary inequality $|a+b+c|^p\leq3^{p-1}(|a^p+b^p+c^p|)$, one has

Recalling Lemma 4.2 (see Eq (4.5)) and using the Hölder inequality result in

In addition, using the Moment inequality and Lemma 4.2 leads to

Substituting Eqs (4.8) and (4.9) into Eq (4.7), we have

where

Next, we continue to consider $m_1(t)$. For any $0\leq t_1\leq t_2$, integrating the third equation of system (2.5) from $t_1$ to $t_2$ yields that

Similar to Eq (4.7), we have from the Hölder inequality and Lemma 4.2 that

where $F_2 = 2^p(t_2-t_1)^{\frac{p}{2}}\sigma_1^pQ_{p}$.

Finally, we investigate $n_1(t)$. For any $0\leq t_1\leq t_2$, by integrating the fifth equation of system (2.5) one has

Similar to Eq (4.11), one obtains

Repeating the same analysis method as above, we obtain that $x_2$, $m_2$ and $n_2$ own similar results as those of Eqs (4.10)-(4.12), respectively. Thus, it follows from Lemma 3.4 in ^[50] that almost every sample path of $X(t)$ is uniformly continuous.

Lemma 4.4. ^[51]. Let $h(t)$ be a nonnegative function defined on $[0, +\infty)$ such that $h(t)$ is integrable on $[0, +\infty)$ and is uniformly continuous on $[0, +\infty)$. Then $\lim_{t\rightarrow +\infty}h(t) = 0$.

Lemma 4.5. If $b_1b_2-c_1c_2 > 0$, then the solution $X(t) = (x_1(t), x_2(t), m_1(t), m_2(t), n_1(t)$, $n_2(t)) > 0$ to system $(2.5)$ is globally attractive, that is, for any solution $\bar{X}(t) = (\bar{x}_1(t), \bar{x}_2(t),$ $\bar{m}_1(t), \bar{m}_2(t), \bar{n}_1(t), \bar{n}_2(t))$ to system $(2.5)$ with $\bar{X}(0) > 0$, there exist $\lim_{t\rightarrow +\infty}|x_i(t)-\bar{x}_i(t)| = 0$, $\lim_{t\rightarrow +\infty}|m_i(t)-\bar{m}_i(t)| = 0$, $\lim_{t\rightarrow +\infty}|n_i(t)-\bar{n}_i(t)| = 0$ a.s., $i = 1, 2$.

Proof. It follows from the last two equations of system (2.5) that

Integrating both sides of Eq (4.13) over the interval $[0, t]$ yields that

As a consequence, one has

Note that

we have

from which we conclude that

Similarly, one has

Assign

A direct calculation of the right differential $D^+M(t)$ of $M(t)$ leads to

from which and Eqs (4.15) and (4.16) one can obtain

Rearranging the above inequality leads to

from which one gets $|x_i(t)-\bar{x}_i(t)|\in L^1[0, +\infty)$. Similarly, it follows from Eqs (4.15) and (4.16) that $|m_i(t)-\bar{m}_i(t)|, |n_i(t)-\bar{n}_i(t)|\in L^1[0, +\infty)$. Thus, we obtain from Lemmas 4.3 and 4.4 that

which confirms Lemma 4.5.

Theorem 4.1. If $a_1-\xi_1-c_1 > 0$, $a_2-\xi_2-c_2 > 0$ and $b_1b_2-c_1c_2 > 0$, then system $(2.5)$ admits a unique stationary distribution.

Proof. To finish this proof, we will consider the following two steps.

Step 1: We first prove the existence of stationary Markov process.

It follows from Lemma 4.1 and Remark 4.1 that we only need to find a nonnegative $C^2$-function $\mathcal{W}(X(x_1, x_2, m_1, m_2, n_1, n_2))$ and a closed set $\mathcal{U}\subset\mathbb{R}_+^6$ such that $L\mathcal{W}(X)\leq-1$ on $X\in\mathbb{R}_+^6\setminus\mathcal{U}$. Let $q > 1$ and

where $\delta_i = \frac{2}{\lambda_i}\max\{2, H_i\}, \lambda_i = a_i-\xi_i-c_i, i = 1, 2$, and the constants $H_i > 0$ will be given later. An application of Itô's formula and Lemma 3.1 gives

Choose $\epsilon > 0$ sufficiently small such that

$i = 1, 2$, where the constant $H_3 > 0$ is supplied later. A bounded closed set is defined by

Assign

To prove that $L\mathcal{W}(X)\leq-1$ for $X\in\mathbb{R}_+^6\setminus\mathcal{U}^\epsilon$, we will consider the following six cases.

Case 1. When $X\in\mathcal{U}_1^\epsilon$, one has from Eq (4.19) that

where

We have from $\delta_1 = \frac{2}{\lambda_1}\max\{2, H_1\}$ that $\delta_1\lambda_1/4\geq 1$. Then

Similarly, for $X\in\mathcal{U}_2^\epsilon$ and $\delta_2 = \frac{2}{\lambda_2}\max\{2, H_2\}$, one has

where

To discuss the following Cases 2-6, we reconsider Eq (4.19) and obtain

where

Case 2. When $X\in\mathcal{U}_3^\epsilon$, it follows from Eq (4.20) that $L\mathcal{W}(X)\leq H_3+5-\frac{1}{2}b_1\epsilon^{-(q+1)}\leq-1$. Similarly, if $X\in\mathcal{U}_4^\epsilon$, then $L\mathcal{W}(X)\leq H_3+5-\frac{1}{2}b_2\epsilon^{-(q+1)}\leq-1$.

Case 3. When $X\in\mathcal{U}_5^\epsilon$ and $X\in\mathcal{U}_6^\epsilon$, one has $L\mathcal{W}(X)\leq H_3+5-\frac{m_i^2}{4} < H_3+5-\frac{1}{4\epsilon^4}\leq-1$.

Case 4. When $X\in\mathcal{U}_7^\epsilon$ and $X\in\mathcal{U}_8^\epsilon$, then $L\mathcal{W}(X)\leq H_3+5-\frac{3n_i^2}{4} < H_3+5-\frac{3}{4\epsilon^4}\leq-1$.

Case 5. When $X\in\mathcal{U}_9^\epsilon$ and $X\in\mathcal{U}_{10}^\epsilon$, we get $L\mathcal{W}(X)\leq H_3+5-\frac{1}{2m_i}n_i < H_3+5-\frac{1}{2\epsilon^3}\epsilon^2\leq-1$.

Case 6. When $X\in\mathcal{U}_{11}^\epsilon$ and $X\in\mathcal{U}_{12}^\epsilon$, then $L\mathcal{W}(X)\leq H_3+5-\frac{1}{n_i}\frac{x_i}{1+x_i} < H_3+5-\frac{1}{\epsilon^2}\frac{\epsilon}{1+\epsilon}\leq-1$.

From the above discussions we know the closed set $\mathcal{U}^\epsilon$ satisfying $\sup_{X\in\mathbb{R}_+^6\setminus\mathcal{U}^\epsilon}L\mathcal{W}(X)\leq -1$.

Step 2: When $b_1b_2-c_1c_2 > 0$, we know from Lemma 4.5 that the solution $X(t)$ is globally attractive.

Combining Step 1 and Step 2, we complete the proof of Theorem 4.1.

5.   Numerical simulations

In this section, we will employ several specific examples to simulate the solutions to system (2.5), and verify the analytical results of the previous section.

For system (2.5), we first fix the parameter values as follows: $a_1 = 0.295$, $a_2 = 0.3$, $b_1 = 0.75$, $b_2 = 0.65$, $c_1 = 0.05$, $c_2 = 0.05$, $\sigma_1 = 0.1$, $\sigma_2 = 0.2$ and initial value $x_1(0) = 0.1$ and $x_2(0) = 0.12$. We will reveal how two coupling noise sources influence the long-time behaviors by choosing different noise intensities $\alpha_{11}^2$, $\alpha_{12}^2$, $\alpha_{21}^2$ and $\alpha_{22}^2$.

Choose $\alpha_{11}^2 = 0.65^2$, $\alpha_{12}^2 = 0.75^2$, $\alpha_{21}^2 = 0.64^2$, $\alpha_{22}^2 = 0.74^2$. A calculation gives that $a_1+c_1 = 0.3 < \xi_1 = 0.4925$ and $a_2+c_2 = 0.35 < \xi_2 = 0.4786$, which surely satisfies Theorem 3.1. Thus, both species will be EE (see Figure 1).

Choose $\alpha_{11}^2 = 0.65^2$, $\alpha_{12}^2 = 0.75^2$, $\alpha_{21}^2 = 0.03^2$, $\alpha_{22}^2 = 0.02^2$. By calculating we have $a_1+c_1 = 0.3 < \xi_1 = 0.4925$ and $a_2 = 0.3 > \xi_2 = 0.00065$, it then follows from Theorem 3.2 that species $x_1$ is EE while $x_2$ is PM (see Figure 2).

Choose $\alpha_{11}^2 = 0.03^2$, $\alpha_{12}^2 = 0.02^2$, $\alpha_{21}^2 = 0.65^2$, $\alpha_{22}^2 = 0.75^2$. Since $a_1 = 0.295 > \xi_1 = 0.00065$ and $a_2+c_2 = 0.35 < \xi_2 = 0.4925$, Theorem 3.3 indicates that species $x_1$ is PM while $x_2$ is EE (see Figure 3).

Choose $\alpha_{11}^2 = 0.03^2$, $\alpha_{12}^2 = 0.02^2$, $\alpha_{21}^2 = 0.02^2$, $\alpha_{22}^2 = 0.01^2$. Together with Theorem 3.4, we have from $a_1 = 0.295 > \xi_1 = 0.00065$ and $a_2 = 0.3 > \xi_2 = 0.00025$ that both species are PTA (see Figure 4).

Furthermore, taking the same noise intensities as in Figure 4, a calculation shows that $a_1-\xi_1-c_1 = 0.24435 > 0$, $a_2-\xi_2-c_2 = 0.24975 > 0$ and $b_1b_2-c_1c_2 = 0.485 > 0$. So we know, by Theorem 4.1, that system (2.5) owns a unique stationary distribution (see Figure 5).

6.   Discussion

This paper is concerned with a stochastic facultative mutualism model with saturation effect and distributed delays, in which strong kernel functions are incorporated (see model (2.3)). By analyzing a corresponding equivalent system (2.5), a set of easily verifiable sufficient conditions for the survival results and stationary distribution of system (2.5) is established. Note that $\xi_1 = 0.5({\alpha_{11}^2+\alpha_{12}^2})$ and $\xi_2 = 0.5({\alpha_{21}^2+\alpha_{22}^2})$ and from the above theoretical results, we have the following conclusions:

● Theorems 3.1-3.4 imply that large coupling noise intensities are harmful for the survival of both species while suitably small coupling noise intensities are advantage for them (see Figures 1-4).

● It follows from Theorems 3.2 and 3.3 that if the intrinsic growth rate of one species is small, the coupling noise intensities are large amplitude and the cooperation from the other species is not enough, then one species will be extinct exponentially (see Figures 2-3). However, if the other species only owns large intrinsic growth rate and relatively small coupling noise intensities, then it will be persistent in the mean (see Figures 2-3).

● Compared with the conditions of Theorems 3.1, 3.4 and 4.1, that is, $a_i+c_i < \xi_i$, $\xi_i < a_i$ and $\xi_i < a_i-c_i$ $(i = 1, 2)$, we find that the intrinsic growth rates $a_i$ are larger and larger while the coupling noise intensity $\xi_i$ is smaller and smaller, then both species go from exponent extinction (see Figure 1) to permanence in time average (see Figure 4) to the existence of a unique stationary distribution (see Figure 5). In addition, $b_1b_2-c_1c_2 > 0$ reveals that the effect of interspecific mutualism is less than that of intra-specific competition.

Acknowledgments

The authors thank the editor and referees for their careful reading and valuable comments. The work is supported by the College Innovation Team Project of Hubei Provincial Department of Education (No.T201812).

Conflict of interest

The authors declare there is no conflict of interest.