1.
Introduction
As a common relationship among species, mutualism has been extensively studied by many experts and scholars. Mutualism models have also received a lot of attention in population dynamics [1,2,3]. For example, the Lotka-Volterra mutualism model, the most common model of interspecific relationships, has the following form [4]
where xi(t) is the population size, ri is the intrinsic growth rate, aii>0 is the intraspecific competition coefficient, and aij>0(j≠i) is the effect of species j on species i. But, in the classical Lotka-Volterra mutualism model, the growth rate of each species is a linear function of the interacting species [5], which is unreasonable in real life. In order to describe the actual problem more accurately, Ayala and Gilpin et al. [5] proposed a nonlinear model in 1973
where θi denotes the positive parameter of the modified Lotka-Volterra mutualism model.
However, in nature, no species is deterministic and will be affected by various environmental factors. To describe these random perturbations in the environment, we consider that the growth rate ri of species in model (1.2) is linearly disturbed by Gaussian white noise [6,7,8,9]
For any time interval [0,t], let ˜ri(t) be the time average of ri(t). Then, we can get
where N(⋅,⋅) is the one-dimensional Gaussian distribution.
However, it is unreasonable to use a linear function of Gaussian white noise to simulate random perturbations in real life [10]. Obviously, the variance of the average growth rate ˜ri tends to ∞ at t→0+. This causes an unreasonable result that the stochastic fluctuations in the growth rate ri(t) can become very large in a small time interval [11]. Therefore, some scholars have begun to consider the use of mean-reverting Ornstein-Uhlenbeck process to simulate random perturbations, that is, the intrinsic growth rate ri of model (1.2) has the form [12,13]
where βi is the reversion rate, σi is the intensity of environmental fluctuation, ˉri is the mean recovery level, and βi,σi>0. The mean reversion of ri(t) to the constant level ˉri when βi>0 can be inferred from (1.3): if ri(t) has diffused above ˉri at some time, then the coefficient of the dt drift term is negative, so ri(t) will tend to move downwards immediately after, with the reverse holding if ri(t) is below ˉri at some time [14,15].
Further, we can get the solution of the OU process (1.3). First, by multiplying eβit on both sides of (1.3) and then sorting, we can get
Then,
Integrating from 0 to t on the both sides of above formula, we get
Thus, we have
where ri(0) is the initial value of the Ornstein-Uhlenbeck process ri(t). Then, we can get the expectation and variance of ri(t) as follows:
Thus, ri(t) obeys the Gaussian distribution N(ˉri+[ri(0)−ˉri]e−βit,σ2i2βi(1−e−2βit)), and σi∫t0e−βi(t−s)dBi(s) obeys the Gaussian distribution N(0,σ2i2βi(1−e−2βit)). From the mean of ri(t), it should be obvious to see the mean reversion feature: When ri(0) deviates from ˉri either upward or downward, the degree of deviation decays at the rate of e−βit and approaches ˉri. When t→+∞, the asymptotic mean and variance are ˉri and σ2i2βi, respectively, which can be understood as stationary, long-run equilibrium mean and variance.
But, in real life, in addition to small environmental disturbances such as white noise, there are also sudden environmental disturbances that cause significant changes in the survival status of species [16], such as earthquakes, hurricanes, epidemics, and so on [17,18]. These phenomena cannot be described by white noise, and the introduction of Lévy jumps in the basic model is a reasonable way to describe these phenomena [17,18]. So, we construct the following stochastic Gilpin-Ayala mutualism model driven by the mean-reverting OU process with Lévy jumps,
where xi(t−),i=1,2,⋯,n is the left limit of xi(t), modified parameter θi≥1,i=1,2,⋯,n, and Bi(t),i=1,2,⋯,n are independent standard Brownian motions defined on the probability space (Ω,F,{F}t≥0,P). N is a Poisson counting measure with characteristic measure v with v(Z)<∞, and Z is a measurable subset of (0,∞). ˜N represents a compensating random measure of Poisson random measure N, defined as ˜N(dt,dz)=N(dt,dz)−v(dz)dt. In order to satisfy the corresponding biological significance, we assume that for all z∈Z, the jump diffusion coefficients γi(z)>−1,i=1,2,⋯,n.
The model studied in this paper is improved on the basis of the classical Lotka-Volterra model, which no longer assumes linear exponential growth of the population and uses the mean reversion OU process to simulate small perturbations in the environment. This is a more reasonable method than assuming that the population parameters are linearly disturbed by Gaussian white noise. Furthermore, we also take into account the sudden disturbance of the population, so we introduce Lévy jumps to construct the model (1.5) studied in this paper. As far as we know, there are relatively few studies on such models, so it is very meaningful to study the properties of model (1.5).
For convenience, the following definitions are taken in this article:
For the sequence cij(1≤i,j≤n), we let
For a symmetric matrix A of order n, we define
2.
Existence and uniqueness of global solution
Assumption 2.1. For any k∈{1,2,...,n}, there exists a constant c>0, and the following inequalities hold:
Assumption 2.2. For matrix A=(0a12⋯a1na210⋯a2n⋮⋮⋮an1an2⋯0), there is
Remark 2.1. Assumption 2.1 indicates that the interference intensity of Lévy noise on the system should not be too large. Assumption 2.2 shows that although system (1.5) is a mutualism system, the intensity of intraspecific competition is still greater than the intensity of interactions between species. Otherwise, if the interference intensity of Lévy noise to the system is too large and the interaction intensity of species is greater than the intraspecific competition intensity, the solution of the system may explode in finite time.
Theorem 2.1. If Assumptions 2.1 and 2.2 hold, for any initial value (x(0),r(0))=(x1(0),⋯,xn(0),r1(0),⋯,rn(0))∈Rn+×Rn, there exists a unique solution (x(t),r(t))=(x1(t),⋯,xn(t),r1(t),⋯,rn(t)) of model (1.5) on t≥0, and it remains in Rn+×Rn with probability one.
Proof. Noting that all the coefficients of model (1.5) satisfy the local Lipschitz condition, for any initial value (x(0),r(0)), the system has a unique local solution (x(t),r(t)) on t∈[0,τe), where τe is the explosion time of the solution. Therefore, to prove the solution (x(t),r(t)) is global, it is needed to prove τe=∞ with probability one only. Hence, we take a sufficiently large p0>0 such that each component of (x(0),er(0)) falls within [1p0,p0]. For each integer p0 greater than p, we define the stopping time
Obviously, τp is monotonically increasing as p increases. For convenience, let τ∞=limp→∞τp, then τ∞≤τe holds with probability one. Therefore, if τ∞=∞, then τe=∞. In the following, we use proof by contradiction to prove τ∞=∞. Suppose τ∞=∞ does not hold with probability one, then there exist constants T>0 and ε∈(0,1) such that P(τ∞≤T)>ε. So, there exists p1≥p0 such that
Defining a C2-function V on Rn+×Rn
When xi>0, we have the inequality xi−1≥lnxi,1≤i≤n, so V is a nonnegative function.
Using the Itˆo formula, we can get
where
Then, there exists a constant N>0 such that
Substituting Eq (2.5) into (2.3), we have
Taking the integral from 0 to τp∧T on both sides of Eq (2.6) and taking the expectation, we obtain
When p≥p1, let Ωp={τp≤T}. From Eq (2.2), we can obtain P(Ωp)≥ε, and from the definition of τp, for each ω∈Ωp such that one of xi(τp,ω),eri(τp,ω)(i=1,2,⋯,n) is equal to p or 1p so that V(x(τp,ω),r(τp,ω)) is not less than (p−1−lnp),(1p−1+lnp), or 14(lnp)4, we have
According to Eq (2.7), we can get
where IΩp(ω) represents the indicator function of Ωp. Let p→∞. Then, ∞>V((x(0),r(0))+NT=∞, and thus we have a contradiction. Therefore, τ∞=∞ holds with probability one. Theorem 2.1 is proved.
3.
Moment boundedness of solution
Assumption 3.1. For any q>0, there is
Remark 3.1. Assumption 3.1 indicates that, in the mutualism system (1.5), for any species in the system, the intensity of intraspecific competition is greater than the sum of the weighted average of interspecific competition intensity, otherwise the system may not have a bounded qth moment.
Theorem 3.1. If Assumptions 2.1 and 3.1 hold, for any initial value (x(0),r(0))=(x1(0),⋯,xn(0),r1(0),⋯,rn(0))∈Rn+×Rn, the solution (x(t),r(t))=(x1(t),⋯,xn(t),r1(t),⋯,rn(t)) of model (1.5) has the property that
for any q>0, where κ(q) is a continuous function with respect to q. That is to say, the qth moment of the solution (x(t),r(t)) is bounded.
Proof. For any q≥2, defining a nonnegative C2-function V : Rn+×Rn→R+
Applying the Itˆo formula to the function V, we obtain
where
Then,
Let η=qmin{β1,β2⋯,βn}. Using the Itˆo formula again, we have
Integrating from 0 to t on both sides of Eq (3.2) and taking the expected value, we obtain
Combining this with Eq (3.1), we have
Substituting Eq (3.4) into (3.3), we get
Then,
Further,
This means E[xqi(t)]≤κ2(q),i=1,2,⋯,n,∀t≥0,q≥2. According to Hölder's inequality, for any ˜q∈(0,2), we obtain
Let κ(q)=max{κ2(q),(κ2(2))˜q2}. Then,
Theorem 3.1 is proved.
Remark 3.1. Similar to the proof of Theorem 3.1, we have E[ri(t)]2q≤Q(q),i=1,2,⋯,n,∀q>0.
4.
Existence of a stationary distribution
In this section, we give sufficient conditions for the existence of the stationary distribution of the solution of model (1.5), which reflects the persistence of species over long periods of time and is an important asymptotic property of population development. Many scholars have also studied the stability of the system. For example, Shao [19,20] studied the asymptotic stability in the distribution of stochastic predator-prey system with S-type distributed time delays, regime switching, and Lévy jumps, and also studied the stationary distribution of predator-prey models with Beddington-DeAngelis function response and multiple delays in a stochastic environment, and used different methods to analyze the stability of the systems according to the different disturbances on the models; Liu et al. [21] gave sufficient conditions for the distribution stability of a two-prey one-predator model with Lévy jumps. Before giving the theorem of the existence of stationary distributions, we give several lemmas.
Assumption 4.1. aii−n∑j=1,j≠iaji>0,βi>1,i=1,2,⋯,n.
Remark 4.1. Assumption 4.1 shows that the impact of intraspecific competition intensity on population density is greater than the sum of the growing-promoting effects of other species on the species, and the reversion rate of the intrinsic growth rate under the interference of OU processes should not be too small. Otherwise, the system may not have a stationary distribution.
Lemma 4.1. Let Xa(t)=(x1(t),⋯,xn(t),r1(t),⋯,rn(t)) and X˜a(t)=(˜x1(t),⋯,˜xn(t),˜r1(t),⋯,˜rn(t)) be solutions of model (1.5) with initial values of a=(x1(0),⋯,xn(0),r1(0),⋯,rn(0))∈D and ˜a=((˜x1(0),⋯,˜xn(0),˜r1(0),⋯,˜rn(0))∈D, where D is any compact subset of Rn+×Rn. If Assumptions 2.1 and 4.1 hold, then the following equation holds:
Proof. Defining a function W
Then, we obtain
Taking the integral on both sides of Eq (4.1) and taking the expectation, we obtain
Noting EW(t)≥0, we then have
Let θi=1,i=1,2,⋯,n. Then,
Thus, according Assumption 4.1, we have
Therefore, according (4.2), we get
Then, we have
According to model (1.5), there are
Therefore, E(xi(t)) and E(ri(t)),i=1,2,⋯,n, are continuously differentiable. According to Theorem 3.1 and Remark 3.1, we have
So, E(xi(t)),E(ri(t)),i=1,2,⋯,n, are uniformly continuous. According to the Barbalat lemma, it can be concluded that limt→+∞E|xi−˜xi|=0,limt→+∞E|ri−˜ri|=0,a.s., and therefore Lemma 4.1 is proven.
Here, in order to prove the following lemma, we introduce the following symbols. Define B(Rn+×Rn) as the set of all probability measures on Rn+×Rn, and for any two measures p1,p2∈B, define the metric dH as
where H={h:Rn+×Rn→R∣|h(x)−h(y)|≤|x−y|,|h(⋅)|≤1}.
Lemma 4.2. If Assumptions 2.1 and 4.1 hold, for any a∈Rn+×Rn, {p(t,a,⋅)∣t≥0} is the Cauchy sequence in the space B(Rn+×Rn) with metric dH.
Proof. For any fixed a∈Rn+×Rn, we only need to prove for any ε>0 that there is a T>0 such that
This is equivalent to prove
For any h∈H,t,s>0, we have
where ˉDR={a∈Rn+×Rn∣|a|≤R},ˉDcR=(Rn+×Rn)−ˉDR. According to Chebyshev's inequality, the transition probability {p(t,a,dz0∣t≥0)} is compact, i.e., for any ε>0, there exists a compact subset D=D(ε,a) over Rn+×Rn such that p(t,a,D)≥1−ε,∀t≥0, where R is sufficiently large and we have
According to Lemma 4.1, there exists T>0 such that
Substituting Eqs (4.5) and (4.6) into (4.4), we have
Since h is arbitrary, inequality (4.3) holds.
Lemma 4.3 [22]. Let M(t),t≥0, be a local martingale with initial value M(0)=0. If limt→+∞ρM(t)<∞, then limt→+∞M(t)t=0 where ρM(t)=∫t0d⟨M,M⟩(s)(1+s)2,t≥0, and ⟨M,M⟩(t) is the quadratic variational process of M(t).
Lemma 4.4. If Assumption 2.1 holds, the solutions of model (1.5) follow that
Proof. Defining a function W(t)=(n∑i=1xi(t))q=w(t)q,q≥1, using the Itˆo formula, we can get
Let θ>0 be sufffciently small and satisfy mθ≤t≤(m+1)θ,m=1,2,.... It follows that
where
Choose θ sufffciently small such that I<h(q). Therefore,
Let ε be an arbitrary positive constant. Based on Chebyshev's inequality, it follows that
By the Borel–Cantelli lemma, there exists an integer-valued random variable m0(ω) such that for almost all ω∈Ω, when m≥m0, we have
Hence, for almost all ω∈Ω, if m≥m0 and mθ≤t≤(m+1)θ, we have
Letting ε→0, we have
then,
Thus,
and it follows that
Lemma 4.5. If Assumption 2.1 holds, ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds>0,i=1,2,⋯,n, then populations xi(t) are weak persistent, a.s..
Proof. According to the definition of weak persistence, we need to prove lim supt→∞xi(t)>0,i=1,2,⋯,n. If the conclusion is not true, then P(U)>0, where U={ω:lim supt→∞xi(t,ω)=0,i=1,2,⋯,n}. Applying the Itˆo formula to lnxi(t) and integrating from 0 to t, we have
where
By Assumption 2.1,
From Lemma 4.3, we obtain
On the one hand, combining the strong law of large numbers [22] and the definition of the Ornstein–Uhlenbeck process, we have
If for all ω∈U, lim supt→∞xi(t,ω)=0,i=1,2,⋯,n, combining with Eq (4.9) we have
As this contradicts the assumption P(U)>0, then lim supt→∞xi(t)>0,i=1,2,⋯,n.
Theorem 4.1. If Assumptions 2.1 and 4.1 hold, ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds>0,i=1,2,⋯,n, and then model (1.5) has a unique ergodic stationary distribution.
Proof. To prove Theorem 4.1, first prove that there is a probability measure η(⋅)∈B such that for any a∈Rn+×Rn, the transition probability p(t,a,⋅) for Xa(t) converges weakly to η(⋅).
According to Proposition 2.5 [23], weak convergence of probability measures is the concept of a metric, i.e., p(t,a,⋅) weakly converging to η(⋅) is equivalent to the existence of a metric d such that limt→+∞d(p(t,a,⋅),η(⋅))=0.
So, we only need to prove that, for any a∈Rn+×Rn, there is
From Lemma 4.2, {p(t,0,⋅∣t≥0)} is the Cauchy sequence in the space B(Rn+×Rn) of the metric dH. So, there is a unique η(⋅)∈B such that
By Lemma 4.1 and the triangle inequality, we have
That is, the distribution of X(t) weakly converges to η.
By the Kolmogorov-Chapman equation, we know that η is constant. From Corollary 3.4.3 [24], it follows that η is strongly mixed. From Theorem 3.2.6 [24], we know that η is ergodic.
5.
Extinction
In this section, we give sufficient conditions for species extinction. For convenience, model (1.5) is written in matrix form as
where
Assumption 5.1. There exists a set of positive constants c1,c2,⋯,cn such that
holds, where C=diag(c1,c2,⋯,cn).
Remark 5.1. In Assumption 5.1, the introduction of the constant ci,i=1,2,⋯,n, indicates that the intraspecific competition intensity of the i-th population and the interspecific interaction intensity of the i-th population to the other n−1 species changes by ci times. If ci≥1, the intraspecific competition intensity and interspecific competition intensity increase by ci times; if ci<1, it is weakened by ci times. Assumption 5.1 means that, under the action of ci, the intraspecific competition intensity of each species is greater than the average of the action intensity of the species on other species and the action intensity of other species on the species. Otherwise, the population might not go extinct.
Theorem 5.1. If Assumptions 2.1 and 5.1 hold, for any initial value (x(0),r(0))∈Rn+×Rn, the solution (x(t),r(t)) of system (5.1) has the property that
In particular, if max1≤i≤nˉri+max1≤i≤n{aii(θi−1)θ−θiθi−1i}+lim supt→∞1t∫t0∫Zln(1+ˇγ(z))v(dz)ds<0, it implies limt→∞|x(t)|=0, and then x(t) is extinct, a.s..
Proof. Define a Lyapunov function
where c=(c1,c2,⋯,cn)T.
Applying the Itˆo formula, we can get
Using the Itˆo formula for lnV(x) again, we have
where
where we use the fact −aiixθii+aiixi≤aii(θi−1)θ−θiθi−1i,i=1,2,⋯,n, and
Substituting the above four inequalities into dlnV(x), we get
Integrating from 0 to t, we have
where
By Assumption 2.1,
From Lemma 4.3, we achieve
On the one hand, combining the strong law of large numbers [22] and the definition of the Ornstein–Uhlenbeck process, we have
Then,
According to Eq (5.2), we obtain
Taking the upper limit on both sides of Eq (5.3), we get
When max1≤i≤nˉri+max1≤i≤n{aii(θi−1)θ−θiθi−1i}+lim supt→∞1t∫t0∫Zln(1+ˇγ(z))v(dz)ds<0, it implies limt→∞|x(t)|=0, then x(t) is extinct, a.s.. Theorem 5.1 is proved.
Remark 5.1. Lemma 4.5 and Theorems 4.1 and 5.1 have very important biological explanations. From the theoretical results obtained, it can be seen that when ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds>0,i=1,2,⋯,n, population xi(t),i=1,2,⋯,n, will be weakly persistent, and if the parameters of model (1.5) satisfy the conditions of Assumption 4.1, the system has a stationary distribution, which indicates the persistence of population growth. When max1≤i≤nˉri+max1≤i≤n{aii(θi−1)θ−θiθi−1i}+lim supt→∞1t∫t0∫Zln(1+ˇγ(z))v(dz)ds<0, and the parameters of model (1.5) satisfy the conditions of Assumption 5.1, population x(t)=(x1(t),⋯,xn(t)) will be extinct. That is, for every 1≤i≤n, when ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds<−aii(θi−1)θ−θiθi−1i, population xi(t),i=1,2,⋯,n, will be extinct. So, the survival and extinction of the biological population of model (1.5) completely depend on the value of ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds.
Remark 5.2. In the following we analyze the effects of white noise simulated by the Ornstein-Uhlenbeck (OU) process on species survival and extinction. Since the OU process acts on the intrinsic growth rate ri,i=1,2,⋯,n, if model (1.5) is not affected by jump noise, the model takes the following form:
Using a similar method as above, it can be proved that when ˉri>0,i=1,2,⋯,n, populations xi(t),i=1,2,⋯,n, are weakly persistent; when max1≤i≤nˉri+max1≤i≤n{aii(θi−1)θ−θiθi−1i}<0, population x(t)=(x1(t),⋯,xn(t)) will be extinct. That is, when ˉri<−aii(θi−1)θ−θiθi−1i,i=1,2,⋯,n, populations xi(t),i=1,2,⋯,n, are extinct. Thus, when the system is only disturbed by OU process, the survival and extinction of the population is only related to the value of the average growth rate ˉri,i=1,2,⋯,n, of the population.
When ˉri>0,i=1,2,⋯,n, the species only disturbed by the OU process are weakly persistent. If the system is affected by jump noise and satisfies ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds<−aii(θi−1)θ−θiθi−1i,i=1,2,⋯,n, the species are extinct. When ˉri<−aii(θi−1)θ−θiθi−1i,i=1,2,⋯,n, the species that are only disturbed by the OU process are extinct, but if there are jump noises such that ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds>0,i=1,2,⋯,n, the species are weakly persistent. Therefore, it can be obtained that jump noise can make the survival system extinct and the extinction system survive.
Remark 5.3. In the following we analyze the effect of the jump diffusion coefficient γi(z),i=1,2,⋯,n, on population survival and extinction. If γi(z)<0,i=1,2,⋯,n, then lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds<0,i=1,2,⋯,n, means that jump noise could accelerate the extinction; if γi(z)>0,i=1,2,⋯,n, then lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds>0,i=1,2,⋯,n, means that jump noise is beneficial to the survival of the population.
6.
Computer simulations
In order to verify the above theoretical results on the stochastic Gilpin-Ayala mutualism model (1.5), we use the Euler-Maruyama method [25] and the R language, and select appropriate parameters for numerical verification. The combination of parameters is shown in Table 1, and the data is from [11,26,27,28,29]. Consider the following stochastic Gilpin-Ayala mutualism model for two populations:
Example 6.1. Letting v(Z)=1, and take the initial value of model (6.1) as x1(0)=0.11,x2(0)=0.2,r1(0)=0.2,r2(0)=0.1, choosing the combination A1 as the parameter values of model (6.1), and using the R language for numerical simulation, Figure 1 is obtained. By calculating, we have
Then, Assumption 2.2 is satisfied. According to Theorem 2.1, the global solution of the stochastic Gilpin-Ayala population model (6.1) exists.
The red lines in Figure 1(a), (b) represent the solutions of populations x1,x2 in a deterministic environment without any disturbance. It can be seen that the development trend of the population is a smooth curve, and the population will not explode due to the limitation of environmental resources. The blue lines in Figure 1(a), (b) show the variation trend of the populations x1,x2 whose growth rate is disturbed by the OU process. The green lines in Figure 1(a), (b) represent the global solution of the population under the disturbance of the OU process and Lévy noise, and since the jump noise values are both positive, it indicates that the jump noise plays a role in promoting the population growth. Combined with the figure, it can be found that, compared with the other two situations, the population number also increases significantly at the same time under the positive Lévy jump interference. Lévy jumps represent some disturbances in the environment that cause sudden changes in the survival condition of the population. For example, when t=16,t=22 in Figure 1(a), (b), we can also see that the population number changes suddenly, which indicates the effect of Lévy jumps on the population.
The red lines in Figure 1(c), (d) represent intrinsic growth rates r1,r2, while the blue lines in Figure 1(c), (d) represent population growth rates disturbed by the OU process, indicating that the interference of random environmental factors will make the growth rate r1(t),r2(t) fluctuate randomly under the interference of the OU process.
Example 6.2. Letting v(Z)=1, taking the initial value of model (6.1) as x1(0)=0.11,x2(0)=0.2,r1(0)=0.2,r2(0)=0.1, choosing the combination A2 as the parameter values of model (6.1), and using the R language for numerical simulation, Figure 2 is obtained. By calculating, we obtain
Then, Assumption 3.1 is satisfied. The numerical simulation results show that E(xq1),E(xq2) are less than κ(q), so E(xq1)≤κ(q),E(xq2)≤κ(q),q>0 hold, and Theorem 3.1 is verified.
From the biological point of view, since the environmental resources are limited, no biological population can grow indefinitely, so we hope that the system solution is ultimately bounded. In Figure 2, letting q=2, we have E(x21)≤κ(2),E(x22)≤κ(2), which indicates that the final second moment of the population is bounded, which conforms to the laws of survival in the real world.
Example 6.3. Letting v(Z)=1, taking the initial value of model (6.1) as x1(0)=0.11,x2(0)=0.2,r1(0)=0.2,r2(0)=0.2, choosing the combination A3 as the parameter values of model (6.1), and using the R language for numerical simulation, Figure 3 is obtained. By calculating, we get
Then, Assumption 4.1 and the conditions of weak persistent are satisfied. Figure 3(a), (c) represent the solution of x1(t),x2(t), and Figure 3(b), (d) represent the histogram of the solution of x1(t),x2(t). According Theorem 4.1, model (6.1) has a stationary distribution η(⋅).
As can be seen from Figure 3(a), (c), the values of population x1(t) are mostly between 1.5–3, and the values of population x2(t) of are mostly between 1.3–2.5, mainly concentrated in the middle region. Figure 3(b), (d) is the frequency histogram of populations x1(t),x2(t), shows a trend that high in the middle and low at both ends, and obeys normal distribution approximately. This indicates that if Assumption 4.1 and ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds>0,i=1,2, hold, the populations will continue to grow steadily over time, the population size will not change dramatically, and the different populations of the system will coexist harmoniously.
Example 6.4. Letting v(Z)=1, taking the initial value of model (6.1) as x1(0)=0.1,x2(0)=0.1,r1(0)=0.2,r2(0)=0.1, choosing the combination A4 as the parameter values of model (6.1), and using the R language for numerical simulation, Figure 4 is obtained.
According to the selected parameters, matrix A is (00.160.120), matrix S is (0.4000.5), and taking C=I∈R2×2, then
Futher,
and then Assumption 5.1 is satisfied. According to Theorem 5.1, the stochastic Gilpin-Ayala population model (6.1) is extinct.
According Remark 5.2, when max1≤i≤nˉri+max1≤i≤n{aii(θi−1)θ−θiθi−1i}<0, the populations x1(t),x2(t) are extinct when the populations disturbed only by the OU process. The red lines in Figure 4(a), (b) show the populations x1(t),x2(t) whose growth rate is disturbed by the OU process. When populations are disturbed only by the OU process, populations x1(t),x2(t) are extinct at t=20. The green lines in Figure 4(a), (b) represent the global solution of the population under the disturbance of the OU process and Lévy noise, population x1(t) is extinct at t=30 and population x2(t) is extinct at t=45. In this example, we let γ1(z)=0.1,γ2(z)=0.2, and according Remark 5.3, this indicates that when the Lévy noise value is greater than 0, the population growth is promoted, and the positive Lévy noise will delay the extinction of the population.
7.
Conclusions
In this paper, we study the dynamic behaviors of a stochastic Gilpin-Ayala mutualism model (1.5) driven by the mean-reverting OU process with Lévy jumps. The existence and uniqueness of the global solution, the moment boundedness of the solution, the existence of the stationary distribution and extinction of the stochastic Gilpin-Ayala mutualism model (1.5) are proved and verified by numerical examples. The existence and uniqueness of the global solution and the moment boundedness of the solution show that, the population shows a fluctuating growth trend under the interference of various random factors, and for any q>0, populations xi(t) (i=1,2,⋯,n) have bounded q-th moments. The existence of the stationary distribution and extinction of the solution show that when ˉri+lim supt→∞1t∫t0∫Zln(1+γi(z))v(dz)ds>0,i=1,2,⋯,n, model (1.5) has a stationary distribution η(⋅), which indicates the persistence of population growth, and the populations x(t) will be extinct when the conditions given by the assumption are satisfied.
However, in model (1.5), only the influence of the OU process and Lévy jumps on the survival of the population were considered. But, in the real world, there are many environmental factors that affect the population, such as rainfall, drought, seasonal changes, etc..These are the questions we will be working on in the future.
Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
This study was supported by National Natural Science Foundation of China (11401085), Heilongjiang Province Postdoctoral Funding Program (LBH-Q21059), Fundamental Research Projects of Chinese Central Universities (2572021DJ04).
Conflict of interest
The authors declare there is no conflict of interest.