% include figures
\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{amssymb}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2013 (2013), No. 162, pp. 1--17.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2013 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2013/162\hfil Asymptotic behavior]
{Asymptotic behavior of stochastic Gilpin-Ayala mutualism model with jumps}

\author[X. Zhang, K. Wang \hfil EJDE-2013/162\hfilneg]
{Xinhong Zhang, Ke Wang}  % in alphabetical order

\address{Xinhong Zhang \newline
Department of Mathematics, Harbin Institute of Technology,
Weihai 264209, China. \newline
College of Science, China University of Petroleum,
Qingdao 266555, China}
\email{zhxinhong@163.com}

\address{Ke Wang \newline
Department of Mathematics, Harbin Institute of Technology,
Weihai 264209, China. \newline
School of Mathematics and Statistics, Northeast Normal University,
 Changchun  130024,  China}
\email{wangke@hitwh.edu.cn, Tel:+86 06315687086}

\thanks{Submitted  May 17, 2013. Published July 19, 2013.}
\subjclass[2000]{34F05, 92D25, 60H10, 60H20}
\keywords{Gilpin-Ayala mutualism model; jumps; moment boundedness;
\hfill\break\indent asymptotic behavior; stochastic permanence}

\begin{abstract}
 This article concerns the study of stochastic Gilpin-Ayala 
 mutualism models with white noise and Poisson jumps. Firstly, 
 an explicit solution for one-dimensional Gilpin-Ayala mutualism model 
 with jumps is obtained  and the asymptotic pathwise behavior is analyzed. 
 Then, sufficient  conditions for the existence of global positive solutions,
 stochastically  ultimate boundedness and stochastic permanence are 
 established for  the n-dimensional model. Asymptotic pathwise behavior 
 of n-dimensional  Gilpin-Ayala mutualism model with jumps is also discussed.
 Finally numerical examples are introduced to illustrate the results developed.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

In nature, mutualism is a usual phenomena. Rhinos and tick birds are
an example of a mutualism relationship. Tick birds eat the ticks on
a rhino while the rhino loses annoying parasites, the relationship is
positive for both the rhino and the tick birds because they both get what
they want. Therefore it is important to study the mutualism models
for multi-species. As is well known now, the most important model among
several  cooperative models is the following Lotka-Volterra mutualism system:
\begin{equation*}
d x_i(t)=x_i(t)\Big[r_i-\sum^{d}_{j=1}a_{ij}x_j(t)\Big]d t, \quad
 1\leq i\leq n,
\end{equation*}
where $x_i(t)$ is the population size of species $i$, $r_i$ is the intrinsic
growth rate of species $i$, $a_{ij}$ ($i\neq j$) represents the effect of
species $j$ upon the growth rate of species $i$,  $a_{ii}$ stands for
the intraspecies interaction, $a_{ii}>0$, $a_{ij}<0$, $i\neq j$.
There is an extensive literature concerned this model, for example,
\cite{m03,m01,m02,gm16,m04,m05,m06,m07}. 
However, in the practical case,
population systems are often subject to various stochastic small perturbation.
The growth rates, interaction coefficients and so on may be influenced by
environmental noise. In recent years, stochastic differential equations
have received much attention, many results have been derived to reveal
how environmental noise affects the population systems.
In particular,  Mao, Marion and Renshaw \cite{gm01}  revealed that
the environmental noise can suppress a potential population explosion.
Suppose the growth rates are perturbed by white noise
$r_i\to r_i+\sum_{j=1}^{m}\sigma_{ij}\dot{B}_j(t)$, here
$\dot{B}_j(t)$ is a white noise, i.e., $(B_1(t),\dots,B_m(t))$ is an
$m$-dimensional Brownian motion defined on a complete probability space
 $(\Omega, \mathbb{F}, \mathbb{P})$, $\sigma_{ij}^2$ stands for the
intensity of the noise. Then the stochastic Lotka-Volterra mutualism
system becomes
\begin{equation*}
d x_i(t)=x_i(t)\Big[(r_i-\sum^{d}_{j=1}a_{ij}x_j(t))dt
+\sum_{j=1}^{m}\sigma_{ij}d B_j(t)\Big], \quad 1\leq i\leq n.
\end{equation*}
Various forms of cooperative Lotka-Volterra system have been extensively
studied and we here mention Hung \cite{m08}, Cheng\cite{m09}, Liu and
Wang \cite{m10,m11} and the references cited therein. The key method
used in our paper is motivated by them.

Unfortunately, in the Lotka-Volterra model, the rate of change of the
size of each species is linear function of sizes of the interacting
species \cite{gm05,gm06}. However in complex ecosystem this is almost impossible.
Therefore, to meet the practical situations, in 1973,
Gilpin and Ayala \cite{gm21}  provided a modification for Lotka-Volterra model,
 called Gilpin-Ayala model. For various forms about the Gilpin-Ayala
system readers can see \cite{ga06,ga07,gm17,gm18} and references therein
for details. But here are few works about stochastic Gilpin-Ayala mutualism model.

On the other hand, the population systems may suffer sudden environmental
perturbations, that is, some jump type stochastic perturbations;
 e.g., earthquakes, hurricanes, epidemics and so on \cite{gs07,gs08}.
These phenomena can not be described by stochastic integrals driven only
 by Brownian motion. So it is feasible to introduce a jump process into
the underlying population system. SDEs with jumps have received considerable
 attention in the past few years. We here mention Applebaum \cite{gs09},
Situ \cite{gs12}, Bao et al \cite{gs07,gs08}. Particularly,
the books by Applebaum \cite{gs09} and Situ \cite{gs12} are good
references in this area. To the best of the author¡¯s knowledge,
to this day, $n$-dimensional Gilpin-Ayala  mutualism model with jumps
 has not been studied. Motivated by these, the following $n$-dimensional
Gilpin-Ayala  mutualism model with jumps is considered in this article:
\begin{equation}\label{m0}
\begin{aligned}
d x_i(t)&=x_i(t^{-})\Big\{\Big[r_i-a_{ii}x_i^{\theta_i}(t^{-})
-\sum^{d}_{j\neq i}a_{ij}x_j(t^{-})\Big]dt
+\sum^{m}_{j=1}\sigma_{ij}d B_j(t)\\
&\quad +\int_ \mathbb{Y}\gamma_i(u)N(d t,d u)\Big\},
\end{aligned}
\end{equation}
for $1\leq i\leq n$,
where $x(t^{-})$ is the left limit of $x(t)$, $\theta_i\geq 1$
is the parameter to modify the classical Lotka-Volterra model,
$\gamma_i(u)>-1$ is a bounded function, $i=1,\dots,n$, $N$ is
a Poisson counting measure with characteristic measure $\nu$ on a
measurable subset $\mathbb{Y}$ of $(0,\infty)$ with $\nu(\mathbb{Y})<\infty$,
and $\widetilde{N}(d t,d u):=N(d t,d u)-\nu(d u)d t$.
We assume Brownian motion and $N$ are independent.

Throughout this article
 $\mathbb{R}_+^{n}:=\{x\in \mathbb{R}^{n}:x_i>0\text{ for  } i=1,\dots,n\}$,
$\bar{A}=(a_{ij})$, $\bar{A}^T$ denotes the transpose of
 $\bar{A}$. If $x\in \mathbb{R}^n$, its norm is denoted by
$|x|=(\sum^{n}_{i=1}x_i^2)^{\frac{1}{2}}$. If $Q$ is a matrix,
 $|Q|=\sqrt{\operatorname{trace}(Q^TQ)}$ represents its trace norm.
If $Q=(q_{ij})_{n\times n}$ is a symmetric matrix, then
 $\lambda_{\max}^+(Q)=\sup_{x\in \mathbb{R}_+^n, |x|=1}x^TQx$.
$K$ is a positive constant and may be different at different places.
We impose the following assumptions:
\begin{itemize}
\item[(A1)] There  is a positive constant $\kappa$ such that
$$
\int_{\mathbb{Y}}(|\ln(1+\gamma (y))|\vee |\ln(1+\gamma (y))|^2)\nu(d y)<\kappa.
$$

\item[(A2)] There are positive constants $p_1,\dots,p_n$ such that
$$
-2\tau := \lambda_{\max}^+(-\bar{P}\bar{A}-\bar{A}^T\bar{P})<0,
$$
where $\bar{P}=\operatorname{diag}(p_1,\dots,p_n)$.

