Survival and stationary distribution of a stochastic facultative mutualism

Abstract: To investigate the roles of both coupling noises and distributed delays with strong kernels, a novel delayed stochastic two-species facultative mutualism model is established, in where the strong kernels indicate that the maximum influence on the growth rate response at some time is due to population densities at the previous time, and the saturation effect is also incorporated because the facultative capacity of each species is finite and their interspecific mutualism should be upper bounded in real life. We first transfer the two-species stochastic model with strong kernels into an equivalent six-dimensional model through a linear chain technique. Later, sufficient conditions for the extinction exponentially, persistence in the mean, permanent in time average and stationary distribution are respectively obtained. Finally, numerical simulations are supplied to support our theoretical results. Our analytical results show that the coupling noise intensities play an important role in the long-time behaviors while the strong kernels are independent of the above asymptotic properties.


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] K(t) = t n σ n+1 e −σt n! , σ > 0, n = 0, 1, 2, · · · , 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 K(t) = σe −σt (n = 0, weak kernel), K(t) = tσ 2 e −σt (n = 1, strong kernel). (1.1) 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.

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 ) 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) = φ i (s) ≥ 0, s ∈ (−∞, 0] and φ i (0) > 0, where φ i are continuous bounded functions on (−∞, 0]. B i (t) are standard independent Brownian motions defined on a complete probability space (Ω, F , {F t } t≥0 , P) with a filtration {F t } t≥0 satisfying the usual conditions. And α 2 1i , α 2 2i 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 1+x 2 (t) − n 2 (t))dt.
(2.5) 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 (2.6) 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 t −∞ σ j e −σ 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 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 (2.7) 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. Signal Survival results of species are defined in Table 2. Table 2. Survival of species.
Cases Conditions

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).
We continue to give the following fundamental lemma on the global existence and uniqueness of positive solution to system (2.5).
Proof. Using Itô's formula, we obtain from system (2.5) that Applying Itô's formula to Eq (3.4) leads to It follows from Lemma 3.1 that m 1 ≤ 1 and m 2 ≤ 1, and then K is bounded when x 1 , x 2 ∈ (0, +∞). 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.
Dividing by t and using Lemma 3.1, one has The strong law of local martingales [47] states lim t→+∞ t −1 B i (t) = 0, and moreover, we derive from Eq (3.6) that lim sup Thus, both species are extinct exponentially by the assumptions of Theorem 3.1.
Proof. If a 1 + c 1 < ξ 1 , then we obtain from Theorem 3.1 that species x 1 will be extinct exponentially and lim sup t→+∞ 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 furthermore, it follows from Lemma 3.1 that lim t→+∞ m i (t)/t = 0 and lim t→+∞ n i (t)/t = 0. And Next, by Eq (3.7) one gets for arbitrarily small ε > 0, there is T > 0 such that for t ≥ T , which together with Eq (3.8) leads to We obtain from Eqs (3.5) and (3.9) that for t ≥ T , An application of Lemma 4 in [48] to the above two inequalities gives that is acquired by the arbitrariness of ε.
Proof. The proof is similar to that of Theorem 3.2, and hence we omit it.
Theorem 3.4. Assume that a 1 > ξ 1 and a 2 > ξ 2 , then both species will be permanent in time average 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.

Stationary distribution
In this section, to discuss the stationary distribution of system (2.5) we make some preliminaries. Consider the following integral equation ς l (s, X(s))dB l (s). (4.1) Lemma 4.1. [49]. Assume that the coefficients of Eq (4.1) are independent of t and satisfy the following conditions for a constant κ: in O R ⊂ R d + and there exists a nonnegative C 2 -function V(x) in R d + satisfying LV(x) ≤ −1 outside some compact set. Then Eq (4.1) exists a solution which has a stationary distribution. 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, Applying Itô's formula to Eq (4.2), one can derive that dV(X(t)) = LV(X(t))dt Obviously, we get, by Young's inequality, that It follows from Lemma 3.1 that m 1 ≤ 1 and m 2 ≤ 1. By Eq (4.3) we have For a constant η > 0, we have L(e ηt V(X(t))) = ηe ηt V(X(t)) + e ηt LV(X(t)) Choosing the above constant η small enough such that b i η/(2σ i ) − b i /2 < 0, and noting that b i η/(2σ i )n q+1 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 −ηt V(X(0)) and G 1 /η, implies that there exists a constant G 2 > 0 such that for all t ≥ 0 We further obtain from Eq (4.2) that E[x q i /q] ≤ E[V(X(t))] ≤ G 2 , and hence Also, it follows from Eq (4.2) that E[m q+1 i ] ≤ 2σ i G 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 The proof is complete.
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. 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 |ds.
Note that from which we conclude that Similarly, one has 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  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 W(X(x 1 , x 2 , m 1 , m 2 , n 1 , n 2 )) and a closed set U ⊂ R 6 + such that LW(X) ≤ −1 on X ∈ R 6 + \ U. Let q > 1 and W(X) = δ 1 (− ln x 1 + c 1 σ 2 n 2 ) + δ 2 (− ln x 2 + c 2 σ 1 n 1 ) where δ i = 2 λ i max{2, H i }, λ i = a i − ξ 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 > 0 sufficiently small such that , To prove that LW(X) ≤ −1 for X ∈ R 6 + \ U , we will consider the following six cases. Case 1. When X ∈ U 1 , one has from Eq (4.19) that where We have from δ 1 = 2 Then Similarly, for X ∈ U 2 and δ 2 = 2 λ 2 max{2, H 2 }, one has where To discuss the following Cases 2-6, we reconsider Eq (4.19) and obtain (4.20) where Case 2. When X ∈ U 3 , it follows from Eq (4.20) that LW( Case 3. When X ∈ U 5 and X ∈ U 6 , one has LW(X) When X ∈ U 9 and X ∈ U 10 , we get LW( Case 6. When X ∈ U 11 and X ∈ U 12 , then LW(X) ≤ H 3 + 5 − 1 From the above discussions we know the closed set U satisfying sup X∈R 6 + \U LW(X) ≤ −1.
Step 2: When b 1 b 2 − c 1 c 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.

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 ξ 1 = 0.5(α 2 11 + α 2 12 ) and ξ 2 = 0.5(α 2 21 + α 2 22 ) 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 < ξ i , ξ i < a i and ξ 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 ξ 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 1 b 2 − c 1 c 2 > 0 reveals that the effect of interspecific mutualism is less than that of intra-specific competition.