A stochastic Gilpin-Ayala mutualism model driven by mean-reverting OU process with Lévy jumps

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 $x_i(t)$ is the population size, $r_i$ is the intrinsic growth rate, $a_{ii} > 0$ is the intraspecific competition coefficient, and $a_{ij} > 0(j\ne 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 $\theta_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 $r_i$ of species in model (1.2) is linearly disturbed by Gaussian white noise ^[6,7,8,9]

For any time interval $[0, t]$, let $\tilde{r} _i(t)$ be the time average of $r_i(t)$. Then, we can get

where $\mathbb{N}(\cdot, \cdot)$ 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 $\tilde{r} _i$ tends to $\infty$ at $t\to 0^+$. This causes an unreasonable result that the stochastic fluctuations in the growth rate $r_i(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 $r_i$ of model (1.2) has the form ^[12,13]

where $\beta_i$ is the reversion rate, $\sigma_i$ is the intensity of environmental fluctuation, $\bar{r}_i$ is the mean recovery level, and $\beta_i, \sigma_i > 0$. The mean reversion of $r_i(t)$ to the constant level $\bar{r}_i$ when $\beta_i > 0$ can be inferred from (1.3): if $r_i(t)$ has diffused above $\bar{r}_i$ at some time, then the coefficient of the $\mathrm{d}t$ drift term is negative, so $r_i(t)$ will tend to move downwards immediately after, with the reverse holding if $r_i(t)$ is below $\bar{r}_i$ at some time ^[14,15].

Further, we can get the solution of the OU process (1.3). First, by multiplying $e^{\beta_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 $r_i(0)$ is the initial value of the Ornstein-Uhlenbeck process $r_i(t)$. Then, we can get the expectation and variance of $r_i(t)$ as follows:

Thus, $r_i(t)$ obeys the Gaussian distribution $N\left (\bar{r}_i+\left [ r_i(0) -\bar{r}_i\right ]e^{-\beta_i t}, \dfrac{\sigma_i ^2}{2\beta_i} \left (1-e^{-2\beta_i t} \right) \right)$, and $\sigma_i\int_{0}^{t} e^{-\beta_i(t-s)} \mathrm{d}B_i (s)$ obeys the Gaussian distribution $N\left (0, \dfrac{\sigma_i ^2}{2\beta_i} \left (1-e^{-2\beta_i t} \right) \right)$. From the mean of $r_i(t)$, it should be obvious to see the mean reversion feature: When $r_i(0)$ deviates from $\bar{r}_i$ either upward or downward, the degree of deviation decays at the rate of $e^{-\beta_i t}$ and approaches $\bar{r}_i$. When $t \to +\infty$, the asymptotic mean and variance are $\bar{r}_i$ and $\dfrac{\sigma_i ^2}{2\beta_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 $x_i(t^-), i = 1, 2, \cdots, n$ is the left limit of $x_i(t)$, modified parameter $\theta_i\ge1, i = 1, 2, \cdots, n$, and $B_i(t), i = 1, 2, \cdots, n$ are independent standard Brownian motions defined on the probability space $\left (\Omega, \mathcal{F}, \left \{ \mathcal{F} \right \} _{t\ge 0}, \mathbb{P} \right)$. $N$ is a Poisson counting measure with characteristic measure $v$ with $v\left (Z \right) < \infty$, and $Z$ is a measurable subset of $(0, \infty)$. $\tilde{N}$ represents a compensating random measure of Poisson random measure $N$, defined as $\tilde{N}(\mathrm{d}t, \mathrm{d}z) = N(\mathrm{d}t, \mathrm{d}z)-v\left (\mathrm{d}z \right) \mathrm{d}t$. In order to satisfy the corresponding biological significance, we assume that for all $z\in Z$, the jump diffusion coefficients $\gamma_i(z) > -1, i = 1, 2, \cdots, 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 $c_{ij}\; (1\le i, j \le 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 \in \{1, 2, ..., n\}$, there exists a constant $c > 0$, and the following inequalities hold:

Assumption 2.2. For matrix $A = \begin{pmatrix} 0 & a_{12} & \cdots & a_{1n}\\ a_{21} & 0 & \cdots & a_{2n}\\ \vdots & \vdots & & \vdots\\ a_{n1} & a_{n2} & \cdots & 0 \end{pmatrix}$, 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)) = (x_1(0), \cdots, x_n(0), r_1(0), \cdots, r_n(0))\in \mathbb{R}_+^n \times \mathbb{R}^n$, there exists a unique solution $(x(t), r(t)) = (x_1(t), \cdots, x_n(t), r_1(t), \cdots, r_n(t))$ of model $(1.5)$ on $t\geq0$, and it remains in $\mathbb{R}_+^n \times \mathbb{R}^n$ 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 \in \left [ 0, \tau _{e} \right)$, where $\tau _{e}$ is the explosion time of the solution. Therefore, to prove the solution $(x(t), r(t))$ is global, it is needed to prove $\tau _{e} = \infty$ with probability one only. Hence, we take a sufficiently large $p_0 > 0$ such that each component of $(x(0), e^{r(0)})$ falls within $[ \dfrac{1}{p_0}, p_0 ]$. For each integer $p_0$ greater than $p$, we define the stopping time

Obviously, $\tau _{p}$ is monotonically increasing as $p$ increases. For convenience, let $\tau _{\infty } = \lim\limits_{p \to \infty} \tau _{p}$, then $\tau _{\infty }\le \tau _{e }$ holds with probability one. Therefore, if $\tau _{\infty } = \infty$, then $\tau _{e } = \infty$. In the following, we use proof by contradiction to prove $\tau _{\infty } = \infty$. Suppose $\tau _{\infty } = \infty$ does not hold with probability one, then there exist constants $T > 0$ and $\varepsilon \in\left (0, 1 \right)$ such that $\mathbb{P} (\tau _{\infty}\le T) > \varepsilon$. So, there exists $p_1\ge p_0$ such that

Defining a $C^2$-function $V$ on $\mathbb{R}_+^n \times \mathbb{R}^n$

When $x_i > 0$, we have the inequality $x_i-1 \ge \ln{x_i}, 1\le i\le n$, so $V$ is a nonnegative function.

Using the $\mathrm{It\hat{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 $\tau _p\wedge T$ on both sides of Eq $(2.6)$ and taking the expectation, we obtain

When $p\ge p_1$, let $\Omega _p = \left \{ \tau _p\le T\right \}$. From Eq (2.2), we can obtain $\mathbb{P}\left (\Omega _p \right) \ge \varepsilon$, and from the definition of $\tau _p$, for each $\omega \in \Omega _p$ such that one of $x_i(\tau _p, \omega), e^{r_i\left (\tau _p, \omega\right)}(i = 1, 2, \cdots, n)$ is equal to $p$ or $\dfrac{1}{p}$ so that $V (x (\tau _p, \omega), r(\tau _p, \omega))$ is not less than $\left (p-1-\ln{p} \right), \left (\dfrac{1}{p} -1+ \ln{p}\right)$, or $\dfrac{1}{4}\left (\ln{p} \right) ^4$, we have

According to Eq (2.7), we can get

where $I_{\Omega _p}(\omega)$ represents the indicator function of ${\Omega _p}$. Let $p\to \infty$. Then, $\infty > V((x(0), r(0)) +NT = \infty$, and thus we have a contradiction. Therefore, $\tau _{\infty } = \infty$ 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 $q$th moment.

Theorem 3.1. If Assumptions 2.1 and 3.1 hold, for any initial value $(x(0), r(0)) = (x_1(0), \cdots, x_n(0), r_1(0), \cdots, r_n(0))\in \mathbb{R}_+^n \times \mathbb{R}^n$, the solution $(x(t), r(t)) = (x_1(t), \cdots, x_n(t), r_1(t), \cdots, r_n(t))$ of model (1.5) has the property that

for any $q > 0$, where $\kappa (q)$ is a continuous function with respect to $q$. That is to say, the $q$th moment of the solution $(x(t), r(t))$ is bounded.

Proof. For any $q\ge 2$, defining a nonnegative $C^2$-function $V$ : $\mathbb{R}_+^n \times \mathbb{R}^n\to \mathbb{R} _+$

Applying the $\mathrm{It\hat{o} }$ formula to the function $V$, we obtain

where

Then,

Let $\eta = q\min\left \{ \beta _1, \beta_2 \cdots, \beta_n\right \}$. Using the $\mathrm{It\hat{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 $\mathbb{E}\left [ x_i^q(t) \right ]\le \kappa_2 (q), i = 1, 2, \cdots, n, \forall t\ge 0, q\ge 2$. According to Hölder's inequality, for any $\tilde{q} \in (0, 2)$, we obtain

Let $\kappa(q) = \max \left \{ \kappa_2(q), \left (\kappa_2(2) \right) ^{\frac{\tilde{q} }{2} } \right \}$. Then,

Theorem 3.1 is proved.

Remark 3.1. Similar to the proof of Theorem 3.1, we have $\mathbb{E} \left [r_i (t)\right ] ^{2q}\le Q(q), i = 1, 2, \cdots, n, \forall 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. $a_{ii}- \sum\limits_{j = 1, j\ne i}^na_{ji} > 0, \beta_{i} > 1, i = 1, 2, \cdots, 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 $X^{a}(t) = (x_1(t), \cdots, x_n(t), r_1(t), \cdots, r_n(t))$ and $X^{\tilde{a} }(t) = (\tilde{x}_1(t), \cdots, \tilde{x}_n(t), \tilde{r}_1(t), \cdots, \tilde{r}_n(t))$ be solutions of model (1.5) with initial values of $a = (x_1(0), \cdots, x_n(0), r_1(0), \cdots, r_n(0)) \in D$ and $\tilde{a} = ((\tilde{x}_1(0), \cdots, \tilde{x}_n(0), \tilde{r}_1(0), \cdots, \tilde{r}_n(0)) \in D$, where $D$ is any compact subset of $\mathbb{R}^n_+ \times \mathbb{R}^n$. 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 $\mathbb{E} W(t) \ge 0$, we then have

Let $\theta_i = 1, i = 1, 2, \cdots, 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, $\mathbb{E}(x_i(t))$ and $\mathbb{E}(r_i(t)), i = 1, 2, \cdots, n,$ are continuously differentiable. According to Theorem 3.1 and Remark 3.1, we have

So, $\mathbb{E}(x_i(t)), \mathbb{E}(r_i(t)), i = 1, 2, \cdots, n,$ are uniformly continuous. According to the Barbalat lemma, it can be concluded that $\lim\limits_{t \to +\infty}\mathbb{E} \left | x_i-\tilde{x} _i \right | = 0, \lim\limits_{t \to +\infty}\mathbb{E} \left | r_i-\tilde{r}_i \right | = 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(\mathbb{R}_+^n \times \mathbb{R}^n)$ as the set of all probability measures on $\mathbb{R}_+^n \times \mathbb{R}^n$, and for any two measures $p_1, p_2\in B$, define the metric $d_H$ as

where $H = \left \{ h:\mathbb{R}_+^n \times \mathbb{R}^n \to\mathbb{R} \mid \left | h(x)-h(y) \right | \le \left | x-y \right |, \left | h(\cdot) \right | \le 1 \right \}.$

Lemma 4.2. If Assumptions 2.1 and 4.1 hold, for any $a\in \mathbb{R}_+^n \times \mathbb{R}^n$, $\left \{ p(t, a, \cdot) \mid t\ge 0 \right \}$ is the Cauchy sequence in the space $B(\mathbb{R}_+^n \times \mathbb{R}^n)$ with metric $d_H$.

Proof. For any fixed $a\in \mathbb{R}_+^n \times \mathbb{R}^n$, we only need to prove for any $\varepsilon > 0$ that there is a $T > 0$ such that

This is equivalent to prove

For any $h \in H, t, s > 0$, we have

where $\bar{D}_\mathbb{R} = \left \{ a\in \mathbb{R}_+^n \times \mathbb{R}^n \mid \left | a \right |\le R \right \}, \bar{D}_\mathbb{R} ^c = (\mathbb{R}_+^n \times \mathbb{R}^n)-\bar{D}_\mathbb{R}$. According to Chebyshev's inequality, the transition probability $\left \{ p(t, a, \mathrm{d}z_0 \mid t\ge 0) \right \}$ is compact, i.e., for any $\varepsilon > 0$, there exists a compact subset $D = D(\varepsilon, a)$ over $\mathbb{R}_+^n \times \mathbb{R}^n$ such that $p(t, a, D)\ge 1-\varepsilon, \forall t\ge 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\ge 0$, be a local martingale with initial value $M(0) = 0$. If $\lim\limits _{t \to +\infty} \rho _{M}(t) < \infty$, then $\lim\limits _{t \to +\infty}\dfrac{M(t)}{t} = 0$ where $\rho _{M}(t) = \int_{0}^{t} \dfrac{\mathrm{d}\left \langle M, M \right \rangle (s) }{(1+s)^2}, t\ge 0$, and $\left \langle M, M \right \rangle (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) = \left(\sum\limits_{i = 1}^n x_i(t)\right)^q = w(t)^q, q\ge 1$, using the $\mathrm{It\hat{o} }$ formula, we can get

Let $\theta > 0$ be sufffciently small and satisfy $m\theta \le t \le (m + 1)\theta, m = 1, 2, ...$. It follows that

where

Choose $\theta$ sufffciently small such that $I < h(q)$. Therefore,

Let $\varepsilon$ 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 $m_0(\omega)$ such that for almost all $\omega \in \Omega$, when $m \ge m_0$, we have

Hence, for almost all $\omega \in \Omega$, if $m \ge m_0$ and $m \theta \le t \le (m+1)\theta$, we have

Letting $\varepsilon \to 0$, we have

then,

Thus,

and it follows that

Lemma 4.5. If Assumption 2.1 holds, $\bar{r}_i+\limsup\limits_{t \to \infty}\dfrac{1}{t}\int_{0}^{t}\int_{Z}\ln(1+\gamma_i(z))v(\mathrm{d} z)\mathrm{d} s > 0, i = 1, 2, \cdots, n$, then populations $x_i(t)$ are weak persistent, a.s..

Proof. According to the definition of weak persistence, we need to prove $\limsup\limits_{t \to \infty}x_i(t) > 0, i = 1, 2, \cdots, n$. If the conclusion is not true, then $\mathbb{P} (U) > 0$, where $U = \left \{ \omega:\limsup\limits_{t \to \infty}x_i(t, \omega) = 0, i = 1, 2, \cdots, n\right \}$. Applying the $\mathrm{It\hat{o} }$ formula to $\ln x_i(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 $\omega \in U$, $\limsup\limits_{t \to \infty}x_i(t, \omega) = 0, i = 1, 2, \cdots, n$, combining with Eq (4.9) we have

As this contradicts the assumption $\mathbb{P} (U) > 0$, then $\limsup\limits_{t \to \infty}x_i(t) > 0, i = 1, 2, \cdots, n$.

Theorem 4.1. If Assumptions 2.1 and 4.1 hold, $\bar{r}_i +\limsup\limits_{t \to \infty}\dfrac{1}{t}\int_{0}^{t}\int_{Z}\ln(1+\gamma_i(z))v(\mathrm{d} z)\mathrm{d} s > 0, i = 1, 2, \cdots, 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 $\eta (\cdot)\in B$ such that for any $a\in \mathbb{R}_+^n \times \mathbb{R}^n$, the transition probability $p(t, a, \cdot)$ for $X^{a}(t)$ converges weakly to $\eta (\cdot)$.

According to Proposition 2.5 ^[23], weak convergence of probability measures is the concept of a metric, i.e., $p(t, a, \cdot)$ weakly converging to $\eta (\cdot)$ is equivalent to the existence of a metric $d$ such that $\lim\limits_{t \to +\infty} d\left (p(t, a, \cdot), \eta (\cdot) \right) = 0$.

So, we only need to prove that, for any $a\in \mathbb{R}_+^n \times \mathbb{R}^n$, there is

From Lemma 4.2, $\left \{ p(t, 0, \cdot \mid t\ge 0) \right \}$ is the Cauchy sequence in the space $B(\mathbb{R}_+^n \times \mathbb{R}^n)$ of the metric $d_H$. So, there is a unique $\eta (\cdot)\in B$ such that

By Lemma 4.1 and the triangle inequality, we have

That is, the distribution of $X(t)$ weakly converges to $\eta$.

By the Kolmogorov-Chapman equation, we know that $\eta$ is constant. From Corollary 3.4.3 ^[24], it follows that $\eta$ is strongly mixed. From Theorem 3.2.6 ^[24], we know that $\eta$ 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 $c_1, c_2, \cdots, c_n$ such that

holds, where $C = diag(c_1, c_2, \cdots, c_n).$

Remark 5.1. In Assumption 5.1, the introduction of the constant $c_i, i = 1, 2, \cdots, 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 $c_i$ times. If $c_i \ge 1$, the intraspecific competition intensity and interspecific competition intensity increase by $c_i$ times; if $c_i < 1$, it is weakened by $c_i$ times. Assumption 5.1 means that, under the action of $c_i$, 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))\in \mathbb{R}_+^n \times \mathbb{R}^n$, the solution $(x(t), r(t))$ of system (5.1) has the property that

In particular, if $\max\limits_{1\le i \le n }\bar{r}_i+\max\limits_{1\le i \le n }\left\{a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }\right \}+\limsup\limits_{t\to \infty} \dfrac{1}{t}\int_{0}^{t}\int_{Z}{\ln}(1+ \check{\gamma}(z))v(\mathrm{d}z)\mathrm{d}s < 0$, it implies $\lim\limits_{t\to \infty}\left | x(t) \right | = 0$, and then $x(t)$ is extinct, a.s..

Proof. Define a Lyapunov function

where $c = (c_1, c_2, \cdots, c_n)^T$.

Applying the $\mathrm{It\hat{o} }$ formula, we can get

Using the $\mathrm{It\hat{o} }$ formula for $\ln V(x)$ again, we have

where

where we use the fact $-a_{ii}x_i^{\theta_i}+a_{ii}x_i \le a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }, i = 1, 2, \cdots, n$, and

Substituting the above four inequalities into $\mathrm{d}\ln V(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 $\max\limits_{1\le i \le n }\bar{r}_i+\max\limits_{1\le i \le n }\left\{a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }\right \}+\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\check{\gamma}(z)) v(\mathrm{d}z)\mathrm{d}s < 0$, it implies $\lim\limits_{t\to \infty} \left | x(t) \right | = 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 $\bar{r}_i+\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\gamma_i(z)) v(\mathrm{d}z)\mathrm{d}s > 0, i = 1, 2, \cdots, n$, population $x_i(t), i = 1, 2, \cdots, 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 $\max\limits_{1\le i \le n }\bar{r}_i+\max\limits_{1\le i \le n }\left\{a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }\right \}+\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\check{\gamma}(z)) v(\mathrm{d}z)\mathrm{d}s < 0$, and the parameters of model (1.5) satisfy the conditions of Assumption 5.1, population $x(t) = (x_1(t), \cdots, x_n(t))$ will be extinct. That is, for every $1\le i \le n$, when $\bar{r}_i+\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\gamma_i(z)) v(\mathrm{d}z)\mathrm{d}s < -a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }$, population $x_i(t), i = 1, 2, \cdots, n,$ will be extinct. So, the survival and extinction of the biological population of model (1.5) completely depend on the value of $\bar{r}_i+\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\gamma_i(z)) v(\mathrm{d}z)\mathrm{d}s$.

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 $r_i, i = 1, 2, \cdots, 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 $\bar{r}_i > 0, i = 1, 2, \cdots, n$, populations $x_i(t), i = 1, 2, \cdots, n,$ are weakly persistent; when $\max\limits_{1\le i \le n }\bar{r}_i+\max\limits_{1\le i \le n }\left\{a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }\right \} < 0$, population $x(t) = (x_1(t), \cdots, x_n(t))$ will be extinct. That is, when $\bar{r}_i < -a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }, i = 1, 2, \cdots, n$, populations $x_i(t), i = 1, 2, \cdots, 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 $\bar{r}_i, i = 1, 2, \cdots, n,$ of the population.