\end{itemize}
From Assumption (A2), it is easy to see that
$\tau\leq a_{ii}p_i$ for $i=1,\dots,n$.

The aim of our work is to study the properties of $n$-dimensional
stochastic Gilpin-Ayala mutualism model with jumps. The significance
of this paper is mainly: (1) Gilpin-Ayala system is more suitable for
the real situations than Lotka-Volterra system, but more complicated;
(2) The white noise and Poisson jumps are taken into account.
The remaining part of this paper is organized as follows.
In section 2, we provide an explicit solution for one-dimensional
Gilpin-Ayala model with jumps and study its asymptotic pathwise behavior.
In section 3,  we show that \eqref{m0} will have a unique global
positive solution under certain conditions. Section 4 and 5 deal
with the asymptotic moment properties and asymptotic pathwise behavior
of the solution, respectively. In section 6, we show that the system
is stochastically permanent if the white noise and Poisson jumps satisfy
our conditions. Finally we introduce some simulation figures to
illustrate our main results.

\section{One-dimensional Gilpin-Ayala model}\label{mean}

As the single population is the basic unit of the whole ecological system,
the establishment and theoretical analysis of the single population model
can help us to understand the overall structure of the complex model.
So we firstly analyze one-dimensional Gilpin-Ayala model with jumps.

\begin{lemma}[\cite{gs07}] \label{lem1}
 Consider the following system of  equations with jumps:
$$
d Y_i(t)=Y_i(t^-)\Big[(a_i-b_{ii}Y_i(t^-))dt
+\sigma_i d B(t)+\int_{\mathbb{Y}}c_i(u)\widetilde{N}(dt,d u)\Big],
$$
where $a_i>0$, $b_{ii}>0$, $c_i(u)>-1$, $B(t)$ is a one-dimensional
Brownian motion. Then for any initial value $Y_i(0)\in \mathbb{R}_+^n$,
this equation admits a unique positive solution $Y_i(t)$, $t\geq 0$,
which is global and admits the explicit formula
$$
Y_i(t)=\frac{\varphi_i(t)}{\frac{1}{Y_i(0)}+\int_{0}^{t}b_{ii}\varphi_i(s)d s},
$$
where
\begin{align*}
\varphi_i(t)
&:=\exp\Big(\int_{0}^{t}\Big[a_i-\frac{1}{2}\sigma_i^2
+\int_\mathbb{Y}(\ln(1+c_i(u))-c_i(u))\nu(d u)\Big]ds
+\int_{0}^{t}\sigma_id B(s)\\
&\quad +\int_{0}^{t}\int_\mathbb{Y}\ln(1+c_i(u))\widetilde{N}(d t,d u)\Big).
\end{align*}
\end{lemma}

\begin{remark} \label{rmk1}\rm
In general, the intrinsic growth rate $a_i$ is positive, but the above
 explicit solution holds for $a_i\leq 0$.
\end{remark}

The one-dimensional Gilpin-Ayala model with jumps is
\begin{equation}\label{m41}
\begin{gathered}
\begin{aligned}
dy_i(t)&=y_i(t^-)(r_i-a_{ii}y_i^{\theta_i}(t^-))d t
+y_i\sum_{j=1}^{m}\sigma_{ij} d B_j(t)\\
&\quad +\int_{\mathbb{Y}}y_i(t^-)\gamma_i(u)N(d t,d u),
\quad  \theta_i>1,
\end{aligned}\\
y_i(0)=x_i(0).
\end{gathered}
\end{equation}
Set $z_i=y_i^{\theta_i}$, by It\^o's formula, we have
\begin{equation}\label{m42}
\begin{gathered}
\begin{aligned}
d z_i(t)
&=z_i(t^-)\Big[\theta_ir_i+\frac{\theta_i(\theta_i-1)}{2}
 \sum_{j=1}^{m}\sigma_{ij}^2
 +\int_{\mathbb{Y}}((1+\gamma_i(u))^{\theta_i}-1) \nu(d u)\\
&\quad-\theta_ia_{ii}z_i(t^-)\Big]d t
 +z_i(t^-)\theta_i\sum_{j=1}^{m}\sigma_{ij} d B_j(t)\\
&\quad +\int_{\mathbb{Y}}z_i(t^-)((1+\gamma_i(u))^{\theta_i}-1)
\widetilde{N}(d t,d u), \quad \theta_i>1,
\end{aligned}\\
z_i(0)=x_i^{\theta_i}(0).
\end{gathered}
\end{equation}
From Lemma \ref{lem1} and Remark \ref{rmk1}, it follows that
 \eqref{m42} has an explicit solution
$$
z_i(t)=\frac{\Phi_i(t)}{\frac{1}{x_i^{\theta_i}(0)}
+\int_{0}^{t}a_{ii}{\theta_i}\Phi_i(s)d s},
$$
where
\begin{align*}
\Phi_i(t)&:=\exp\Big\{\int_{0}^{t}\theta_i
\Big[r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2
 +\int_\mathbb{Y}\ln(1+\gamma_i(u))\nu(d u)\Big]d s\\
&\quad +\int_{0}^{t}\theta_i\sum_{j=1}^{m}\sigma_{ij}d B_j(s)
+\int_{0}^{t}\int_\mathbb{Y}\theta_i\ln(1+\gamma_i(u))\widetilde{N}(d t,d u)
\Big\}
\end{align*}
Combining \cite[Lemma 4.4 and Theorem 4.4]{gs07}, we can deduce the
following lemma.

\begin{lemma} \label{lem2}
 Assume {\rm (A1)} holds. If
$c_i:=r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2
+\int_{\mathbb{Y}}\ln(1+\gamma_i(u))\nu(d u)\geq 0$, then for
$i=1,\dots,n$, we have
\begin{equation}\label{m43}
\lim_{t\to\infty}\frac{\ln z_i(t)}{t}=0\ \ a.s.
\end{equation}
\end{lemma}

\begin{remark} \label{rmk2} \rm
Based on the above analysis,  \eqref{m41} has a unique positive solution
$y_i(t)$ for any value $y_i(0)=x_i(0)>0$ which is global and represented by
$$
y_i(t)=\Big(\frac{\Phi_i(t)}{\frac{1}{x_i^{\theta_i}(0)}
+\int_{0}^{t}a_{ii}{\theta_i}\Phi_i(s)d s}\Big)^{1/\theta_i},
$$
where $\Phi_i(t)$ is defined as above.
Under the conditions of Lemma \ref{lem2}, we obtain
$\lim_{t\to\infty}\frac{\ln y_i(t)}{t}=0$ a.s.
\end{remark}

\begin{theorem} \label{thm1}
Suppose that $y_i(t)$ is a positive solution of \eqref{m41}.
If $c_i\geq0$, then
$$
\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}y_i^{\theta_i}(s)d s
=\frac{r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2
+\int_{\mathbb{Y}}\ln(1+\gamma_i(u))\nu(d u)}{a_{ii}} \quad\text{a.s.}
$$
\end{theorem}

\begin{proof}
Applying It\^o's formula to $\ln y_i(t)$ results in
\begin{align*}
d \ln y_i(t)
&=(r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2+\int_{\mathbb{Y}}
 \ln(1+\gamma_i(u))\nu(d u)-a_{ii}y_i^{\theta_i})d t\\
&\quad +\sum_{j=1}^{m}\sigma_{ij} d B_j(t)+\int_{\mathbb{Y}}
 \ln(1+\gamma_i(u))\widetilde{N}(d t,d u).
\end{align*}
Integrating from 0 to $t$ yields
\begin{align*}
\ln y_i(t)-\ln y_i(0)
&=\Big(r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2
+\int_{\mathbb{Y}}\ln(1+\gamma_i(u))\nu(d u)\Big)t\\
&\quad -a_{ii} \int_{0}^{t}y_i^{\theta_i}(s)d s+M_1(t)+M_2(t),
\end{align*}
where $M_1(t)=\int_{0}^{t}\sum_{j=1}^{m}\sigma_{ij} d B_j(s)$,
 $M_2(t)=\int_{0}^{t}\int_{\mathbb{Y}}\ln(1+\gamma_i(u))\widetilde{N}(d s,d u)$
 are real valued local martingales vanishing at $t=0$.
