1. Introduction

This paper considers the long-term behavior of a system that results from impulsive and stochastic perturbations of a deterministic dynamical system with delay. This type of dynamics occurs naturally for example in the modeling of biological systems such as a growing bacterial colony, modeled deterministically, from which samples of random size are drawn regularly in an experiment or a marine ecosystem from which fish are harvested with a net [2], or the modeling of stochasticity in gene regulatory networks with feedback through a cell's signaling system [19], etc.

The associated mathematical model with stochastic perturbations

$dx(t)=x(t)(r(t)-a(t)x(t))dt+\sigma(t)x(t)dB(t),$

(1)

where $x(t)$ is the population size, $B(t)$ is a standard Brownian motion, $r(t), \ a(t)$ and $\sigma(t)$ are continuous bounded functions on $[0, \infty)$, has been considered by many authors [3,6,11,12,15,16]. Recently, Liu and Wang [17] studied a stochastic logistic system with impulsive effects which takes the following form

$\left\{\begin{array}{l} dx(t)=x(t)(r(t)-a(t)x(t))dt+\sigma(t)x(t)dB(t), ~t\neq t_k, \\ x(t_k^+)-x(t_k)=b_kx(t_k), ~t=t_k, k\in N, \end{array}\right.$

(2)

where $N$ denotes the set of positive integers. The time sequence $\{t_k\}$ is strictly increasing such that $\lim_{k\rightarrow\infty}t_k=+\infty$. The term $b_k$ is impulsive perturbations at the moments of time $t_k$, and satisfies that $1+b_k>0$ for all $k\in N$. The authors obtained the sufficient conditions for extinction, non-persistence in the mean, weak persistence, persistence in the mean and stochastic permanence of the system.

As is well known, time delay always exists in the evolutionary processes of the population (e.g., resource regeneration times, maturation periods, feeding times, reaction times, etc) and it can cause the oscillatory or even unstable phenomena (see monographs [8,14]). Therefore, it seems more realistic to consider model (2) with delay. As far as we know, few investigations have been made on such type system. Based on model (2), we propose the following model:

$\left\{\begin{array}{l} dx(t)=x(t)\big(r-ax(t)-c\int_0^\infty f(s)x(t-s)ds\big)dt+\sigma(t)x(t)dB(t), ~t\neq t_k, \\ x(t_k^+)-x(t_k)=b_kx(t_k), ~t=t_k, k\in N, \end{array}\right.$

(3)

where $c$ stands for the effect power of the past history. $f(s)$, called the delay kernel, is a weighting factor which indicates how much emphasis should be given to the size of the population at earlier times to determine the present effect on resource availability[24]. Usually, we use the Gamma distribution delay kernel

$f_\alpha^n(s)=\frac{\alpha^n s^{n-1}e^{-\alpha s}}{(n-1)!},$

where $\alpha>0$ is a constant, $n$ an integer, with the average delay $T=n/\alpha$. The corresponding version of model (3) with a Gamma distribution delay kernel can be written as

$\left\{\begin{array}{l} dx(t)=x(t)\left[r-ax(t)-c\int_0^\infty f_\alpha^n(s)x(t-s)ds\right]dt+\sigma x(t)dB(t), ~t\neq t_k, \\ x(t_k^+)-x(t_k)=b_kx(t_k), ~t=t_k, ~k\in N. \end{array}\right.$

(4)

To perform a thorough analysis on model (3) with a general delay kernel function is very difficult. So in this paper we mainly devote our attention to the investigation on the dynamics of its special case model (4).

Notice that in the absence of impulsive and random perturbations, it is obvious that model (4) has a stable positive equilibrium $x^*=\frac{r}{a+c}$. However, under impulsive and random perturbations, the equilibrium does not exist. That is, the steady state of model (4) can no longer be represented by a single point. Therefore, we turn our attention to the study of the stationary distribution of the system. To the best of our knowledge, though some significant progress has been made in the techniques and methods of determining the existence of the stationary distribution and the stability of the density for stochastic differential equations [20,9], unfortunately, few works have been performed for the corresponding problems of impulsive stochastic differential equations up to now. So, our work in this paper may be considered as a first attempt to the investigation on the stationary distribution of the impulsive stochastic differential equations.

The organization of this paper is as follows. In the next section, we show the existence and uniqueness of a global positive solution of model (4). Then, in section 3, we present our main results: we first carry out the global attractivity analysis of the model; then using the Hasminskii's method and constructing Lyapunov function, we prove the existence of the stationary distribution of the model; we also perform an extinction analysis of the model. Finally, some discussions and numerical simulations are presented in section 4.

---

2. Existence and uniqueness of a global positive solution

In this section, we show the uniqueness and global existence of a positive solution of model (4) for any positive initial value. Using the linear chain trick, if we define

$y_n(t)=\int^{\infty}_0f_\alpha^n(s)x(t-s)ds,$

then model (4) is equivalent to the following system

$\left\{\begin{array}{l}\left.\begin{array}{l} dx(t)=x(t)(r-ax(t)-c y_n(t))dt+\sigma x(t)dB(t), \\ dy_1(t)=-\alpha (y_1(t)-x(t))dt, \\ dy_i(t)=-\alpha (y_i(t)-y_{i-1}(t))dt, \ i=2, 3, \ldots, n, \end{array}\right\} ~t\neq t_k, \\ x(t_k^+)-x(t_k)=b_kx(t_k), ~t=t_k, k\in N. \end{array}\right.$

(5)

Therefore, to study the dynamics of model (4), we need only to consider system (5).

To begin with, we consider a general $d$-dimensional impulsive stochastic differential equation:

$\left\{\begin{array}{l} dX(t)=F(t, X(t))dt+G(t, X(t))dB(t), ~t\neq t_k, ~ k\in N, \\ X(t_k^+)-X(t_k)=b_kX(t_k), ~k\in N \end{array}\right.$

(6)

with initial condition $X(0)=X_0\in \mathbb{R}^d$. $B(t)$ denotes $m$-dimensional standard Brownian motions defined on the probability space $(\Omega, \mathscr{F}, \{\mathscr{F}_t\}_{t\geq0}, P)$, which is a complete probability space with a filtration $\{\mathscr{F}_t\}_{t\geq0}$ satisfying the usual conditions (i.e. it is right continuous and $\mathscr{F}_0$ contains all P-null sets). Denote the differential operator $L$ associated with Eq. (6) by

$L=\frac{\partial}{\partial t}+\sum\limits_{i=1}^dF_i(t, X)\frac{\partial}{\partial X_i}+\frac{1}{2}\sum\limits_{i, j=1}^d\big[G^T(t, X)G(t, X)\big]_{ij}\frac{\partial^2}{\partial X_i\partial X_j}$

(7)

for $t\in[t_{k-1}, t_k)$. Let $PC^{1, 2}([t_{k-1}, t_k)\times\mathbb{R}^d; \mathbb{R}_+)$ denote the family of all nonnegative functions $V(t, X)$ on $[t_{k-1}, t_k)\times\mathbb{R}^d$ which are continuous once differentiable in $t$ and twice differentiable in $X$. If $L$ acts on such $V(t, X)$, then

$LV(t, X)=V_t(t, X)+V_X(t, X)F(t, X)+\frac{1}{2}trace\big[G^T(t, X)V_{XX}(t, X)G(t, X)\big].$

Definition 2.1. (See [17]) A stochastic process $X(t)=(X_1(t), X_2(t)..., X_d(t))^T, t\in R_+$, is said to be a solution of Eq. (6) if

(ⅰ) $X(t)$ is $\mathcal {F}_t$-adapted and is continuous on $(0, t_1)$ and each interval $(t_k, t_{k+1})\subset R_+, k\in N$; $F(t, X(t))\in \mathcal {L}^1(R_+;R^d), G(t, X(t))\in \mathcal {L}^2(R_+;R^d)$, where $\mathcal {L}^k(R_+;R^d)$ is all $R^n$ valued measurable $\{\mathcal {F}_t\}$-adapted processes $f(t)$ satisfying $\int_0^T|f(t)|^kdt<\infty$ almost surely for every $T>0$;

(ⅱ) for each $t_k, k\in N, X(t_k^+)=\lim_{t\rightarrow t_k^+}X(t)$ and $X(t_k^-)=\lim_{t\rightarrow t_k^-}X(t)$ exist and $X(t_k)=X(t_k^-)$ with probability one;

(ⅲ) for almost all $t\in[0, t_1]$, $X(t)$ obeys the integral equation

$X(t)=X(0)+\int_0^tF(s, X(s))ds+\int_0^tG(s, X(s))dB(s)$

and for almost all $t\in(t_k, t_{k+1}], k\in N, X(t)$ obeys the integral equation

$X(t)=X(t_k^+)+\int_0^tF(s, X(s))ds+\int_0^tG(s, X(s))dB(s).$

Moreover, $X(t)$ satisfies the impulsive conditions at each $t=t_k, k\in N$ with probability one.

We now establish a fundamental lemma to reduce the dynamics of nonlinear stochastic differential system under impulsive perturbation to the corresponding problems of a nonlinear stochastic differential system without impulses.

Lemma 2.1. Consider the following stochastic differential equation:

$dY(t)=f(t, Y(t))dt+g(t, Y(t))dB(t),$

(8)

where

$\begin{align} f(t, Y(t))=\prod\limits_{0<t_k<t}(1+b_k)^{-1}F\Big(t, \prod\limits_{0<t_k<t}(1+b_k)Y(t)\Big), \notag \\ g(t, Y(t))=\prod\limits_{0<t_k<t}(1+b_k)^{-1}G\Big(t, \prod\limits_{0<t_k<t}(1+b_k)Y(t)\Big) \end{align}$

with the initial $Y(0)=X_0$.

( $i$ ) If $Y(t)$ is a solution of Eq. (8), then $X(t)=\prod_{0<t_k<t}(1+b_k)Y(t)$ is a solution of Eq. (6) on $[0, \infty)$;

( $ii$ ) if $X(t)$ is a solution of Eq. (6), then $Y(t)=\prod_{0<t_k<t}(1+b_k)^{-1}X(t)$ is a solution of Eq. (8) on $[0, \infty).$

Proof. The proof is inspired by the method of Yan [26]. Fist, we prove (i). Let $Y(t)$ be a possible solution to Eq. (8), it is easy to see that $X(t)=\prod_{0<t_k<t}(1+b_k)Y(t)$ is continuous on each interval $(t_k, t_{k+1})\subset [0, \infty), k\in N$, and for any $t\neq t_k$,

$\begin{align} dX(t)&=d\Big(\prod\limits_{0<t_k<t}(1+b_k)Y(t)\Big) \notag \\ &=\prod\limits_{0<t_k<t}(1+b_k)f(t, Y(t))dt+\prod\limits_{0<t_k<t}(1+b_k)g(t, Y(t))dB(t) \notag \\ &=F\Big(t, \prod\limits_{0<t_k<t}(1+b_k)Y(t)\Big)dt+G\Big(t, \prod\limits_{0<t_k<t}(1+b_k)Y(t)\Big)dB(t) \notag \\ &=F(t, X(t))dt+G(t, X(t))dB(t). \end{align}$

So, $X(t)$ satisfies Eq. (6) for almost everywhere on $[0, \infty)/\{t_k\}, k=1, 2, ...$

On the other hand, for each $k\in N$ and $t_k\in [0, \infty)$,

$\begin{align} X(t_k^+) &=\lim\limits_{t\rightarrow t_k^+}\prod\limits_{0<t_j<t}(1+b_j)Y(t)\notag \\ &=\prod\limits_{0<t_j\leq t_k}(1+b_j)Y(t_k^+) \notag \\ &=(1+b_k)\prod\limits_{0<t_j< t_k}(1+b_j)Y(t_k) \notag \\ & =(1+b_k)X(t_k) \notag \end{align}$

and

$\begin{align} X(t_k^-)&=\lim\limits_{t\rightarrow t_k^-}\prod\limits_{0<t_j<t}(1+b_j)Y(t) \notag\\ &=\prod\limits_{0<t_j< t_k}(1+b_j)Y(t_k^-) \notag\\ &=\prod\limits_{0<t_j< t_k}(1+b_j)Y(t_k) \notag\\ &=X(t_k), \end{align}$

which implies that $X(t)$ is the solution of Eq. (6).

We now prove (ⅱ). Since $X(t)$ is a solution of Eq. (6) and is continuous on each interval $(t_k, t_{k+1})\subset [0, +\infty)$, thus $Y(t)=\prod_{0<t_k<t}(1+b_k)^{-1}X(t)$ is continuous on $(t_k, t_{k+1})\subset [0, +\infty)$. What's more, we have that, for any $k\in N$ and $t_k\in [0, +\infty)$,

$\begin{align} Y(t_k^+) & =\lim\limits_{t\rightarrow t_k^+}\prod\limits_{0<t_j<t}(1+b_j)^{-1}X(t) =\prod\limits_{0<t_j\leq t_k}(1+b_j)^{-1}X(t_k^+) \notag\\ & =\prod\limits_{0<t_j<t_k}(1+b_j)^{-1}X(t_k) =Y(t_k) \notag \end{align}$

and

$Y(t_k^-)=\lim\limits_{t\rightarrow t_k^-}\prod\limits_{0<t_j<t}(1+b_j)^{-1}X(t) =\prod\limits_{0<t_j<t_k}(1+b_j)^{-1}X(t_k^-)=Y(t_k).$

Therefore, $Y(t)$ is continuous on $[0, +\infty)$. Similar to the proof in the case (i), we can easily check that $Y(t)$ is the solution of Eq. (8) on $[0, +\infty)$. This completes the proof of the Lemma.

It should be emphasized that we consider Eqs. (6) and (8) on the same probability space. For Eq. (8), the associated differential operator $L$ is the same as that defined in (7) for $t\in[0, \infty)$ rather than for $t\in [t_{k-1, t_k})$.

Now we turn to show the existence and uniqueness of the solution of system (5) by employing Lemma 2.1. Let $x(t)=\prod_{0<t_k<t}(1+b_k)u(t)$ and $y_i(t)=v_i(t), \ i=1, 2, \ldots, n$, then system (5) becomes

$\left\{\begin{array}{l} du(t)=u(t)\Big(r-a\prod_{0<t_k<t}(1+b_k)u(t)-c v_n(t)\Big)dt+\sigma u(t)dB(t), \\ dv_1(t)=-\alpha\Big(v_1(t)-\prod_{0<t_k<t}(1+b_k)u(t)\Big)dt, \\ dv_i(t)=-\alpha(v_i(t)-v_{i-1}(t))dt, \ i=2, 3, \ldots, n. \end{array}\right.$

(9)

By Lemma 2.1, if $(u(t), v_1(t), \ldots, v_n(t))$ is a solution of system (9), then $(x(t), y_1(t),$ $\ldots, y_n(t))$ is a solution of system (5). So, in the following, we only need to show the global existence of the positive solution of system (9).

Lemma 2.2. There is a unique positive global solution $(u(t), v_1(t), \ldots, v_n(t))$ of system (9) a.s. for any initial value $(u_0, v_{10}, \ldots, v_{n0})\in R^{n+1}_+$.

Proof. Consider the system

$\left\{\begin{array}{l} d\bar{u}(t)=\Big(r-\frac{1}{2}\sigma^2-a\prod\limits_{0<t_k<t}(1+b_k)e^{\bar{u}(t)}-c e^{\bar{v}_n(t)}\Big)dt+\sigma dB(t), \\ d\bar{v}_1(t)=-\alpha\Big(1-\frac{\prod\limits_{0<t_k<t}(1+b_k) e^{\bar{u}(t)}}{e^{\bar{v}_1(t)}}\Big)dt, \\ d\bar{v}_i(t)=-\alpha\Big(1-\frac{ e^{\bar{v}_{i-1}(t)}}{e^{\bar{v}_i(t)}}\Big)dt, \ i=2, 3, \ldots, n \end{array}\right.$

(10)

with initial value $\bar{u}(0)=\ln u_0, ~\bar{v}_i(0)=\ln v_{i0}, \ i=1, 2, \ldots, n$. Obviously, the coefficients of system (10) satisfy the local Lipschitz condition, then there is a unique local solution $(\bar{u}(t), \bar{v}_1(t), \ldots, \bar{v}_n(t))$ on $t\in [0, \tau_e)$, where $\tau_e$ is the explosion time. By Itô's formula, it is easy to see

$(u(t), v_1(t), \ldots, v_n(t))=(e^{\bar{u}(t)}, e^{\bar{v}_1(t)}, \ldots, e^{\bar{v}_n(t)})$

is the unique positive local solution of system (9) with initial value $(u_0, v_{10}, \ldots, v_{n0}) \in R^{n+1}_+$.

We now show the solution of system (9) is global, i.e., $\tau_e=\infty$. Since the solution is positive for $t\in[0, \tau_e)$, we have

$du(t)\leq u(t)\Big(r-a\prod\limits_{0<t_k<t}(1+b_k)u(t)\Big)dt+\sigma u(t)dB(t).$

Let

$\phi(t)=\frac{e^{(r-\frac{\sigma^2}{2})t+\sigma B(t)}} {\frac{1}{u_0}+a\prod_{0<t_k<t} (1+b_k)\int^t_0e^{(r-\frac{\sigma^2}{2})s+\sigma B(s)}ds},$

then $\phi(t)$ is the unique solution of the equation

$\left\{\begin{array}{l} d\phi(t)=\phi(t)\Big(r-a\prod_{0<t_k<t}(1+b_k)\phi(t)\Big)dt+\sigma u(t)dB(t), \\ \phi(0)=u_0, \end{array}\right.$

and by the comparison theorem for the stochastic equation, yields

$u(t)\leq\phi(t), t\in[0, \tau_e), \ a.s.$

Besides, we can get

$dv_1(t)\leq\Big(-\alpha v_1(t)+\alpha\prod\limits_{0<t_k<t}(1+b_k)\phi(t)\Big)dt.$

Obviously,

$\varphi_1(t)=v_{10}e^{-\alpha t}+\alpha\prod\limits_{0<t_k<t}(1+b_k)\int_0^t\phi(s)e^{-\alpha(t-s)}ds$

is the solution to the equation

$\left\{\begin{array}{l} d\varphi_1(t)=\Big(-\alpha\varphi_1(t)+\alpha\prod_{0<t_k<t}(1+b_k)\phi(t)\Big)dt, \\ \varphi_1(0)=v_{10} \end{array}\right.$

and $v_1(t)\leq\varphi_1(t), t\in[0, \tau_e), \ a.s.$ Similarly, let $\varphi_i, \ i=2, 3, \ldots, n$ be the solution to the equation

$\left\{\begin{array}{l} d\varphi_i(t)=-\alpha(\varphi_i(t)-\varphi_{i-1}(t))dt, \\ \varphi_i(0)=v_{i0}. \end{array}\right.$

Then we have $v_i(t)\leq \varphi_i(t)$, where $\varphi_i(t)=v_{i0}e^{-\alpha t}+\alpha\int_0^t\varphi_{i-1}(s)e^{-\alpha(t-s)}ds$.

On the other hand,

$du(t)\geq u(t)\Big(r-a\prod\limits_{0<t_k<t}(1+b_k)u(t)-c\varphi_n(t)\Big)dt+\sigma u(t)dB(t).$

It follows that

$u(t)\geq\tilde{\phi}(t), \ t\in[0, \tau_e), \ a.s.$

where

$\tilde{\phi}(t)=\frac{e^{(r-\frac{\sigma^2}{2})t-c\int_0^t\varphi_n(s)ds+\sigma B(t)}} {\frac{1}{u_0}+a\prod_{0<t_k<t}(1+b_k)\int^t_0e^{(r-\frac{\sigma^2}{2})s-c\int_0^s\varphi_n(\tau)d\tau+\sigma B(s)}ds}.$

Then one has

$dv_1(t)\geq\Big(-\alpha v_1(t)+\alpha\prod\limits_{0<t_k<t}(1+b_k)\tilde{\phi}(t)\Big)dt.$

Arguing as above, we can get

$v_1(t)\geq v_{10}e^{-\alpha t}+\alpha\prod\limits_{0<t_k<t}(1+b_k)\int_0^t\tilde{\phi}(s)e^{-\alpha(t-s)}ds\triangleq\tilde{\varphi}_1(t), t\in[0, \tau_e), a.s.$

Similarly, we have $v_i(t)\geq \tilde{\varphi}_i(t), \ i=2, 3, \ldots, n$, where

$\tilde{\varphi}_i(t)=v_{i0}e^{-\alpha t}+\alpha\int_0^t\tilde{\varphi}_{i-1}(s)e^{-\alpha(t-s)}ds.$

To sum up, we have that

$\tilde{\phi}(t)\leq u(t)\leq\phi(t), \ \ \tilde{\varphi}_i(t)\leq v(t)\leq \varphi_i(t), \ i=1, 2, \ldots, n, \ t\in(0, \infty), \ a.s.$

(11)

This completes the proof of the lemma.

By Lemma 2.2, we have the following result on the global existence of the positive solution to system (5).

Theorem 2.1. There is a unique positive global solution $(x(t), y_1(t), \ldots, y_n(t))$ to system (5) a.s. for any initial value $(x_0, y_{10}, \ldots, y_{n0})\in R_+^{n+1}$.

---

3. Asymptotic behavior of system (5)

In this section, we will first investigate the global attractivity of system (5), then prove the existence of its stationary distribution, and finally, we perform an extinction analysis of the system. We first give the following fundamental assumptions on the impulsive perturbations $b_k$.

• Assumption 3.1. $1+b_k>0$ for all $k\in N$.

• Assumption 3.2. There are two positive constants $\theta_1$ and $\theta_2$ such that for all $t>0$,

$\theta_1\leq \prod\limits_{0<t_k<t}(1+b_k)\leq \theta_2.$

(12)

---

3.1. Global attractivity of system (5)

Before proving the global attractivity of system (5), we first prepare some lemmas.

Lemma 3.1. Let $(u(t), v_1(t), \ldots, v_n(t))$ be a solution of system (9) with any initial value $(u(0), v_1(0), \ldots, v_n(0))\in R^{n+1}_+$. Then there are $q>1$, $K_1(q)$ and $K_2(q)$ such that

$\limsup\limits_{t\rightarrow\infty}E(u^q(t))\leq K_1(q)\ \ \ and \ \ \ \limsup\limits_{t\rightarrow\infty}E(v_i^q(t))\leq K_{2}(q), \ i=1, 2, \ldots, n.$

Proof. Applying Itô's formula to system (9), we compute

$E(e^tu^q(t))=u^q(0)+E\int_0^te^sV(s)ds,$

where

$\begin{align} V(t)&=-aq\prod\limits_{0<t_k<t}(1+b_k)u^{q+1}(t)+\big(rq+\frac{1}{2}q(q-1)\sigma^2+1\big)u^q(t)-cqu^q(t)v_n(t)\notag\\ &\leq-aq\theta_1u^{q+1}(t)+\big(rq+\frac{1}{2}q(q-1)\sigma^2+1\big)u^q(t)\notag\\ &\leq K_1(q).\notag \end{align}$

Therefore,

$E(e^tu^q(t))\leq u^q(0)+ K_1(q)e^t.$

So

$\limsup\limits_{t\rightarrow\infty}E(u^q(t))\leq K_1(q).$

(13)

On the other hand, one can compute that

$dv_1^q(t)=qv_1^{q-1}(t)\Big(-\alpha v_1(t)+\alpha\prod\limits_{0<t_k<t}(1+b_k)u(t)\Big)dt.$

It follows that

$\begin{align} \frac{dE(v_1^q(t))}{dt}&= q\alpha E\Big[v_1^{q-1}(t)\Big(- v_1(t)+\prod\limits_{0<t_k<t}(1+b_k)u(t)\Big)\Big] \notag\\ &\leq -q\alpha E(v_1^q(t))+q\alpha\theta_2 E(u(t)v_1^{q-1}(t)). \notag \end{align}$

Using Hölder inequality, one can then derive that

$\frac{dE(v_1^q(t))}{dt}\leq q\alpha E(v_1^q(t))\Big[-1+\theta_2(Eu^q(t))^{\frac{1}{q}}(E(v_1^q(t)))^{-\frac{1}{q}}\Big].$

(14)

From (13) we know that for any $\varepsilon>0$, there exists a $T>0$ such that when $t>T$,

$E(u^q(t))\leq K_1(q)+\varepsilon.$

It then follows from (14) that

$\frac{dE(v_1^q(t))}{dt}\leq q\alpha E(v_1^q(t))\Big[-1+\theta_2(K_1(q)+\varepsilon)^{\frac{1}{q}}(E(v_1^q(t)))^{-\frac{1}{q}}\Big].$

Using comparison theorem and noting also that $\varepsilon$ is arbitrary, we obtain

$\limsup\limits_{t\rightarrow+\infty} E(v_1^q(t))\leq\theta_2^qK_1(q)\triangleq K_2(q).$

(15)

Similarly, for $i=2, 3, \ldots, n$, we compute that

$\begin{align} \frac{dE(v_i^q(t))}{dt}&= q\alpha [-E(v_i^q(t))+E(v_i^{q-1}(t)v_{i-1}(t))] \notag\\ &\leq q\alpha E(v_i^q(t))[-1+(Ev_{i-1}^q(t))^{\frac{1}{q}}(E(v_i^q(t)))^{-\frac{1}{q}}]. \notag \end{align}$

Then we can deduce that

$\limsup\limits_{t\rightarrow+\infty} E(v_i^q(t))\leq K_2(q).$

This completes the proof of Lemma 3.1.

Definition 3.1. (See [18]) Let $X(t), \bar{X}(t)$ be two arbitrary solutions of Eq. (5) with initial values $X(0)>0$ and $\bar{X}(0)>0$, respectively. If

$\lim\limits_{t\rightarrow+\infty}|X(t)-\bar{X}(t)|=0 \ a.s.,$

then we say Eq. (5) is globally attractive.

Lemma 3.2. (See [13]) Suppose that a stochastic process $X(t)$ on $t\geq0$ satisfies the condition

$E|X(t)-X(s)|^{\alpha_1}\leq c_1|t-s|^{1+\beta}, \ 0\leq s, \ t<\infty$

for some positive constants $\alpha_1, \beta$ and $c_1$. Then there exists a continuous modification $\tilde{X}(t)$ of $X(t)$ which has the property that for every $\gamma\in (0, \beta/\alpha_1)$, there is a positive random variable $h(\omega)$ such that

$\mathscr{P}\left\{\omega:\sup\limits_{0<|t-s|<h(\omega), 0\leq s, t<\infty} \frac{|\tilde{X}(t, \omega)-\tilde{X}(s, \omega)|}{|t-s|^\gamma}\leq\frac{2}{1-2^{-\gamma}}\right\}=1.$

In other words, almost every sample path of $\tilde{X}(t)$ is locally but uniformly Hölder continuous with exponent $\gamma$.

Lemma 3.3. Let $(u(t), v_1(t), \ldots, v_n(t))$ be a solution of system (9) for any positive initial value $(u(0), v_1(0), \ldots, v_n(0))$, then almost every sample path of $(u(t), v_1(t),$ $\ldots, v_n(t))$ is uniformly continuous for $t\geq 0$.

Proof. From (13) and the continuity of $E(u^q(t))$ we have that there is a $K^*_1(q)$ such that for $t\geq 0$,

$E(u^q(t))\leq K_1^*(q).$

Similarly, by (15) and the continuity of $E(v_i^q(t))$, there is a $K^*_2(q)>0$ such that for $t\geq 0$,

$E(v_i^q(t))\leq K^*_2(q).$

The first equation of system (9) is equivalent to the following stochastic integral equation

$u(t)=u_0+\int_0^tu(s)\Big[r-a\prod\limits_{0<t_k<s}(1+b_k)u(s)-c v_n(s)\Big]ds+\int_0^t\sigma u(s)dB(s).$

Then

$\begin{align} &E\Big|u(t)\Big[r-a\prod\limits_{0<t_k<t}(1+b_k)u(t)-c v_n(t)\Big]\Big|^q \notag\\ \leq& 0.5 E|(u(t))|^{2q}+0.5E\Big|r-a\prod\limits_{0<t_k<t}(1+b_k)u(t)-c v_n(t)\Big|^{2q} \notag\\ \leq& 0.5E|u(t)|^{2q}+0.5\times3^{2q-1}\Big[r^{2q}+\Big(a\prod\limits_{0<t_k<t}(1+b_k)\Big)^{2q}E|u(t)|^{2q}+c^{2q}E|v_n(t)|^{2q}\Big] \notag\\ \leq& 0.5 K_1^*(2q)+0.5\times3^{2q-1}\Big[r^{2q}+(a\theta_2)^{2q}K_1^*(2q)+c^{2q}K_2^*(2q)\Big] \notag\\ \triangleq& L_1(q). \notag \end{align}$

By the moment inequality for stochastic integrals, we know that for $0\leq t_1\leq t_2$ and $q>2$,

$E\Big|\int_{t_1}^{t_2}\sigma u(s)dB(s)\Big|^q\leq \sigma^{q}\Big[\frac{q(q-1)}{2}\Big]^{\frac{q}{2}} (t_2-t_1)^{\frac{q}{2}-1}\int_{t_1}^{t_2}E|u(s)|^qds.$

For $0<t_1<t_2<\infty$, $t_2-t_1\leq1$, and $\frac{1}{p}+\frac{1}{q}=1$, we have

$\begin{align} &E(|u(t_2)-u(t_1)|^q) \notag\\ = &E\Big|\int_{t_1}^{t_2}u(s)[r-a\prod\limits_{0<t_k<s}(1+b_k)u(s)-c v_n(s)]ds+\int_{t_1}^{t_2}\sigma u(s)dB(s)\Big|^q \notag\\ \leq &2^{q-1}E\Big|\int_{t_1}^{t_2}u(s)[r-a\prod\limits_{0<t_k<s}(1+b_k)u(s)-c v_n(s)]ds\Big|^q \notag\\ &+2^{q-1}E\Big|\int_{t_1}^{t_2}\sigma u(s)dB(s)\Big|^q \notag\\ \leq&2^{q-1}(t_2-t_1)^{q-1}\int_{t_1}^{t_2}E\Big|u(s)\Big[r-a\prod\limits_{0<t_k<s}(1+b_k)u(s)-c v_n(s)\Big]\Big|^qds \notag\\ & +2^{q-1}\sigma^q\Big[\frac{q(q-1)}{2}\Big]^{\frac{q}{2}} (t_2-t_1)^{\frac{q}{2}-1}\int_{t_1}^{t_2}E|u(s)|^qds \notag\\ \leq &2^{q-1}(t_2-t_1)^{q}L_1(q) +2^{q-1}\sigma^{q}\Big[\frac{q(q-1)}{2}\Big]^{\frac{q}{2}} (t_2-t_1)^{\frac{q}{2}}K_1^*(q) \notag\\ \leq&2^{q-1}(t_2-t_1)^{\frac{q}{2}}\Big[(t_2-t_1)^{\frac{q}{2}}L_1(q) +\sigma^{q}\Big[\frac{q(q-1)}{2}\Big]^{\frac{q}{2}}K_1^*(q)\Big] \notag\\ \leq&(t_2-t_1)^{\frac{q}{2}}L_2(q), \notag \end{align}$

where $L_2(q)=2^{q-1}\Big[L_1(q) +\sigma^{q}\Big[\frac{q(q-1)}{2}\Big]^{\frac{q}{2}}K_1^*(q)\Big]$. Then it follows from Lemma 3.2 that almost every sample path of $u(t)$ is locally but uniformly Hölder-continuous with exponent $\gamma_1$ for every $\gamma_1 \in(0, \frac{q-2}{2q})$.

From the second equation of system (9), we have

$v_1(t)=v_{10}+\int_0^t\Big(-\alpha v_1(s)+\alpha\prod\limits_{0<t_k<s}(1+b_k)u(s)\Big)ds.$

Then for $0<t_1<t_2<\infty$ and $t_2-t_1\leq1$, we compute that

$\begin{align} &E|v_1(t_2)-v_1(t_1)|^q \notag\\ =& E\Big|\int_{t_1}^{t_2}(-\alpha v_1(s)+\alpha\prod\limits_{0<t_k<s}(1+b_k)u(s))ds\Big|^q \notag\\ \leq& (t_2-t_1)^{q-1}\int_{t_1}^{t_2}E\Big|-\alpha v_1(s)+\alpha\prod\limits_{0<t_k<s}(1+b_k)u(s)\Big|^qds \notag\\ \leq & 2^{q-1}(t_2-t_1)^{q-1}\int_{t_1}^{t_2}\Big[\alpha E|v_1(s)|^q+\alpha\prod\limits_{0<t_k<s}(1+b_k)E|u(s)|^q\Big]ds \notag\\ \leq& 2^{q-1}(t_2-t_1)^{q}\Big[\alpha K_2^*(q)+\alpha\prod\limits_{0<t_k<t}(1+b_k)K_1^*(q)\Big] \notag\\ \triangleq&2^{q-1}(t_2-t_1)^{q}L_3(q). \notag \end{align}$

Similarly, for $i=2, \ldots, n$, we can compute that

$\begin{align} E|v_i(t_2)-v_i(t_1)|^q &=E\Big|\int_{t_1}^{t_2}(-\alpha v_i(s)+\alpha v_{i-1}(s))ds\Big|^q \\& \leq 2^{q-1}(t_2-t_1)^{q-1}\int_{t_1}^{t_2}\Big[\alpha E|v_i(s)|^q+\alpha E|v_{i-1}(s)|^q\Big]ds \notag\\ &\leq 2^{q-1}(t_2-t_1)^{q}\alpha K_2^*(q). \notag \end{align}$

In view of Lemma 3.2, almost every sample path of $v_i(t)(i=1, 2, \ldots, n)$ is locally but uniformly Hölder-continuous with exponent $\gamma_2$ for every $\gamma_2 \in(0, \frac{q-1}{q})$.

Therefore, almost every sample path of $(u(t), v_1(t), \ldots, v_n(t))$ is uniformly continuous on $t\geq 0$.

Lemma 3.4. (See [4]) Let $f$ be a non-negative function defined on $R_+=[0, +\infty)$ such that $f$ is integrable on $R_+$, and is uniformly continuous on $R_+$, then

$\lim\limits_{t\rightarrow\infty}f(t)=0.$

Now, we are in a position to prove the global attractivity of system (5). We have the following theorem.

Theorem 3.1. Under Assumptions 3.1 and 3.2, if $c<a$, then system (5) is globally attractive.

Proof. Let $(x(t), y_1(t), \ldots, y_n(t))$ and $(\bar{x}(t), \bar{y}_1(t), \ldots, \bar{y}_n(t))$ be two arbitrary solutions of system (5) with positive initial values, and suppose $x(t)=\prod_{0<t_k<t}(1+b_k)u(t)$, $y_i=v_i(t)$, $\bar{x}(t)=\prod_{0<t_k<t}(1+b_k)\bar{u}(t)$, $\bar{y}_i=\bar{v}_i(t)$. Then for $V_1(t)=|\ln u(t)-\ln \bar{u}(t)|$, the right differential $d^+V_1(t)$ of $V_1(t)$ is

$\begin{align} d^+V_1(t)&=\textrm{sgn}(u(t)-\bar{u}(t))d(\ln u(t)-\ln \bar{u}(t)) \notag\\ &=\textrm{sgn}(u(t)-\bar{u}(t))\Big[-a\prod\limits_{0<t_k<t}(1+b_k)(u(t)-\bar{u}(t))-c(v_n(t)-\bar{v}_n(t))\Big]dt \notag\\ &\leq\Big[-a\prod\limits_{0<t_k<t}(1+b_k)|u(t)-\bar{u}(t)|+c |v_n(t)-\bar{v}_n(t)|\Big]dt. \notag \end{align}$

Define $V_2=\sum_{i=1}^n|v_i-\bar{v}_i|$, by directly calculating the right differential $D^+V_2(t)$ of $V_2(t)$ and then making use of Itô's formula, we have

$\begin{align} d^+V_2(t)=&\sum\limits_{i=1}^n\textrm{sgn}(v_i(t)-\bar{v}_i(t))d(v_i(t)-\bar{v}_i(t)) \notag\\ \leq&\Big[-\alpha|v_1(t)-\bar{v}_1(t)|+\alpha\prod\limits_{0<t_k<t}(1+b_k)|u(t)-\bar{u}(t)|\Big]dt \notag\\ &+\Big[-\alpha\sum\limits_{i=2}^n|v_i(t)-\bar{v}_i(t)|+\alpha \sum\limits_{i=2}^n|v_{i-1}(t)-\bar{v}_{i-1}(t)|\Big]dt\notag\\ =&\Big[-\alpha|v_n(t)-\bar{v}_n(t)|+\alpha\prod\limits_{0<t_k<t}(1+b_k)|u(t)-\bar{u}(t)|\Big]dt \notag \end{align}$

Then for $V(t)=V_1(t)+\lambda V_2(t), \ \lambda\in(c/\alpha, a/\alpha)$, one has

$d^+V(t)\leq\Big[-(a-\lambda\alpha) \prod\limits_{0<t_k<t}(1+b_k)|u(t)-\bar{u}(t)|-(\alpha\lambda-c)|v_n(t)-\bar{v}_n(t)|\Big]dt.$

(16)

Integrating both sides of (16) from $0$ to $t$ yields

$V(t)+\int_0^t\Big[(a-\lambda\alpha) \prod\limits_{0<t_k<s}(1+b_k)|u(s)-\bar{u}(s)|+(\alpha\lambda-c)|v_n(s)-\bar{v}_n(s)|\Big]ds\leq V(0)<\infty,$

which leads to

$|u(t)-\bar{u}(t)|\in L^1[0, \infty)\ \textrm{and} \ |v_i(t)-\bar{v}_i(t)|\in L^1[0, \infty), \ i=1, 2, \ldots, n.$

It follows from Lemmas 3.3 and 3.4 that

$\lim\limits_{t\rightarrow\infty}|u(t)-\bar{u}(t)|=\lim\limits_{t\rightarrow\infty}|v_i(t)-\bar{v}_i(t)|=0.$

Consequently, by Assumption 3.2, we obtain

$\lim\limits_{t\rightarrow\infty}|x(t)-\bar{x}(t)|=\lim\limits_{t\rightarrow\infty}|y_i(t)-\bar{y}_i(t)|=0, \ i=1, 2, \ldots, n.$

---

3.2. Stationary distribution of system (5)

In this subsection, we will prove the existence of a stationary distribution of system (5). Due to the complexity of system (5) and the lake of effective mathematical techniques available, we only consider the special case for Gamma distribution delay kernel with $n=1$, i.e., $f(s)=\alpha e^{-\alpha s}$, which is called weak delay kernel, indicating that the maximum weighted response of the growth rate is due to current population density while past densities have (exponentially) decreasing influence [24].

For convenience, we first study system (9) since systems (5) and (9) are equivalent when we set $x(t)=\prod_{0<t_k<t}(1+b_k)u(t)$ and $y(t)=v(t)$. Furthermore, we assume that

$\lim\limits_{t\rightarrow\infty}\prod\limits_{0<t_k<t}(1+b_k)=\theta$

(17)

holds and let $(\tilde{u}(t), \tilde{v}(t))$ be the solution of

$\left\{\begin{array}{l} d\tilde{u}(t)=\tilde{u}(t)\left(r-a\theta \tilde{u}(t)-c \tilde{v}(t)\right)dt+\sigma \tilde{u}(t)dB(t), \\ d\tilde{v}(t)=\left(-\alpha \tilde{v}(t)+\alpha\theta \tilde{u}(t)\right)dt. \end{array}\right.$

(18)

Notice that system (18) is the limit system of (9), then by the global attractivity of (9), we only need to study the existence of stationary distribution of system (18).

Let $\tilde{u}(t)=e^{\tilde{\xi}(t)}$ and $\tilde{v}(t)=e^{\tilde{\eta}(t)}$, then system (18) becomes

$\left\{\begin{array}{l} d\tilde{\xi}(t)=(r-\frac{\sigma^2}{2}-a\theta e^{\tilde{\xi}(t)}-c e^{\tilde{\eta}(t)})dt+\sigma dB(t), \\ d\tilde{\eta}(t)=(-\alpha+\alpha\theta e^{\tilde{\xi}(t)-\tilde{\eta}(t)})dt. \end{array}\right.$

(19)

We now study the existence of stationary distribution of the equivalent system (19) of (18). Let $X=R^2$, $\Sigma$ be the $\sigma-$algebra of Borel subsets of $X$, and $m$ be the Lebesgue measure on $(X, \Sigma)$. We denote by $\mathscr{P}(t, \tilde{x}, \tilde{y}, A)$ the transition probability function for the diffusion process $(\tilde{\xi}, \tilde{\eta})$, i.e. $\mathscr{P}(t, \tilde{x}, \tilde{y}, A)$ $=\text{Prob}((\tilde{\xi}, \tilde{\eta})\in A)$ and $(\tilde{\xi}, \tilde{\eta})$ is the solution of system (19) with the initial condition $(\tilde{\xi}(0), \tilde{\eta}(0))=(\tilde{x}, \tilde{y})$. If the distribution of $(\tilde{\xi}, \tilde{\eta})$ is absolutely continuous with respect to the Lebesgue measure with the density $U(t, \tilde{x}, \tilde{y})$, then $U(t, \tilde{x}, \tilde{y})$ satisfies the following Fokker-Planck equation:

$\frac{\partial U}{\partial t}=\frac{1}{2}\sigma^2\frac{\partial^2U}{\partial \tilde{x}^2}-\frac{\partial(f_1(\tilde{x}, \tilde{y})U)}{\partial \tilde{x}}-\frac{\partial(f_2(\tilde{x}, \tilde{y})U)}{\partial \tilde{y}},$

(20)

where

$f_1(\tilde{x}, \tilde{y})=r-\frac{\sigma^2}{2}-a\theta e^{\tilde{x}}-c e^{\tilde{y}}, f_2(\tilde{x}, \tilde{y})=-\alpha+\alpha\theta e^{\tilde{x}-\tilde{y}}.$

(21)

Then for the asymptotical stability of system (19), we have the following theorem.

Theorem 3.2. Let $(\tilde{\xi}, \tilde{\eta})$ be the solution of system (19). Then for every $t (>0)$ the distribution of $(\tilde{\xi}, \tilde{\eta})$ has the density $U(t, \tilde{x}, \tilde{y})$ satisfying (20). Furthermore, if $r>\frac{1}{2}\sigma^2$, then there exists a unique density $U^*(\tilde{x}, \tilde{y})$ such that

$\lim\limits_{t\rightarrow\infty}\iint_{R^2}|U(t, \tilde{x}, \tilde{y})-U^*(\tilde{x}, \tilde{y})|d\tilde{x}d\tilde{y}=0$

(22)

and

$supp U^*=R^2.$

Remark 1. By the support of a measurable function $f$ we simply mean the set

$\textrm{ supp} f=\{(\tilde{x}, \tilde{y})\in R^2:f(\tilde{x}, \tilde{y})\neq0\}.$

Besides, note that the Fokker-Planck equation corresponding to system (19) is of a degenerate type, thus the asymptotic stability of the system can't follow directly from the known results from Hasminskii [9]. We will show it by using the theory of integral Markov semigroups (see Appendix A). Now we introduce an integral Markov semigroup connected with system (19).

Denote by $k(t, \tilde{x}, \tilde{y};\tilde{x}_0, \tilde{y}_0)$ the density of $\mathscr{P}(t, \tilde{x}_0, \tilde{y}_0, \cdot)$. Then

$\mathcal{P}(t)f(\tilde{x}, \tilde{y})=\iint_{R^2} k(t, \tilde{x}, \tilde{y};u, v)f(u, v)dudv$

and consequently $\{\mathcal{P}(t)\}_{t\geq 0}$ is an integral Markov semigroup. Thus the asymptotic stability of the semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ implies that all the densities of the process $(\tilde{\xi}(t), \tilde{\eta}(t))$ convergence to an invariant density in $\mathcal{L}^1$. Therefore, for Theorem 3.2, we only need to show the asymptotic stability of the semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$. The outline of our proof is as follows:

(ⅰ) First, using the Hörmander condition [21], we show that the transition function of the process $(\tilde{\xi}(t), \tilde{\eta}(t))$ is absolutely continuous;

(ⅱ) Then according to support theorems [5,1], we find a set $E$ on which the density of the transition function is positive;

(ⅲ) Next we show that the set $E$ is an attractor and the semigroup satisfies the Foguel alternative [25,22,23];

(ⅳ) Finally, we exclude sweeping by showing that there exists a Khasminski$\breve{\text{i}}$ function.

In the following, we give the proof of Theorem 3.2 through four lemmas in succession, which correspond respectively to (ⅰ)-(ⅳ) above.

Lemma 3.5. The transition probability function $\mathscr{P}(t, \tilde{x}_0, \tilde{y}_0, A)$ has a continuous density $k(t, \tilde{x}, \tilde{y};$ $\tilde{x}_0, \tilde{y}_0)$.

Proof. The proof is based on the Hörmander theorem for the existence of smooth densities of the transition probability function for degenerate diffusion processes. If $a(x)$ and $b(x)$ are vector fields on $R^d$, then the Lie bracket $[a, b]$ is a vector field given by

$[a, b]_j(X)=\sum\limits_{k=1}^d\Big(a_k\frac{\partial b_j}{\partial x_k}(X)-b_k\frac{\partial a_j}{\partial x_k}(X)\Big)^T, \ \ \ j=1, 2, \ldots, d.$

Let

$a(\tilde{x}, \tilde{y})=(r-\frac{\sigma^2}{2}-a\theta e^{\tilde{x}}-c e^{\tilde{y}}, -\alpha+\alpha\theta e^{\tilde{x}-\tilde{y}})^T$

and

$b(\tilde{x}, \tilde{y})=(\sigma, 0)^T.$

Then

$[a, b](\tilde{x}, \tilde{y})=\big(\sigma a\theta e^{\tilde{x}}, -\sigma\alpha\theta e^{\tilde{x}-\tilde{y}}\big)^T.$

It follows that

$\left|\begin{array}{cc} \sigma & 0\\ \sigma a\theta e^{\tilde{x}} & -\sigma\theta e^{\tilde{x}-\tilde{y}} \end{array}\right| =-\sigma^2\alpha\theta e^{\tilde{x}-\tilde{y}}.$

Thus for every $(\tilde{x}, \tilde{y})\in R^2$, vectors $b(\tilde{x}, \tilde{y})$ and $[a, b](\tilde{x}, \tilde{y})$ span the space $R^2$. This implies that the transition probability function $\mathscr{P}(t, \tilde{x}_0, \tilde{y}_0, A)$ has a density $k(t, \tilde{x}, \tilde{y};\tilde{x}_0, \tilde{y}_0)$ and $k\in C^\infty \big((0, \infty)\times R^2\times R^2\big)$. This completes the proof of lemma 3.5.

Lemma 3.6. Let $E=R^2$. Then for each $(\tilde{x}_0, \tilde{y}_0)\in E$ and $(\tilde{x}, \tilde{y})\in E$, there exists $T>0$ such that $k(T, \tilde{x}, \tilde{y};\tilde{x}_0, \tilde{y}_0)>0$.

Proof. We now use support theorems to check that the kernel $k$ is positive. Fix a point $(\tilde{x}_0, \tilde{y}_0)\in R^2$ and a continuous function $\phi\in \mathcal{L}^2([0, T];R)$. Consider the following system

$\tilde{x}_\phi(t)=\tilde{x}_0+\int_0^t [f_1(\tilde{x}_\phi(s), \tilde{y}_\phi(s))+\sigma \phi]ds,$

(23)

$\tilde{y}_\phi(t)=\tilde{y}_0+\int_0^t f_2(\tilde{x}_\phi(s), \tilde{y}_\phi(s))ds,$

(24)

where $f_1(\tilde{x}_\phi, \tilde{y}_\phi)=r-\frac{\sigma^2}{2}-a\theta e^{\tilde{x}_\phi}-c e^{\tilde{y}_\phi}$ and $f_2(\tilde{x}_\phi, \tilde{y}_\phi)=-\alpha+\alpha\theta e^{\tilde{x}_\phi-\tilde{y}_\phi}$.

Let $D_{\tilde{x}_0, \tilde{y}_0;\phi}$ be the Frechét derivative of the function $h\mapsto X_{\phi+h}(T)$ with $X_{\phi+h}=[x_{\phi+h}, y_{\phi+h}]^T$. If for some $\phi$ the derivative $D_{\tilde{x}_0, \tilde{y}_0;\phi}$ has rank $2$, then we get $k(T, \tilde{x}, \tilde{y};\tilde{x}_0, \tilde{y}_0)>0$ for $x=\tilde{x}_\phi(T)$ and $y=\tilde{y}_\phi(T)$. The derivative $D_{x_0, \tilde{y}_0;\phi}$ can be found by means of the perturbation method for ordinary differential equations. In other words, let $\Gamma(t)={\bf f}'(\tilde{x}_\phi(t), \tilde{y}_\phi(t))$, where ${\bf f'}$ is the Jacobians of ${\bf f}=[f_1(\tilde{x}, \tilde{y}), f_2(\tilde{x}, \tilde{y})]^T$. Let $Q(t, t_0)$ be a matrix function such that $Q(t_0, t_0)=I$, $\frac{\partial Q(t, t_0)}{\partial t}=\Gamma (t)Q(t, t_0)$ for $T\geq t\geq t_0\geq 0$ and ${\bf v}=(\sigma, 0)^T$. Then

$D_{\tilde{x}_0, \tilde{y}_0;\phi}h=\int_0^t Q(T, s){\bf v}h(s)ds.$

We first check that the rank of $D_{\tilde{x}_0, \tilde{y}_0;\phi}$ is $2$. Let $\varepsilon \in(0, T)$ and $h={\bf 1}_{[T-\varepsilon, T]}(t)$, where $t\in [0.T]$ and ${\bf 1}_{[T-\varepsilon, T]}$ is the characteristic function of interval $[T-\varepsilon, T]$. Since $Q(T, s)=I+\Gamma(T)(T-s)+o(T-s)$, we obtain

$D_{\tilde{x}_0, \tilde{y}_0;\phi}h=\varepsilon {\bf v}+\frac{1}{2}\varepsilon^2 \Gamma(T){\bf v}+o(\varepsilon^2).$

Then, we have

$\begin{array}{lll} \Gamma(T){\bf v} = \left [ \begin{array}{cc} -a\theta e^{\tilde{x}} & -c e^{\tilde{y}}\\ \alpha\theta e^{\tilde{x}-\tilde{y}} & -\alpha\theta e^{\tilde{x}-\tilde{y}} \end{array} \right] \left [ \begin{array}{l} \sigma\\ 0 \end{array} \right] = \sigma e^{\tilde{x}}\left[ \begin{array}{l} -a\theta\\ \alpha\theta e^{-\tilde{y}} \end{array} \right]. \end{array}$

Hence, vectors ${\bf v}$ and $\Gamma(T){\bf v}$ are linearly independent. Thus $D_{\tilde{x}_0, \tilde{y}_0;\phi}$ has rank $2$.

Next, we prove that for any two points $(\tilde{x}_0, \tilde{y}_0)\in E$ and $(x_T, y_T)\in E$, there exists a control function $\phi$ and $T>0$ such that $\tilde{x}_\phi(0)=\tilde{x}_0, \ \tilde{y}_\phi(0)=\tilde{y}_0, \ \tilde{x}_\phi(T)=x_T$ and $\tilde{y}_\phi(T)=y_T$. Taking derivatives of systems (23) and (24) yield

$\tilde{x}'_\phi(t)=f_1(\tilde{x}_\phi(t), \tilde{y}_\phi(t))+\sigma\phi,$

(25)

$\tilde{y}'_\phi(t)=f_2(\tilde{x}_\phi(t), \tilde{y}_\phi(t)).$

(26)

We construct the function $\phi$ in the following way. First, we find a positive constant $T$ and a differential function $\tilde{x}_\phi:[0, T]\rightarrow \mathbb{R_+}$ such that $\tilde{y}_\phi(0)=\tilde{y}_0$, $\tilde{y}_\phi(T)=\tilde{y}_T$,

$\tilde{y}'_\phi(0)=-\alpha+\alpha\theta e^{\tilde{x}_0-\tilde{y}_0}, \ \ \tilde{y}'_\phi(T)=-\alpha+\alpha\theta e^{\tilde{x}_T-\tilde{y}_T}$

(27)

and

$\tilde{y}'_\phi(t)+\alpha=\alpha\theta e^{\tilde{x}_t-\tilde{y}_t}>0 \ \ \text{for} \ \ t\in[0, T].$

(28)

Denote constants $\tilde{y}'_\phi(0)$ and $\tilde{y}'_\phi(T)$ by $a_0$ and $a_T$. We split the construction of the function $\tilde{y}_\phi$ on three intervals $[0, \varepsilon]$, $[\varepsilon, T-\varepsilon]$ and $[T-\varepsilon, T]$, where $0<\varepsilon<T/2$. By (27), we can construct a $C^2$ function $\tilde{y}_\phi:[0, \varepsilon]\rightarrow R$ such that $\tilde{y}_\phi(0)=\tilde{y}_0$, $\tilde{y}'_\phi(0)=a_0$, $\tilde{y}'_\phi(\varepsilon)=0$ and $\tilde{y}_\phi$ satisfies inequality (28) for $t\in[0, \varepsilon]$. Similarly, we construct a function $\tilde{y}_\phi:[T-\varepsilon, T]\rightarrow R$ such that $\tilde{y}_\phi(T)=\tilde{y}_T$, $\tilde{y}'_\phi(T)=a_T$, $\tilde{y}'_\phi(T-\varepsilon)=0$ and $\tilde{y}_\phi$ satisfies inequality (28) for $t\in[T-\varepsilon, T]$. Let $T$ be sufficiently large, we can extend the function $\tilde{y}_\phi: [0, \varepsilon]\cup [T-\varepsilon, T]\rightarrow R$ to a $C^2$ function $\tilde{x}_\phi$ defined on the whole interval $[0, T]$ such that $\tilde{y}_\phi$ satisfies inequality (28). From (28), we can find a $C^1$ function $\tilde{x}_\phi$ which satisfies $(26)$ and finally we can determine a continuous function $\phi$ from (25). This completes the proof of Lemma 3.6.

Notice that for every density $f$, one has

$\lim\limits_{t\rightarrow\infty}\iint_{R^2} \mathcal{P}(t)f(\tilde{x}, \tilde{y})d\tilde{x}d\tilde{y}=1.$

Lemma 3.7. The semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is asymptotically stable or is sweeping with respect to compact sets.

Proof. By Lemma 3.5, the semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is an integral Markov semigroup with a continuous kernel $k(t, \tilde{x}, \tilde{y})$ for $t>0$. Let $E=R^2$, then it is sufficient to investigate the restriction of the semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ to the space $\mathcal{L}^1(\bar{E})$, where $\bar{E}$ denotes the closure of the set $E$. In view of Lemma 3.6, for every $f\in \mathcal{D}$, we have

$\int_0^\infty \mathcal{P}(t)fdt>0 \ \ a.e.\ \ \text{on}\ E,$

where $\mathcal{D}$ is defined in Appendix A. Therefore, according to Lemma A.1 in Appendix A, it follows that the semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is asymptotically stable or is sweeping with respect to compact sets. This completes the proof of Lemma 3.7.

Lemma 3.8. The semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is asymptotically stable.

Proof. In order to exclude sweeping, we now construct a nonnegative $C^2-$function $V$ and a closed set $\Gamma\in \Sigma$ such that

$\sup\limits_{\tilde{x}, \tilde{y}\notin \Gamma}\mathscr{A}^*V(\tilde{x}, \tilde{y})<0,$

where $\mathscr{A}^*$ is the adjoint operator of the infinitesimal generator $\mathscr{A}$ of the semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$, which is of the form

$\mathscr{A}^*V=\frac{1}{2}\sigma^2 \frac{\partial^2V}{\partial \tilde{x}^2}+f_1(\tilde{x}, \tilde{y})\frac{\partial V}{\partial \tilde{x}}+f_2(\tilde{x}, \tilde{y})\frac{\partial V}{\partial \tilde{y }},$

(29)

where $f_i(\tilde{x}, \tilde{y}), \ i=1, 2$ are defined in (21). Define

$V(\tilde{x}, \tilde{y})=\theta(e^{\tilde{x}}-e^{\tilde{x}^*}-e^{\tilde{x}^*}(\tilde{x}-\tilde{x}^*)) +\frac{1}{2}c(e^{\tilde{y}}-e^{\tilde{y}^*})^2,$

where

$\tilde{x}^*=\ln\frac{r}{a\theta+c \theta}, \ \tilde{y}^*=\ln \frac{ r}{a+c}.$

By (29), one has

$\begin{align} &\mathscr{A}^*V(\tilde{x}, \tilde{y})\notag \\ =&\theta(e^{\tilde{x}}-e^{\tilde{y}^*})(r-\frac{\sigma^2}{2}-a\theta e^{\tilde{x}}-c e^{\tilde{y}})+\frac{1}{2}\sigma^2\theta e^{\tilde{x}}+c(e^{\tilde{y}}-e^{\tilde{y}^*})(-\alpha e^{\tilde{y}}+\alpha\theta e^{\tilde{x}}) \notag\\ =&-a\theta^2(e^{\tilde{x}}-e^{\tilde{x}^*})^2-\alpha c(e^{\tilde{y}}-e^{\tilde{y}^*})^2+\frac{1}{2}\sigma^2\theta e^{\tilde{x}^*} \notag \end{align}$

Thus, the ellipsoid

$a\theta^2(e^{\tilde{x}}-e^{\tilde{x}^*})^2+\alpha c(e^{\tilde{y}}-e^{\tilde{y}^*})^2=\frac{1}{2}\sigma^2\theta e^{\tilde{x}^*}$

lies entirely in $R^2$. Thus, there exist a closed set $\Gamma\in \Sigma$ which contains this ellipsoid and $c>0$ such that

$\sup\limits_{\tilde{x}, \tilde{y}\notin \Gamma}\mathscr{A}^*V(\tilde{x}, \tilde{y}) \leq -c<0.$

The function $V$ is called a Khasminski$\breve{\text{i}}$ function. By means of standard arguments similar to those in [22], one can check that the semigroup is not sweeping from the set $\Gamma$ due to the existence of a Khasminski$\breve{\text{i}}$ function. Therefore, according to Lemma A.1 in Appendix A, we conclude that the semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is asymptotically stable.

In view of Theorem 3.2, for system (18) there exists a unique positive invariant density with the support set $R^2_+$ since it is equivalent to system (19). Notice that system (18) is the limit system of (9) as $\lim_{t\rightarrow\infty}\prod_{0<t_k<t}(1+b_k)=\theta$, then by the global attractivity of (9), we conclude that there exists a unique positive invariant density for system (9). Further notice that system (5) is equivalent to system (9). Therefore, we have the following result for system (5).

Theorem 3.3. Assume that condition (17) holds. If $r>\frac{1}{2}\sigma^2$ and $c<a$, then for system (5) there exists a unique positive invariant density with the support set $R^2_+$.

---

3.3. Extinction of system (5)

In this subsection, we show that large noise can lead to the extinction of system (5).

Theorem 3.4. Let $(x(t), y(t))$ be a solution of system (5). Then

$\lim\limits_{t\rightarrow\infty}\frac{\ln x(t)}{t}\leq r-\frac{1}{2}\sigma^2.$

In particular, if $r<\frac{1}{2}\sigma^2$, then $\lim_{t\rightarrow \infty}x(t)=0$ a.s.

Proof. Applying Itô's formula to system (9), one has

$d\ln u(t)=\Big(r-\frac{1}{2}\sigma^2-a\prod\limits_{0<t_k<t}(1+b_k)u(t)-c v(t)\Big)dt+\sigma dB(t).$

Integrating both sides from $0$ to $t$, we have

$\ln u(t)\leq \ln u(0)+(r-\frac{1}{2}\sigma^2)t+M_1(t),$

(30)

where $M_1(t)=\int_0^t\sigma dB(s)$. Note that $M_1(t)$ is a local martingale, whose quadratic variation is $\langle M_1(t), M_1(t)\rangle=\sigma^2t.$ Making use of the strong law of large numbers for local martingales leads to

$\lim\limits_{t\rightarrow\infty}M_1(t)/t=0, ~ a.s.$

From (30), we know

$\ln x(t)=\sum\limits_{0<t_k<t}\ln(1+b_k)+\ln u(t) \leq \ln u(0) +\ln \theta+(r-\frac{1}{2}\sigma^2)t+M_1(t).$

(31)

It then follows from (31) that

$\lim\limits_{t\rightarrow\infty}\frac{\ln x(t)}{t}\leq r-\frac{1}{2}\sigma^2.$

Remark 2. From Theorems 3.3 and 3.4, one can see that if the intensity of the noise is small, there exists an invariant and asymptotically stable density of the system, while large noise will make the population extinct eventually. Noting further that if $c<a$ and condition (17) hold, then $r-\frac{1}{2}\sigma^2$ can be considered as a threshold determining the asymptotical stability and the extinction of system.

---

4. Simulations and discussions

Impulsive and uncertain variability together with time delay are always present in a natural system, which should be accounted for in its mathematical model. The research performed in this paper is an attempt in this direction, using an impulsive stochastic model with delay. More specifically, the impulse is introduced at fixed moments, the stochastic perturbation is of white noise type and is assumed to be proportional to the population density, and the delay takes the distributed type with a weak delay kernel. To perform a detailed analysis on the dynamics of model (4), we first transform it to an equivalent stochastic system (5) with impulsive effects using the linear chain trick. Then based on Lemma 2.1, it can be further reduced to the problem of a nonlinear stochastic differential system without impulses, i.e., system (9).

For system (5), we first carry out the analysis of its global attractivity. Theorem 3.1 shows that the system is globally attractive provided that $c<a$, Assumptions 3.1 and 3.2 hold. Note that Assumptions 3.1 and 3.2 are only relative to the impulse, so we see that the noise and the delay do not affect the global attractivity of the system. Then we study the existence of the stationary distribution of system (5). It is difficult to directly study its distribution, so we turn to study its equivalent system (9). Notice that (9) is non-autonomous and it has been showed to be globally attractive in Theorem 3.1. So we only need to study the limit system (18) of (9). Due to the Fokker-Planck equation corresponding to system (18) is of a degenerate type, thus the existence of the stationary distribution can't follow directly from the known result from Hasminskii [9]. We show it based on the theory of integral Markov semigroups. Theorem 3.3 shows that for system (5) there exists a unique invariant density with the support $R^2_+$ provided that $r>\frac{1}{2}\sigma^2$, $c<a$ and condition (17) hold. Moreover, it is shown in Theorem 3.4 that if $r<\frac{1}{2}\sigma^2$, the population will be extinct eventually. Thus, under $c<a$ and condition (17), $r-\frac{1}{2}\sigma^2$ can be considered as a threshold determining the asymptotical stability and the extinction of system (5).

To illustrate the results obtained above, some numerical simulations are carried out by using Milstein scheme [10]. Consider the discretization of model (5) for $t=0, \Delta t, 2\Delta t, \ldots, n \Delta t$:

$\left\{\begin{array}{l}\left.\begin{array}{l} x_{i+1}=x_i+x_i(r-ax_i-c y_i)\Delta t+\sigma x_i\sqrt{\Delta t}\xi_i, \\ y_{i+1}=y_i+(-\alpha y_i+\alpha x_i)\Delta t, \end{array}\right\} ~t\neq t_k, \\ x(t_k^+)-x(t_k)=b_kx(t_k), ~t=t_k, k\in N. \end{array}\right.$

where time increment $\Delta t>0$ and $\xi_i$ are $N(0, 1)-$distributed independent random variables which can be generated numerically by pseudorandom number generators.

First, we demonstrate system (5) is asymptotical stability. For this purpose, let $r=0.8, \ \alpha=0.6, \ a=0.4, \ c=0.33, \ t_k=k \ (k\in N)$, initial value $(x(0), y(0))=(1, 1)$, $b_k$ vary in the range $[-0.3/k, 0.3/k]$ and $\sigma$ vary in the range $[0, 0.03]$. Simple computations show that $\sigma^2<2r=1.6$ and $\lim_{t\rightarrow\infty}\prod_{0<t_k<t}(1+b_k)$ exists.

By Theorem 3.3, system (5) is asymptotical stability. Our simulation supports this conclusion as shown in Fig. 1, where we show the effect from different noise intensities $\sigma$ and different impulse intensities $b_k$, respectively.

It is seen in Fig. 1 (a) that the steady state of the system is point $E^*(x^*, y^*)=(1.0909, 1.8182)$ in the absence of impulsive and random perturbations (i.e. $\sigma=b_k=0$). Increasing the value of $\sigma$ and keeping $b_k$ the same, we can see that the steady state of the system can no longer be represented by a single point; instead, it is represented by a region around $E^*$. Shown in Fig. 1 (a)-(c), the larger the noise intensity is, the more diffusive of the system state is. When $b_k$ changes and $\sigma$ stays $0$, we can see in Fig. 1 (d)-(e) that there is an attractor for the system, but $E^*$ is not included in the attraction, which is above the attraction in the case of $b_k<0$ and is below the attraction in the case of $b_k>0$. When we choose $\sigma=0.03$ and $b_k=0.3/k$, the system is still stable, please see Fig. 1 (f).

Fig. 1 shows the trajectories of system (5) under different values of parameters. However, we should point out that the trajectory in each subgraph is drawn for a single sample, which is stochastic for each sample under the same parameters, that is the outcome for a single trajectory is not predictable. But, the probability distribution of all possible outcomes can be determined. Our simulation supports this conclusion as shown in Fig. 2, where $\sigma=0.03, \ b_k=0.3/k$ and the densities are drawn based on $10000$ sample pathes, computed with differential initial value and different iterative times, respectively. We can see from the figure the distribution is stationary.

Next, we show system (5) is extinct. To this end, we set initial value $(x(0), y(0))=(1, 1)$, $r=0.3, \ a=0.4, \ c=0.33, \ \alpha=0.6, \ \ b_k=0, \ \ t_k=5k$ and $\sigma=0.8$. Since $r<\frac{1}{2}\sigma^2=0.32$, by Theorem 3.4, the population is extinct eventually. Our simulation supports this conclusion as shown in Figs. 3 (a) and (b). Decreasing the value of $\sigma$ from $0.8$ to $0.1$, we obtain $r>\frac{1}{2}\sigma^2=0.005$. Then the population is persistent and system (5) is asymptotically stable, please see Figs. 3 (c) and (d). Thus, as we mentioned before, under condition (17), $r-\frac{1}{2}\sigma^2$ can be considered as a threshold determining the asymptotical stability and the extinction of system (5).

To sum up, this paper presents an investigation on the dynamics of a impulsive stochastic system with delay. Our findings are useful for better understanding of the effects of impulses, stochastic perturbations and delay on the dynamics of a system. We should point out there are still some other interesting topics meriting further investigation, for example, the long term behavior of multi-population system with impulsive and stochastic perturbations. We leave these for future considerations.

---

Appendix A.

Let the triple $(X, \Sigma, m)$ be a $\sigma-$finite measure space. Denote by $\mathcal{D}$ be the subset of the space $\mathcal{L}^1$ which contains all densities, i.e.

$\mathcal{D}=\{f\in \mathcal{L}^1: f\geq 0, \|f\|=1\}.$

A linear mapping $\mathcal{P}: \mathcal{L}^1\rightarrow\mathcal{L}^1$ is called a Markov operator if $\mathcal{P}(\mathcal{D})\subset \mathcal{D}$.

The Markov operator $\mathcal{P}$ is called an integral or kernel operator if there exists a measurable function $k:X\times X\rightarrow [0, \infty)$ such that

$\int_X k(\tilde{x}, \tilde{y})m(d\tilde{x})=1$

(A.1)

for all $\tilde{y}\in X$ and

$\mathcal{P}f(x)=\int_X k(\tilde{x}, \tilde{y})f(\tilde{y})m(d\tilde{y})$

for every density $f$.

A family $\{\mathcal{P}(t)\}_{t\geq 0}$ of Markov operators which satisfies conditions:

(ⅰ) $\mathcal{P}(0)=Id$,

(ⅱ) $\mathcal{P}(t+s)=\mathcal{P}(t)\mathcal{P}(s)$ for $s, t\geq 0$,

(ⅲ) for each $f\in \mathcal{L}^1$ the function $t\mapsto \mathcal{P}(t) f$ is continuous with respect to the $\mathcal{L}^1$ norm is called a Markov semigroup. A Markov semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is called integral, if for each $t>0$, the operator $\mathcal{P}(t)$ is an integral Markov operator.

A density $f_*$ is called invariant if $\mathcal{P}(t)f_*=f_*$ for each $t>0$. The Markov semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is called asymptotically stable if there is an invariant density $f_*$ such that

$\lim\limits_{t\rightarrow\infty}\|\mathcal{P}(t)f-f_*\|=0 \ \ \text{for} \ \ f\in \mathcal{D}.$

A Markov semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is called sweeping with respect to a set $A\in \Sigma$ if for every $f\in \mathcal{D}$

$\lim\limits_{t\rightarrow \infty}\int_A \mathcal{P}(t)f(\tilde{x})m(d\tilde{x})=0.$

We need some result concerning asymptotic stability and sweeping which can be found in [7].

Lemma A.1. Let $X$ be a metric space and $\Sigma$ be the $\sigma-$algebra of Borel sets. Let $\{\mathcal{P}(t)\}_{t\geq 0}$ be an integral Markov semigroup with a continuous kernel $k(t, \tilde{x}, \tilde{y})$ for $t>0$, which satisfies (A.1) for all $\tilde{y}\in X$. We assume that for every $f\in \mathcal{D}$, we have

$\int_0^\infty \mathcal{P}(t)fdt>0 \ \ a.e.$

Then this semigroup is asymptotically stable or is sweeping with respect to compact sets.

The property that a Markov semigroup $\{\mathcal{P}(t)\}_{t\geq 0}$ is asymptotically stable or sweeping for a sufficiently large family of sets (e.g. for all compact sets) is called the Foguel alternative.

---

Acknowledgments

We would like to thank anonymous reviewers for helpful suggestions which greatly improved this paper.

---