When $\bar{r}_i > 0, i = 1, 2, \cdots, n$, the species only disturbed by the OU process are weakly persistent. If the system is affected by jump noise and satisfies $\bar{r}_i+\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\gamma_i(z)) v(\mathrm{d}z)\mathrm{d}s < -a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }, i = 1, 2, \cdots, n$, the species are extinct. When $\bar{r}_i < -a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }, i = 1, 2, \cdots, n$, the species that are only disturbed by the OU process are extinct, but if there are jump noises such that $\bar{r}_i+\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\gamma_i(z)) v(\mathrm{d}z)\mathrm{d}s > 0, i = 1, 2, \cdots, 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 $\gamma_i(z), i = 1, 2, \cdots, n,$ on population survival and extinction. If $\gamma_i(z) < \; 0, i = 1, 2, \cdots, n$, then $\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\gamma_i(z)) v(\mathrm{d}z)\mathrm{d}s < 0, i = 1, 2, \cdots, n$, means that jump noise could accelerate the extinction; if $\gamma_i(z) > 0, i = 1, 2, \cdots, n$, then $\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\gamma_i(z)) v(\mathrm{d}z)\mathrm{d}s > 0, i = 1, 2, \cdots, 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 $x_1(0) = 0.11, x_2(0) = 0.2, r_1(0) = 0.2, r_2(0) = 0.1$, choosing the combination $\mathcal{A}_1$ 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 $x_1, x_2$ 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 $x_1, x_2$ 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 $r_1, r_2$, 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 $r_1(t), r_2(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 $x_1(0) = 0.11, x_2(0) = 0.2, r_1(0) = 0.2, r_2(0) = 0.1$, choosing the combination $\mathcal{A}_2$ 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 $\mathbb{E} (x_1^q), \mathbb{E}(x_2^q)$ are less than $\kappa (q)$, so $\mathbb{E} (x_1^q)\le \kappa (q), \mathbb{E}(x_2^q)\le \kappa(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 $\mathbb{E} (x^2_1)\le \kappa (2), \mathbb{E}(x^2_2)\le \kappa(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 $x_1(0) = 0.11, x_2(0) = 0.2, r_1(0) = 0.2, r_2(0) = 0.2$, choosing the combination $\mathcal{A}_3$ 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 $x_1(t), x_2(t)$, and Figure 3(b), (d) represent the histogram of the solution of $x_1(t), x_2(t)$. According Theorem 4.1, model (6.1) has a stationary distribution $\eta(\cdot)$.

As can be seen from Figure 3(a), (c), the values of population $x_1(t)$ are mostly between 1.5–3, and the values of population $x_2(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 $x_1(t), x_2(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 $\bar{r}_i+\limsup\limits_{t\to \infty}\dfrac{1}{t}\int_{0}^{t}\int_Z \ln (1+\gamma_i(z)) v(\mathrm{d}z)\mathrm{d}s > 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 $x_1(0) = 0.1, x_2(0) = 0.1, r_1(0) = 0.2, r_2(0) = 0.1$, choosing the combination $\mathcal{A}_4$ 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 $\begin{pmatrix} 0 & 0.16\\ 0.12 & 0 \end{pmatrix}$, matrix $S$ is $\begin{pmatrix} 0.4 & 0\\ 0 & 0.5 \end{pmatrix}$, and taking $C = I \in \mathbb{R}^{2\times 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 $\max\limits_{1\le i \le n }\bar{r}_i+\max\limits_{1\le i \le n }\left\{a_{ii}(\theta _i-1)\theta _i^{-\frac{\theta _i}{\theta _i-1} }\right \} < 0$, the populations $x_1(t), x_2(t)$ are extinct when the populations disturbed only by the OU process. The red lines in Figure 4(a), (b) show the populations $x_1(t), x_2(t)$ whose growth rate is disturbed by the OU process. When populations are disturbed only by the OU process, populations $x_1(t), x_2(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 $x_1(t)$ is extinct at $t = 30$ and population $x_2(t)$ is extinct at $t = 45$. In this example, we let $\gamma_1(z) = 0.1, \gamma_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 $x_i(t)$ $(i = 1, 2, \cdots, n)$ have bounded $q$-th moments. The existence of the stationary distribution and extinction of the solution show that when $\bar{r}_i +\limsup\limits_{t \to \infty}\dfrac{1}{t}\int_{0}^{t}\int_{Z}\ln(1+\gamma_i(z))v(\mathrm{d} z)\mathrm{d} s > 0, i = 1, 2, \cdots, n$, model (1.5) has a stationary distribution $\eta (\cdot)$, 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.