Hence
\begin{equation}\label{m44}
\begin{aligned}
\frac{\ln y_i(t)}{t}-\frac{\ln y_i(0)}{t}
&=\Big(r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2
+\int_{\mathbb{Y}}\ln(1+\gamma_i(u))\nu(d u)\Big)\\
&\quad -\frac{a_{ii}}{t}\int_{0}^{t}y_i^{\theta_i}(s)d s+\frac{M_1(t)}{t}
+\frac{M_2(t)}{t}.
\end{aligned}
\end{equation}
Then by \cite[Proposition 2.4]{gs13},
\begin{gather*}
\langle M_1\rangle (t)=\int_{0}^{t}\sum_{j=1}^{m}\sigma_{ij}^2 d s
=\sum_{j=1}^{m}\sigma_{ij}^2 t,
\\
\langle M_2\rangle (t)
=\int_{0}^{t}\int_{\mathbb{Y}}[\ln(1+\gamma_i(u))]^2 \nu(d u)d s
= t\int_{\mathbb{Y}}[\ln(1+\gamma_i(u))]^2 \nu(d u),
\end{gather*}
where $\langle M\rangle (t):=\langle M,M\rangle$ is Meyer's angle bracket
 process.
 We have
\begin{equation*}
\int_{0}^{t}\frac{1}{(1+s)^2}d s=\frac{t}{t+1}<\infty,
\end{equation*}
by the strong law of large numbers for local martingales \cite{gs14},
we then obtain
\begin{equation*}
\lim_{t\to\infty}\frac{M_1(t)}{t}=0 \text{ a.s.}, \quad
\lim_{t\to\infty}\frac{M_2(t)}{t}=0 \text{ a.s.}
\end{equation*}
Taking limits on both sides of \eqref{m44} and combining
Remark \ref{rmk2} lead to
$$
\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}y_i^{\theta_i}(s)d s
=\frac{r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2
+\int_{\mathbb{Y}}\ln(1+\gamma_i(u))\nu(d u)}{a_{ii}}\quad\text{a.s. }
$$
This completes the proof.
\end{proof}


\section{Global positive solutions of \eqref{m0}}\label{mean2}

As $x_i(t)$ in \eqref{m0} denotes the size of species $i$,
it should be nonnegative.  
To guarantee that the stochastic differential
equations (SDEs) have a unique global solution for any given 
initial value, the coefficients of the equation are generally 
required to satisfy both the linear growth
condition and the local Lipschitz condition (see e.g.\cite{gm12,gm19}).
But we can find that the coefficients of \eqref{m0} are locally Lipschitz
continuous, and they do not satisfy the linear growth condition.
So the solution of \eqref{m0} may explode in a finite time.
The following theorem gives sufficient condition for global positive solutions.

\begin{theorem} \label{thm2}
Let {\rm (A1), (A2)} hold and $\theta_i\geq 1$, $i=1,\dots,n$.
 Then for any initial value $x_0\in \mathbb{R}_+^n$, Equation \eqref{m0}
has a unique global solution $x(t)\in \mathbb{R}_+^n$ for all
$t\geq 0$ almost surely.
\end{theorem}

\begin{proof}
Our proof is motivated by Mao, Marion and Renshaw \cite{gm01}.
Clearly, the coefficients of \eqref{m0} are locally Lipschitz continuous,
so for any initial value $x_0\in \mathbb{R}_+^n$ Equation \eqref{m0} has
 a unique maximal local solution $x(t)$ on $t\in [0,\tau_e)$, where $\tau_e$
is the explosion time. If we show that $\tau_e=\infty$ a.s., then the
solution is global. Now let $k_0$ be big enough for every component of
$x_0$ lying within the interval $[1/k_0,k_0]$. For any integer $k\geq k_0$,
define the stopping time
$$
\tau_k=\inf\{t\in [0,\tau_e):x_i(t)\notin(1/k,k)\text{ for some }
i=1,\dots,n\},
$$
where we set $\inf\emptyset=\infty$. Obviously,
 $\tau _k$ is increasing as $k\to\infty$. Set
$\tau_\infty=\lim_{k\to\infty}\tau_k$, whence $\tau_\infty\leq \tau_e$ a.s.
Now all we need to show is $\tau_\infty=\infty$ a.s.
If this assertion is false, then there is a pair of constants
$T>0$ and $\epsilon\in (0,1)$ such that
$\mathbb{P}\{\tau_\infty\leq T\}>\epsilon$.
Therefore, there is an integer $k_1\geq k_0$ such that
\begin{equation}\label{m21}
\mathbb{P}\{\tau_k\leq T\}\geq\epsilon \quad\text{for \ all } k\geq k_1.
 \end{equation}
Define a $C^2$-function $V:\mathbb{R}_+^n\to \mathbb{R}_+$ by
 $$
V(x)=\sum_{i=1}^{n}p_i(x_i-1-\ln x_i).
$$
If $x(t)\in \mathbb{R}_+^n$, It\^{o}'s  formula shows that
\begin{align*}
&d V(x(t))\\
&=\sum_{i=1}^{n}\Big\{p_i\Big[(x_i-1)(r_i-a_{ii}x_i^{\theta_i}
-\sum_{j\neq i}^{n}a_{ij}x_j)+0.5\sum_{j=1}^{m}\sigma_{ij}^2\Big]\Big\}d t\\
&\quad +\sum_{i=1}^{n}p_i(x_i-1)\sum_{j=1}^{m}\sigma_{ij}d B_j(t)
+\sum_{i=1}^{n}\int_\mathbb{Y}p_i[x_i\gamma_i(u)-\ln(1+\gamma_i(u))]N(d t,d u)\\
&=\sum_{i=1}^{n}\Big\{p_i\Big[(x_i-1)(r_i-a_{ii}x_i^{\theta_i}
 -\sum_{j\neq i}^{n}a_{ij}x_j)+0.5\sum_{j=1}^{m}\sigma_{ij}^2\\
&\quad +\int_\mathbb{Y}(x_i\gamma_i(u)-\ln(1+\gamma_i(u)))\nu(d u)\Big]\Big\}d t
 +\sum_{i=1}^{n}p_i(x_i-1)\sum_{j=1}^{m}\sigma_{ij}d B_j(t)\\
&\quad +\sum_{i=1}^{n}\int_\mathbb{Y}p_i[x_i\gamma_i(u)-\ln(1+\gamma_i(u))]
\widetilde{N}(d t,d u),
\end{align*}
where we drop $t^{-}$ from $x(t^{-})$.
Using (A1) and (A2), we get that there exists a positive constant $K$ 
such that
\begin{align*}
&\sum_{i=1}^{n}p_i\Big[(x_i-1)(r_i-a_{ii}x_i^{\theta_i}
 -\sum_{j\neq i}^{n}a_{ij}x_j)+0.5\sum_{j=1}^{m}\sigma_{ij}^2\\
& +\int_\mathbb{Y}(x_i\gamma_i(u)-\ln(1+\gamma_i(u)))\nu(d u)\Big]
\\
&\leq \sum_{i=1}^{n}p_i\Big[(r_i+\int_\mathbb{Y}\gamma_i(u)\nu(d u))x_i
 +a_{ii}x_i^2-a_{ii}x_i^{\theta+1}+a_{ii}x_i^{\theta}-r_i
 +0.5\sum_{j=1}^{m}\sigma_{ij}^2\\
&\quad -\int_\mathbb{Y}\ln(1+\gamma_i(u))\nu(d u)\Big]
+0.5x^T(-\bar{P}\bar{A}-\bar{A}^T\bar{P})x
\\
&\leq \sum_{i=1}^{n}\Big\{-a_{ii}p_ix_i^{\theta_i+1}
 -(\tau-a_{ii}p_i) x_i^2+p_ia_{ii}x_i^{\theta_i}
 +p_i\Big(r_i+\int_\mathbb{Y}\gamma_i(u)\nu(d u)\Big)x_i \\
&\quad +p_i\Big(-r_i+0.5\sum_{j=1}^{m}\sigma_{ij}^2
-\int_\mathbb{Y}\ln(1+\gamma_i(u))\nu(d u)\Big)\Big\}
\leq  K.
\end{align*}
Therefore,
\begin{align*}
\int_{0}^{\tau_k\wedge T}d V(x(t))
&\leq \int_{0}^{\tau_k\wedge T}K d t
 +\int_{0}^{\tau_k\wedge T}\sum_{i=1}^{n}p_i(x_i(t)-1)
 \sum_{j=1}^{m}\sigma_{ij}d B_j(t)\\
&\quad +\int_{0}^{\tau_k\wedge T}\int_\mathbb{Y}\sum_{i=1}^{n}
\int_\mathbb{Y}p_i[x_i\gamma_i(u)-\ln(1+\gamma_i(u))]\widetilde{N}(d t,d u).
\end{align*}
Taking expectations on both sides results in
\begin{equation}\label{m22}
\mathbb{E}V(x(\tau_k\wedge T))
\leq V(x_0)+K\mathbb{E}(\tau_k\wedge T)\leq V(x_0)+KT.
 \end{equation}
Set $\Omega_k=\{\tau_k\leq T\}$ for $k\geq k_1$, from \eqref{m21}
 we have $\mathbb{P}(\Omega_k)\geq \epsilon$.
 Note that for every $\omega\in \Omega_k$, there is some $i$ such that
$x_i(\tau_k,\omega)$ equals either $k$ or $1/k$, hence
$$
V(x(\tau_k,\omega))\geq p_i[k-1-\ln(k)]\wedge p_i[1/k-1-\ln(1/k)].
$$
Using \eqref{m22}, yields
 $$
V(x_0)+KT\geq \mathbb{E}\left(I_{\Omega_k}V(x(\tau_k,\omega))\right)
\geq \epsilon \left( p_i[k-1-\ln(k)]\wedge p_i[1/k-1-\ln(1/k)]\right),
$$
 where $I_{\Omega_k}$ is the indicator function of $\Omega_k$.
When $k\to\infty$ we obtain
 $$
\infty > V(x_0)+KT=\infty,
$$
it results in $\tau_\infty=\infty$ a.s. The proof is complete.
\end{proof}

\section{Ultimate boundedness}

In the previous section, we saw that \eqref{m0} has a unique global
solution $x(t)\in \mathbb{R}_+^n$ for any $t\geq 0$ almost surely.
Based on this fundamental theorem, we discuss the ultimate boundedness
and asymptotic boundedness in any $p$th moment of the solutions.

\begin{theorem} \label{thm3}
Under the assumptions of Theorem \ref{thm2},
for any initial value $x_0\in \mathbb{R}_+^n$, the solution of \eqref{m0}
satisfies
\begin{equation*}
\limsup_{t\to\infty}\mathbb{E}|x(t)|\leq K.
\end{equation*}
\end{theorem}

\begin{proof}
Define the Lyapunov function
$$
V(x):=\sum_{i=1}^{n}p_ix_i,\ \ x\in \mathbb{R}_+^n.
$$
Applying It\^o's formula, we obtain
\begin{equation}\label{m31}
\begin{split}
d V(x(t))
&=\sum_{i=1}^{n}p_ix_i(r_i-a_{ii}x_i^{\theta_i}
 -\sum_{j\neq i}^{n}a_{ij}x_j)d t
 +\sum_{i=1}^{n}p_ix_i\sum_{j=1}^{m}\sigma_{ij}d B_j(t)\\
&\quad +\int_\mathbb{Y}\sum_{i=1}^{n}p_ix_i\gamma_i(u)N(d t,d u)\\
&=L V(x)d t+\sum_{i=1}^{n}p_ix_i\sum_{j=1}^{m}\sigma_{ij}d B_j(t)
+\int_\mathbb{Y}\sum_{i=1}^{n}p_ix_i\gamma_i(u)\widetilde{N}(d t,d u),
\end{split}
\end{equation}
where we write $x(t^-)=x$, and
\begin{align*}
L V(x)
&=\sum_{i=1}^{n}p_i\Big(r_i-a_{ii}x_i^{\theta_i}-\sum_{j\neq i}^{n}a_{ij}x_j
 +\int_\mathbb{Y}\gamma_i(u)\nu(d u)\Big)x_i\\
&=\sum_{i=1}^{n}p_i \Big(r_i-a_{ii}x_i^{\theta_i}+a_{ii}x_i
 +\int_\mathbb{Y}\gamma_i(u)\nu(d u)\Big)x_i+x^T(-\bar{P}\bar{A})x\\
&=\sum_{i=1}^{n}p_i\Big(r_i-a_{ii}x_i^{\theta_i}+a_{ii}x_i
 +\int_\mathbb{Y}\gamma_i(u)\nu(d u)\Big)x_i+0.5x^T(-\bar{P}\bar{A}
 -\bar{A}^T\bar{P})x\\
&\leq \sum_{i=1}^{n}\Big[-a_{ii}p_ix_i^{\theta_i+1}-(\tau-a_{ii}p_i) x_i^2
 +p_i\Big(r_i+\int_\mathbb{Y}\gamma_i(u)\nu(d u)\Big)x_i\Big].
\end{align*}
For arbitrary $\alpha>0$, making use of the conditions of this Theorem,
applying It\^o's formula once again yields
\begin{align*}
&d (e^{\alpha t}V(x(t)))\\
&= \alpha e^{\alpha t}V(x(t))d t+e^{\alpha t}d V(x(t))\\
&\leq e^{\alpha t}\sum_{i=1}^{n}[-a_{ii}p_ix_i^{\theta_i+1}
 -(\tau-a_{ii}p_i) x_i^2+p_i(\alpha +r_i
 +\int_\mathbb{Y}\gamma_i(u)\nu(d u))x_i]d t \\
&\quad +e^{\alpha t}\sum_{i=1}^{n}p_ix_i\sum_{j=1}^{m}\sigma_{ij}d B_j(t)
 + e^{\alpha t}\int_\mathbb{Y}\sum_{i=1}^{n}p_ix_i\gamma_i(u)
 \widetilde{N}(d t,d u)\\
&\leq  K_0 e^{\alpha t}d t+e^{\alpha t}
 \sum_{i=1}^{n}p_ix_i\sum_{j=1}^{m}\sigma_{ij}d B_j(t)
 +e^{\alpha t}\int_\mathbb{Y}\sum_{i=1}^{n}p_ix_i
 \gamma_i(u)\widetilde{N}(d t,d u),
\end{align*}
where $K_0$ is a positive constant.
Therefore,
$$
\mathbb{E}(e^{\alpha t}V(x(t)))\leq V(x_0)+\frac{K_0}{\alpha}(e^{\alpha t}-1);
$$
that is to say
\begin{equation}\label{m32}
\limsup_{t\to\infty}\mathbb{E}V(x(t))\leq \frac{K_0}{\alpha}.
\end{equation}
Noting that
$|x(t)|\leq \sum_{i=1}^{n}x_i(t)\leq \frac{V(x(t))}{\min_{1\leq i\leq n}p_i}$,
 we obtain
$$
\limsup_{t\to\infty}\mathbb{E}|x(t)|\leq \frac{K_0}{\alpha\min_{1\leq i\leq n}p_i}
=:K.
$$
This completes the proof.
\end{proof}

\begin{definition}[\cite{gm03}]  \label{def1} \rm
The solution of\eqref{m0} is said to be stochastically ultimately bounded
if for any $\epsilon\in(0,1)$, there is a constant $H=H(\epsilon)$
such that for any $x_0\in \mathbb{R}_+^n$,
$$
\limsup_{t\to\infty}\mathbb{P}\left\{|x(t)|>H\right\}<\epsilon.
$$
\end{definition}

As an application of Theorem \ref{thm3}, together with the Chebyshev inequality,
we have the following corollary.

\begin{corollary} \label{coro1}
Under the conditions of Theorem \ref{thm3},
the solution of \eqref{m0} is stochastically ultimately bounded.
\end{corollary}

Furthermore, we can get the following property.

\begin{theorem} \label{thm4}
Assume {\rm (A1), (A2)} hold and $\theta_i>1$, $i=1,\dots,n$.
Then for $p>0$, there exists a positive constant $K=K(p)$,
for any initial value $x_0\in \mathbb{R}_+^n$, the solution of \eqref{m0}
has the property
\begin{equation*}
\limsup_{t\to\infty}\mathbb{E}x_i^p(t)\leq K(p) ,\quad t\geq 0,\; i=1,\dots,n.
\end{equation*}
\end{theorem}

\begin{proof}
Define the Lyapunov function
$$
V(x,t):=\sum_{i=1}^{n}e^tx_i^p,\quad x\in \mathbb{R}_+^n.
$$
Applying It\^o's formula, we obtain
\begin{align*}
d V(x(t),t)
&=L V(x(t))d t+e^tp\sum_{i=1}^{n}x_i^p\sum_{j=1}^{m}\sigma_{ij}d B_j(t)\\
&\quad +e^t\sum_{i=1}^{n}x_i^p\int_\mathbb{Y}[(1+\gamma_i(u))^p-1]
\widetilde{N}(d t,d u),
\end{align*}
where we write $x(t^-)=x$, and
\begin{align*}
L V(x)&=e^t\sum_{i=1}^{n}p\Big\{\Big[\frac{1}{p}+r_i
 +\frac{p-1}{2}\sum_{j=1}^{m}\sigma_{ij}^2+\frac{1}{p}\int_\mathbb{Y}
 [(1+\gamma_i(u))^p-1]\nu (d u)\Big]x_i^p\\
&\quad -a_{ii}x_i^{p+\theta_i}-\sum_{j\neq i}^{n}a_{ij}x_i^px_j\Big\}
\\
&\leq e^t\sum_{i=1}^{n}p\Big\{\Big[\frac{1}{p}+r_i+\frac{p-1}{2}
 \sum_{j=1}^{m}\sigma_{ij}^2+\frac{1}{p}\int_\mathbb{Y}[(1+\gamma_i(u))^p-1]
 \nu (d u)\Big]x_i^p
\\
&\quad -a_{ii}x_i^{p+\theta_i}-\sum_{j\neq i}^{n}a_{ij}
 \Big[\frac{px_i^{p+1}}{p+1}+\frac{x_j^{p+1}}{p+1}\Big]\Big\}
\\
&\leq e^t\sum_{i=1}^{n}p\Big\{\Big[\frac{1}{p}+r_i
+\frac{p-1}{2}\sum_{j=1}^{m}\sigma_{ij}^2+\frac{1}{p}
 \int_\mathbb{Y}[(1+\gamma_i(u))^p-1]\nu (d u)\Big]x_i^p
\\
&\quad -a_{ii}x_i^{p+\theta_i}
 -\Big[\sum_{j\neq i}^{n}\Big(\frac{p}{p+1}(a_{ij})+\frac{1}{p+1}(a_{ji})\Big)
 \Big]x_i^{p+1}\Big\}
\\
&\leq e^t\sum_{i=1}^{n}K_i(p),
\end{align*}
where $K_i(p)$ is a positive constant.
Hence
\begin{equation*}
e^t\mathbb{E}[\sum_{i=1}^{n}x_i^p(t)]\leq \sum_{i=1}^{n}x_i^p(0)
+\mathbb{E}\int_{0}^{t}e^s\sum_{i=1}^{n}K_i(p)d s
=\sum_{i=1}^{n}x_i^p(0)+\sum_{i=1}^{n}K_i(p)(e^t-1).
\end{equation*}
It is not difficult to derive that
\begin{equation*}
\limsup_{t\to\infty}\mathbb{E}[\sum_{i=1}^{n}x_i^p(t)]
\leq \sum_{i=1}^{n}K_i(p)=:K(p).
\end{equation*}
The required assertion follows immediately.
\end{proof}

\section{Pathwise estimation}

In this section we consider the asymptotic pathwise estimation of the
solution to \eqref{m0}.

\begin{theorem} \label{thm5}
For $\theta_i>1$, $i=1,\dots,n$, under Assumptions {\rm (A1), (A2)},
for any initial value $x_0\in \mathbb{R}_+^n$, the solution of \eqref{m0}
has the property
\begin{equation*}
\limsup_{t\to\infty}\frac{\ln x_i(t)}{\ln t}\leq 1 \quad\text{a.s., }
i=1,\dots,n.
\end{equation*}
\end{theorem}

\begin{proof}
Here we adopt the same notation as in the proof of Theorem \ref{thm3}.
From \eqref{m31}, by simple manipulation, one has
\begin{align}
&\mathbb{E}\Big(\sup_{t\leq u\leq t+1}V(x(u))\Big) \nonumber \\
&\leq  \mathbb{E}(V(x(t)))+\mathbb{E}
\Big(\sup_{t\leq u\leq t+1}\int_{t}^{u}\sum_{i=1}^{n}
\Big[-a_{ii}p_ix_i^{\theta_i+1}(s)-(\tau-a_{ii}p_i) x_i^2(s)
\nonumber  \\
&\quad +p_i\Big(r_i+\int_\mathbb{Y}\gamma_i(u)\nu(d u)\Big)x_i(s)\Big]d s\Big)
+\mathbb{E}\Big(\sup_{t\leq u\leq t+1}\int_{t}^{u}\sum_{i=1}^{n}p_ix_i(s)
\sum_{j=1}^{m}\sigma_{ij}d B_j(s)\Big) \nonumber \\
&\quad +\mathbb{E}\Big(\sup_{t\leq u\leq t+1}
 \int_{t}^{u}\int_\mathbb{Y}\sum_{i=1}^{n}p_ix_i(s)\gamma_i(u)\widetilde{N}
(d s,d u)\Big)
\nonumber  \\
&\leq  \mathbb{E}(V(x(t)))
+\mathbb{E}\Big(\sup_{t\leq u\leq t+1}
 \int_{t}^{u}\sum_{i=1}^{n}\Big[a_{ii}p_i x_i^2(s) \nonumber \\
&\quad +p_i\Big(r_i+\int_\mathbb{Y}|\gamma_i(u)|\nu(d u)\Big)x_i(s)\Big]d s\Big)
\nonumber \\
&\quad +\mathbb{E}\Big(\sup_{t\leq u\leq t+1}\int_{t}^{u}
 \sum_{i=1}^{n}p_ix_i(s)\sum_{j=1}^{m}\sigma_{ij}d B_j(s)\Big)\nonumber \\
&\quad +\mathbb{E}\Big(\sup_{t\leq u\leq t+1}|\int_{t}^{u}\int_\mathbb{Y}
 \sum_{i=1}^{n}p_ix_i(s)\gamma_i(u)\widetilde{N}(d s,d u)|\Big)
\nonumber \\
&\leq \mathbb{E}(V(x(t)))+q_1\int_{t}^{t+1}\mathbb{E}|x(s)|d s
 +q_2\int_{t}^{t+1}\mathbb{E}|x(s)|^2d s\nonumber \\
&\quad +\mathbb{E}\Big(\sup_{t\leq u\leq t+1}\int_{t}^{u}
 \sum_{i=1}^{n}p_ix_i(s)\sum_{j=1}^{m}\sigma_{ij}d B_j(s)\Big)
\nonumber \\
&\quad +\mathbb{E}\Big(\sup_{t\leq u\leq t+1}|\int_{t}^{u}\int_\mathbb{Y}
 \sum_{i=1}^{n}p_ix_i(s)\gamma_i(u)\widetilde{N}(d s,d u)|\Big), \label{m33}
\end{align}
where $q_1=\sqrt{n}\max_{1\leq i\leq n}\{p_i(r_i+\int_\mathbb{Y}|
\gamma_i(u)|\nu(d u))\}$, $q_2=\max_{1\leq i\leq n}\{a_{ii}p_i\}$.
By the Burkholder-Davis-Gundy inequality for local martingale
(see, e.g., \cite{gm12,gs09}) and the H\"older inequality, we obtain that
\begin{align*}
\mathbb{E}\Big(\sup_{t\leq u\leq t+1}\int_{t}^{u}
 \sum_{i=1}^{n}p_ix_i(s)\sum_{j=1}^{m}\sigma_{ij}d B_j(s)\Big)
&\leq 3\mathbb{E}\Big(\int_{t}^{t+1}|x^T\bar{P}\sigma|^2d s\Big)^{1/2}\\
&\leq 3|\bar{P}\sigma|\Big(\mathbb{E}\int_{t}^{t+1}|x(s)|^2d s\Big)^{1/2},
\end{align*}
where $\sigma=(\sigma_{ij})$,
and
\begin{align*}
&\mathbb{E}\Big(\sup_{t\leq u\leq t+1}|\int_{t}^{u}
 \int_\mathbb{Y}\sum_{i=1}^{n}p_ix_i(s)\gamma_i(u)\widetilde{N}(d s,d u)|\Big)\\
&\leq \sum_{i=1}^{n}p_iJ\mathbb{E}
 \Big(\int_{t}^{t+1}\int_\mathbb{Y}x_i^2(s)\gamma_i^2(u)N(d s,d u)\Big)^{1/2}\\
&\leq \sum_{i=1}^{n}p_iJ\Big(\mathbb{E}\int_{t}^{t+1}
 \int_\mathbb{Y}x_i^2(s)\gamma_i^2(u)N(d s,d u)\Big)^{1/2} \\
&=\sum_{i=1}^{n}p_iJ\Big(\int_\mathbb{Y}\gamma_i^2(u)\nu(d u)
 \mathbb{E}\int_{t}^{t+1}x_i^2(s)d s\Big)^{1/2}\\
&\leq  J\sum_{i=1}^{n}p_i\Big(\int_\mathbb{Y}\gamma_i^2(u)\nu(d u)\Big)^{1/2}
\Big(\mathbb{E}\int_{t}^{t+1}|x(s)|^2d s\Big)^{1/2},
\end{align*}
where $J$ is a positive constant. Moreover, we can derive from
Theorem \ref{thm4} that there exists  positive constants $K_1$, $K_2$ such
that $\limsup_{t\to\infty}\mathbb{E}\int_{t}^{t+1}|x(s)|d s\leq K_1$
and $\limsup_{t\to\infty}\mathbb{E}\int_{t}^{t+1}|x(s)|^2d s\leq K_2$.
Substituting the above inequalities into \eqref{m33} and combining \eqref{m32},
we can see that
\begin{align*}
&\limsup_{t\to\infty}\mathbb{E}\Big(\sup_{t\leq u\leq t+1}V(x(u))\Big)\\
&\leq \frac{K_0}{\alpha}+q_1K_1+q_2K_2
+\Big(3|\bar{p}\sigma|+J\sum_{i=1}^{n}p_i
 \Big(\int_\mathbb{Y}\gamma_i^2(u)\nu(d u)\Big)^{1/2}\Big)K_2^{1/2}.
\end{align*}
Hence there is a positive constant $M$ such that
\[
\mathbb{E}\Big(\sup_{n\leq t\leq n+1}|x(t)|\Big)\leq M ,\quad n=1,2,\dots.
\]
Let $\varepsilon>0$ be arbitrary, by the Chebyshev inequality, we have
\begin{equation*}
\mathbb{P}\big\{\sup_{n\leq t\leq n+1}|x(t)|>n^{1+\varepsilon}\big\}
\leq \frac{M}{n^{1+\varepsilon}} \quad  n=1,2,\dots.
\end{equation*}
Since the series $\sum_{n=1}^{\infty}\frac{M}{n^{1+\varepsilon}}$ converges,
then from the Borel-Cantella lemma \cite{gm12} that there exists a
$n_0:=n_0(\omega)$ such that for almost all $\omega\in\Omega$, whenever
$n\geq n_0$ and $n\leq t\leq n+1$, we have
\begin{equation*}
\sup_{n\leq t\leq n+1}|x(t)|\leq n^{1+\varepsilon}.
\end{equation*}
So
\begin{equation*}
\limsup_{t\to\infty}\frac{\ln|x(t)|}{\ln t}\leq 1+\varepsilon.
\end{equation*}
Letting $\varepsilon\to 0$ leads to the desired assertion.
\end{proof}

\begin{remark} \label{rmk3} \rm
Noting that the limit $\lim_{t\to\infty}\frac{\ln t}{t}=0$,
under the conditions of Theorem \ref{thm5},  we obtain
$\limsup_{t\to\infty}\frac{\ln x_i(t) }{t}\leq 0$, a.s., $i=1,\dots,n$.
\end{remark}

On the other hand, by the positivity of solution of \eqref{m0} and the
 comparison theorem \cite[Theorem 3.1]{m12},  we obtain that
$$
x_i(t)\geq y_i(t),\quad i=1,\dots,n,
$$
where $y_i(t)$ is the solution of \eqref{m41}. According to the analysis
for \eqref{m41} in section 2, we obtain the following results.

\begin{theorem} \label{thm6}
Let {\rm(A1), (A2)} hold and $\theta_i>1$, $i=1,\dots,n$.
If $r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2
+\int_{\mathbb{Y}}\ln(1+\gamma_i(u))\nu(d u)\geq 0$, then for any initial
value $x_0\in \mathbb{R}_+^n$, the solution of \eqref{m0} satisfies
\begin{gather*}
\begin{aligned}
\liminf_{t\to\infty}\frac{1}{t}\int_{0}^{t}x_i^{\theta_i}(s)d s
&\geq \lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}y_i^{\theta_i}(s)d s\\
&=\frac{r_i-\frac{1}{2}\sum_{j=1}^{m}\sigma_{ij}^2
 +\int_{\mathbb{Y}}\ln(1+\gamma_i(u))\nu(d u)}{a_{ii}}, 
\end{aligned}\\
\liminf_{t\to\infty}\frac{\ln x_i(t)}{t}\geq 0 \quad\text{a.s., } i=1,\dots,n.
\end{gather*}
\end{theorem}

Now combining Remark \ref{rmk3} and Theorem \ref{thm6} leads to the following theorem.

\begin{theorem} \label{thm7}
Under the conditions of Theorem \ref{thm6}, for each $i=1,\dots,n,$
$$
\lim_{t\to\infty}\frac{\ln x_i(t)}{t}=0\quad\text{a.s. }
$$
\end{theorem}

\section{Stochastic permanence}

Stochastic permanence is one of the most interesting and important topics.
 In this section, stochastic permanence is studied based
on the results in Section 4. We firstly introduce the definition of
stochastic permanence.

\begin{definition}[\cite{m11}] \label{def2} \rm
If for arbitrary $\varepsilon\in (0,1)$, there are two positive constants
$\beta_1$ and $\beta_2$ such that for any initial data $x_0\in\mathbb{R}_+^n$,
 the solution $x(t)$ of Eq.\eqref{m0} has the property that
$$
\liminf_{t\to\infty}\mathbb{P}\{x_i(t)\geq \beta_1\}\geq 1-\varepsilon,\quad
\liminf_{t\to\infty}\mathbb{P}\{x_i(t)\leq \beta_2\}
\geq 1-\varepsilon, \quad 1\leq i\leq n,
$$
then \eqref{m0} is said to be stochastically permanent.
\end{definition}

\begin{theorem} \label{thm8}
Under the conditions of Theorem \ref{thm2}, if there exists a positive constant $\alpha$,
such that
$$
r_i-\frac{3+\alpha}{2}\sum_{j=1}^{m}\sigma_{ij}^2-\int_{\mathbb{Y}}
\Big[\frac{1}{(2+\alpha) (1+\gamma_i(u))^{2+\alpha}}-\frac{1}{2+\alpha}\Big]
\nu(d u)>0,
$$
then, for the case $1\leq\theta_i\leq 2+\alpha$, Equation \eqref{m0}
is stochastically permanent.
\end{theorem}
\begin{proof}

Define $y_i=1/x_i$ for $x_i>0$, applying It\^o's formula, we obtain
\begin{align*}
d y_i^{2+\alpha}
&=d (\frac{1}{x_i})^{2+\alpha}\\
&=-(2+\alpha)(\frac{1}{x_i})^{\alpha+3}x_i(r_i-a_{ii}x_i^{\theta_i}
 -\sum_{j\neq i}^{n}a_{ij}x_j)d t\\
&\quad +\frac{2+\alpha}{2}(\alpha+3)(\frac{1}{x_i})^{\alpha+4}x_i^2
 \sum_{j=1}^{m}\sigma_{ij}^2d t
 -(2+\alpha)(\frac{1}{x_i})^{\alpha+3}x_i\sum_{j=1}^{m}
 \sigma_{ij}d B_j
\\
&\quad+\int_\mathbb{Y}[\frac{1}{(x_i+x_i\gamma_i(u))^{2+\alpha}}
 -\frac{1}{x_i^{2+\alpha}}]N(d t, d u)
\\
&\leq  (2+\alpha)(\frac{1}{x_i})^{\alpha}\Big\{-\frac{1}{x_i^2}\Big[r_i
 -\frac{\alpha+3}{2}\sum_{j=1}^{m}\sigma_{ij}^2
 -\int_\mathbb{Y}(\frac{1}{(2+\alpha)(1+\gamma_i(u))^{2+\alpha}}
\\
&\quad -\frac{1}{2+\alpha})\nu(d u)\Big]+\frac{a_{ii}}{x_i^{2-\theta_i}}\Big\}d t
-(2+\alpha)(\frac{1}{x_i})^{\alpha+2}\sum_{j=1}^{m}\sigma_{ij}d B_j
\\
&\quad +\int_\mathbb{Y}[\frac{1}{(x_i+x_i\gamma_i(u))^{2+\alpha}}
 -\frac{1}{x_i^{2+\alpha}}]\widetilde{N}(d t, d u)
\\
&=  (2+\alpha)y_i^{\alpha}\Big\{-y_i^2\Big[r_i-\frac{\alpha+3}{2}
 \sum_{j=1}^{m}\sigma_{ij}^2
  -\int_\mathbb{Y}(\frac{1}{(2+\alpha)(1+\gamma_i(u))^{2+\alpha}}\\
&\quad  -\frac{1}{2+\alpha})\nu(d u)\Big]+a_{ii}y_i^{2-\theta_i}\Big\}d t
 -(2+\alpha)y_i^{\alpha+2}\sum_{j=1}^{m}\sigma_{ij}d B_j
\\
&\quad +\int_\mathbb{Y}y_i^{2+\alpha}[\frac{1}{(1+\gamma_i(u))^{2+\alpha}}-1]
 \widetilde{N}(d t, d u).
\end{align*}
Choose a sufficiently small positive $\zeta$  such that
$$
r_i-\frac{\alpha+3}{2}\sum_{j=1}^{m}\sigma_{ij}^2
-\int_\mathbb{Y}\Big(\frac{1}{(2+\alpha)(1+\gamma_i(u))^{2+\alpha}}
-\frac{1}{2+\alpha}\Big)\nu(d u)>\frac{\zeta}{2+\alpha}.
$$
Define $V=e^{\zeta t}y_i^{2+\alpha}$, using It\^o's formula results in
\begin{align*}
d V
&\leq (2+\alpha) e^{\zeta t}y_i^{\alpha}
\Big\{-y_i^2\Big[r_i-\frac{\alpha+3}{2}\sum_{j=1}^{m}\sigma_{ij}^2
 -\int_\mathbb{Y}(\frac{1}{(2+\alpha)(1+\gamma_i(u))^{2+\alpha}}
\\
&\quad  -\frac{1}{2+\alpha})\nu(d u)\Big]+a_{ii}y_i^{2-\theta_i}\Big\}d t
+\zeta e^{\zeta t}y_i^{2+\alpha} d t
\\
&\quad -(2+\alpha)e^{\zeta t}y_i^{\alpha+2}\sum_{j=1}^{m}
 \sigma_{ij}d B_j+e^{\zeta t}\int_\mathbb{Y}y_i^{2+\alpha}
 [\frac{1}{(1+\gamma_i(u))^{2+\alpha}}-1]\widetilde{N}(d t, d u)
\\
&=: e^{\zeta t}F(y_i)d t-(2+\alpha) e^{\zeta t} y_i^{2+\alpha}
 \sum_{j=1}^{m}\sigma_{ij}d B_j\\
&\quad +e^{\zeta t}\int_\mathbb{Y}y_i^{2+\alpha}
[\frac{1}{(1+\gamma_i(u))^{2+\alpha}}-1]
 \widetilde{N}(d t,d u),
\end{align*}
where
\begin{align*}
F(y_i)&=(2+\alpha) y_i^{\alpha}\Big\{-y_i^2\Big[r_i
-\frac{\alpha+3}{2}\sum_{j=1}^{m}\sigma_{ij}^2
 -\int_\mathbb{Y}(\frac{1}{(2+\alpha)(1+\gamma_i(u))^{2+\alpha}}\\
&\quad -\frac{1}{2+\alpha})\nu(d u)-\frac{\zeta}{2+\alpha}\Big]
+a_{ii}y_i^{2-\theta_i}\Big\}
\end{align*}
has an upper positive bound, say $K$.
Integrating from 0 to $t$ and taking expectations, we obtain
\begin{equation*}
\mathbb{E}[e^{\zeta t}y_i^{2+\alpha}]\leq y_i^{2+\alpha}(0)
+\frac{K}{\zeta}(e^{\zeta t}-1).
\end{equation*}
Therefore
\begin{equation*}
\limsup_{t\to\infty}\mathbb{E}[x_i^{-(\alpha+2)}(t)]
\leq \frac{K}{\zeta}=:K_0.
\end{equation*}
For any given $\varepsilon>0$, set
$\beta_1=\varepsilon^{1/(2+\alpha)}/K_0^{1/(2+\alpha)}$,
by the Chebyshev inequality, we obtain
\begin{align*}
\limsup_{t\to\infty}\mathbb{P}\big\{x_i(t)<\beta_1\big\}
&= \limsup_{t\to\infty}\mathbb{P}\big\{x_i^{-(2+\alpha)}(t)>\beta_1^{-(2+\alpha)}
\big\}\\
&\leq\limsup_{t\to\infty}\beta_i^{2+\alpha}
\mathbb{E}[x_i^{-(2+\alpha)}(t)]=\varepsilon.
\end{align*}
Hence $\liminf_{t\to\infty}\mathbb{P}\{x_i(t)\geq\beta_1\}\geq 1-\varepsilon $,
$i=1,\dots,n$.

On the other hand, as an application of Theorem \ref{thm3}, and
the Chebyshev inequality, we can easily show that for arbitrary
$\varepsilon\in (0,1)$, there is a positive constant $\beta_2$
such that for any initial data $x_0\in\mathbb{R}_+^n$,
$\liminf_{t\to\infty}\mathbb{P}\{x_i(t)\leq \beta_2\}\geq 1-\varepsilon$,
$ 1\leq i\leq n$.
Therefore, \eqref{m0} is stochastically permanent.
\end{proof}

\section{Example and numerical simulations}

Consider the two-species stochastic Gilpin-Ayala mutualism system with jumps
\begin{equation}\label{m7}
\begin{gathered}
d x_1=x_1(r_1-a_{11}x_1^{\theta_1}-a_{12}x_2)d t
 +\sigma_{11}d B_1(t)+\sigma_{12}d B_2(t)
 +\int_{\mathbb{Y}}\gamma_1(u)x_1N(d t,d u),\\
d x_2=x_2(r_2-a_{22}x_2^{\theta_2}-a_{21}x_1)d t
 +\sigma_{21}d B_1(t)+\sigma_{22}d B_2(t)
 +\int_{\mathbb{Y}}\gamma_2(u)x_2N(d t,d u).
\end{gathered}
\end{equation}

In Figure \ref{fig1}, we choose $r_1 = 0.06$, $r_2 = 0.05$, $a_{11} = 0.08$,
$a_{22} = 0.04$, $a_{12} =a_{21}=-0.005$, $\theta_1=\theta_2=1.01$,
$\sigma_{11}=0.2$, $\sigma_{22}=0.05$, $\sigma_{12}=\sigma_{21}=0$,
$\gamma_1(u)=0.2$, $\gamma_2(u)=0.24$, $x_1(0)=1.1$, $x_2(0)=1.5$,
$\mathbb{Y}=(0,\infty)$, $\lambda(\mathbb{Y})=1$.
Since $a_{11}a_{22}-a_{12}a_{21}>0$, then (A1) and (A2) hold, so
\eqref{m7} has a unique global positive solution for any positive initial
value by Theorem \ref{thm2} and pth moment of the solution of \eqref{m7} is
asymptotic bounded, see Figure \ref{fig1}. Moreover,
\begin{gather*}
r_1-\frac{\sigma_{11}^2+\sigma_{12}^2}{2}
+\int_{\mathbb{Y}}\ln(1+\gamma_1(u))\nu(d u)=0.22>0,
\\
r_2-\frac{\sigma_{21}^2+\sigma_{22}^2}{2}
 +\int_{\mathbb{Y}}\ln(1+\gamma_2(u))\nu(d u)=0.26>0.
\end{gather*}
Then in view of Theorem \ref{thm6}, we obtain $\frac{\ln x_i(t)}{t}\to 0$,
 $i=1,2$, Figure \ref{fig1} confirms these.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig1} % zxhliyi.eps 
\end{center}
\caption{Solutions of \eqref{m7} for $r_1 = 0.06$,
$r_2 = 0.05$, $a_{11} = 0.08$, $a_{22} = 0.04$,
$a_{12} =a_{21}=-0.005$, $\theta_1=\theta_2=1.01$,
$\sigma_{11}=0.2$, $\sigma_{22}=0.05$,
$\sigma_{12}=\sigma_{21}=0$ , $\gamma_1(u)=0.2$,
$\gamma_2(u)=0.24$, $x_1(0)=1.1$, $x_2(0)=1.5$,
$\mathbb{Y}=(0,\infty)$, $\lambda(\mathbb{Y})=1$}
\label{fig1}
\end{figure}

In Figure \ref{fig2}, we choose $r_1 = 0.5$, $r_2 = 0.2$, $a_{11} =a_{22} = 0.9$,
$a_{12} =a_{21}=-0.05$, $\theta_1=\theta_2=1.01$, $\sigma_{11}=0.05$,
$\sigma_{22}=0.1$, $\sigma_{12}=\sigma_{21}=0$, $\gamma_1(u)=0.2$,
$\gamma_2(u)=0.12$, $x_1(0)=0.1$, $x_2(0)=0.8$, $\mathbb{Y}=(0,\infty)$,
 $\lambda(\mathbb{Y})=1$. Since $a_{11}a_{22}-a_{12}a_{21}>0$,
then  (A1) and (A2) hold.
In Theorem \ref{thm8}, choose $\alpha=1$, by simple calculation, we have
\begin{gather*}
r_1-\frac{3+\alpha}{2}\sum_{j=1}^{m}\sigma_{1j}^2
-\int_{\mathbb{Y}}\Big[\frac{1}{(2+\alpha) (1+\gamma_1(u))^{2+\alpha}}
-\frac{1}{2+\alpha}\Big]\nu(d u)=0.64>0,\\
r_2-\frac{3+\alpha}{2}\sum_{j=1}^{m}\sigma_{2j}^2
 -\int_{\mathbb{Y}}\Big[\frac{1}{(2+\alpha) (1+\gamma_2(u))^{2+\alpha}}
 -\frac{1}{2+\alpha}\Big]\nu(d u)=0.28>0\,.
\end{gather*}
Theorem \ref{thm8} tells us that \eqref{m7} is stochastically permanent,
and Figure \ref{fig2} confirms this.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig2} % lier.eps
\end{center}
\caption{Solutions of \eqref{m7} for $r_1 = 0.5$, $r_2 = 0.2$,
$a_{11} =a_{22} = 0.9$, $a_{12} =a_{21}=-0.05$, $\theta_1=\theta_2=1.01$,
$\sigma_{11}=0.05$, $\sigma_{22}=0.1$, $\sigma_{12}=\sigma_{21}=0$,
$\gamma_1(u)=0.2$, $\gamma_2(u)=0.12$, $x_1(0)=0.1$, $x_2(0)=0.8$,
$\mathbb{Y}=(0,\infty)$, $\lambda(\mathbb{Y})=1$}
\label{fig2}
\end{figure}

\subsection*{Conclusions}
An stochastic Gilpin-Ayala mutualism model with jumps has been
studied in this article. The high nonlinearity of Gilpin-Ayala model
and Poisson jumps make the problem difficult. Sufficient criteria
for the existence of global positive solution, stochastically
ultimate boundedness and stochastic permanence are derived for
the n-dimensional model by analysis of Lyapunov functions which has been
used by many authors. We also investigate asymptotic pathwise estimation.
The simulation results verify the effectiveness of the proposed results.


\subsection*{Acknowledgements}
The authors want thank the editor and referee for their important
and valuable comments. The authors also want thank the National Natural
Science Foundation of  China for their economic support through the
grants: 11171081, 11171056, 11126219, 11001032, 11226254,
 and 11101183.

\begin{thebibliography}{00}

\bibitem{m03} J. F. Addicott, H. I. Freedman;
On the structure and stability of mutualistic systems:
analysis of predatorprey and competition models as modified by
the action of a slow-growing mutualist, Theor. Popul. Biol. 26 (1984), 320-339.

\bibitem{gs09} D. Applebaum;
L\'evy Processes and Stochastic Calculus, 2nd ed., Cambridge
University Press, 2009.

\bibitem{gs07} J. Bao, X. Mao, G. Yin, C. Yuan;
Competitive Lotka-Volterra population dynamics with jumps,
Nonlinear Anal. 74 (2011), 6601-6616.

\bibitem{gs08} J. Bao, C. Yuan;
Stochastic population dynamics driven by L\'evy noise,
J. Math. Anal. Appl. 391 (2012), 363-375.

\bibitem{ga06} F. Chen, C. Shi;
 Global attractivity in an almost periodic multi-species nonlinear
ecological model, Appl. Math. Comput. 180 (2006), 376-392.

\bibitem{m09} S. Cheng;
Stochastic population systems, Stoch. Anal. Appl. 27 (2009), 854-874.

\bibitem{gm06} M. Fan, K. Wang;
 Global periodic solutions of a generalized n-species Gilpin-Ayala
 competition model, Comput. Math. Appl. 40 (2000), 1141-1151.
\bibitem{gm21} M. E. Gilpin, F. J. Ayala;
Global models of growth and competition, Proc. Nat. Acad. Sci.
 USA 70 (1973), 3590-3593.

\bibitem{m01} B. S. Goh;
 Stability in models of mutualism, Amer. Natural. 113 (1979), 261-275.

\bibitem{m02} K. Golpalsamy;
 Stability and Oscillations in Delay Differential Equations of Population
 Dynamics, Kluwer Academic, Dordrecht, 1992.

\bibitem{m08} L. Hung;
Stochastic delay population systems, Appl. Anal. 88 (2009), 1303-1320.

\bibitem{gm19} R. Z. Khasminskii;
 Stochastic Stability of Differential Equations, Sijthoff and Noordhoff,
 Alphen aan den Rijn, 1980.

\bibitem{gm16} Y. Kuang;
 Delay Differential Equations with Applications in Population Dynamics,
  Academic Press, Boston, 1993.

\bibitem{gs13} H. Kunita;
 It\^o's stochastic calculus: Its surprising power for applications,
Stochastic Processses and their Applications, 120 (2010), 622-652.

\bibitem{gm03} X. Li, X. Mao;
 Population dynamical behavior of non-autonomous Lotka-Volterra
competitive system with random perturbation, Discrete Contin. Dyn. Syst.
24 (2009), 523-545.

\bibitem{gm05} B. Lian, S. Hu;
 Asymptotic behaviour of the stochastic Gilpin-Ayala competition models,
J. Math. Anal. Appl. 339 (2008), 419-428.

\bibitem{ga07} X. Liao, S. Zhou, Y. Chen;
 On permanence and global stability in a general Gilpin-Ayala
competition predator-prey discrete system,
Appl. Math. Comput. 190 (2007), 500-509.

\bibitem{gs14} R. Lipster;
A strong law of large numbers for local martingales,
Stochastics, 3 (1980), 217-228.

\bibitem{m10} M. Liu, K. Wang;
 Population dynamical behavior of Lotka-Volterra cooperative systems with
random perturbations, Discrete Contin. Dyn. Syst. 33 (2013), 2495-2522.

\bibitem{m11} M. Liu, K. Wang;
Analysis of a stochastic autonomous mutualism model,
J. Math. Anal. Appl. 402 (2013), 392-403.

\bibitem{m04} Z. Lu, Y. Takeuchi;
Permanence and global stability for cooperative Lotka-Volterra diffusion
systems, Nonlinear Anal. 19 (1992), 963-975.

\bibitem{gm01} X. Mao, G. Marion, E. Renshaw;
 Environmental Brownian noise suppresses explosions in population dynamics,
 Stoch. Proc. Appl. 97 (2002), 95-110.

\bibitem{gm12} X. Mao, C. Yuan;
Stochastic Differential Equations With Markovian Switching,
Imperial College Press, London, 2006.

\bibitem{m05} J. Pan, Z. Jin, Z. Ma;
Thresholds of survival for an n-dimensional Volterra mutualistic system
in a polluted environment, J. Math. Anal. Appl. 252 (2000), 519-531.

\bibitem{m12} S. Peng, X. Zhu;
Necessary and sufficient condition for comparison theorem of
1-dimensional stochastic differential equations,
Stoch. Proc. Appl. 116 (2006), 370-380.

\bibitem{m06} H. L. Smith;
On the asymptotic behavior of a class of deterministic models of
cooperating species, SIAM J. Math. Anal. 46 (1986), 368-375.

\bibitem{gs12} R. Situ; Theory of Stochastic Differential Equations with
Jumps and Applications, Springer, 2005.

\bibitem{m07} A. Tineo;
On the asymptotic behavior of some population models, II,
 J. Math. Anal. Appl. 197 (1996), 249-258.

\bibitem{gm17} M. Vasilova, M. Jovanovi\'c;
Dynamics of Gilpin-Ayala competition model with random perturbation,
Filomat 24 (2010), 101-113.

\bibitem{gm18} M. Vasilova, M. Jovanovi\'c;
Stochastic Gilpin-Ayala competition model with infinite delay,
Appl. Math. Comput. 217 (2011), 4944-4959.

\end{thebibliography}

\end{document}