Global Asymptotic Stability of Stochastic Nonautonomous Lotka-Volterra Models with Infinite Delay

Abstract

A kind of general stochastic nonautonomous Lotka-Volterra models with infinite delay is investigated in this paper. By constructing several suitable Lyapunov functions, the existence and uniqueness of global positive solution and global asymptotic stability are obtained. Further, the solution asymptotically follows a normal distribution by means of linearizing stochastic differential equation. Moment estimations in time average are derived to improve the approximation distribution. Finally, numerical simulations are given to illustrate our conclusions.

1. Introduction

The impact of random factors cannot be neglected in the real world. Different kinds of random perturbations of stochastic models have been investigated in many pieces of literature. Bahar and Mao [1] discussed a stochastic delay Lotka-Volterra model, and they showed that environmental noise would suppress a potential population explosion and also made the solutions to be stochastically ultimately bounded. Almost at the same time, Mao [2] revealed that different types of environmental noise had different effects on delay population models. Meanwhile, Jiang and Shi [3] considered a randomized nonautonomous Logistic equation and represented the unique continuous global positive solution and positive 𝑇-periodic solution.

Recently, stochastic models with delay are paid more attention by many researchers. Shen et al. [4] studied stochastic Lotka-Volterra competitive models with variable delay, and they obtained the unique global positive solution, stochastically ultimate boundedness, and moment average in time of the solutions. Liu and Wang [5] considered stochastic Lotka-Volterra models with infinite delay; they replaced the intrinsic growth rate by a random perturbation term which was dependent on the difference between the population size and the equilibrium state in their model, and then the sufficient criteria for global asymptotic stability of the solution were established. Huang et al. [6] extended the conclusion of Liu and Wang [5] to a general case. Based on a general phase space 𝐶_𝑟∶=𝜑∈𝐶(−∞,0];𝑅^𝑛₊₊∶𝜑_𝐶𝑟=sup_{−∞<𝑠≤0}𝑒^𝑟𝑠𝜑(𝑠())<∞,𝑟>0,H1

Xu et al. [7] and Xu [8] investigated an autonomous stochastic Lotka-Volterra model with infinite delay 𝑑𝑥(𝑡)=diag𝑥₁(𝑡),…,𝑥_𝑛(𝑡)×𝑏+𝐴𝑥(𝑡)+𝐵⁰_−∞𝑥()(𝜃)(𝑡)(𝑡;𝑡+𝜃𝑑𝜇𝑑𝑡+𝜎𝑥𝑑𝐵)(1) they established the asymptotic pathwise properties of the solution to model (1), where 𝑥=(𝑥₁,…,𝑥_𝑛)^𝑇, 𝑏=(𝑏₁,…,𝑏_𝑛)^𝑇, 𝐴=(𝑎_𝑖𝑗)_𝑛×𝑛, 𝐵=(𝑏_𝑖𝑗)_𝑛×𝑛, and 𝜎=(𝜎_𝑖𝑗)_𝑛×𝑛; 𝐶_𝑟 was a Banach space (see [7, 9, 10]); 𝜇 was probability measure defined on (−∞,0] such that𝜇_𝑟∶=⁰_−∞𝑒^−2𝑟𝜃𝑑𝜇(𝜃)<∞,𝑟>0.(H2)

In particular, when 𝜇(𝜃)=𝑒^𝑘𝑟𝜃(𝑘>2) for 𝜃≤0, assumption (H2) was satisfied.

Motivated by the works mentioned previously, we always assume that the intensity of white noise is dependent on the difference between the population size and the equilibrium state and also is dependent on time 𝑡 in this paper. Now, we consider a more general stochastic nonautonomous Lotka-Volterra model with infinite delay 𝑑𝑥_𝑖(𝑡)=𝑥_𝑖(𝑡)ℎ_𝑖(𝑡)+^𝑛_𝑗=1𝑎_𝑖𝑗(𝑡)𝑥_𝑗(𝑡)+^𝑛_𝑗=1𝑏_𝑖𝑗(𝑡)𝑥_𝑗((𝑡+𝑡−𝛿))_𝑛𝑗=1𝑐_𝑖𝑗(𝑡)⁰_−∞𝑥_𝑗(𝑡+𝜃)𝑑𝜇(𝜃)𝑑𝑡+𝑥_𝑖(𝑡)𝜎_𝑖(𝑡)𝑥_𝑖(𝑡)−𝑥^∗_𝑖+^𝑛_𝑗=1𝛼_𝑖𝑗(𝑡)𝑥_𝑗(𝑡)𝑥_𝑖(𝑡)−𝑥^∗_𝑖+^𝑛_𝑗=1𝛽_𝑖𝑗(𝑡)𝑥_𝑗(𝑡−𝛿(𝑡))𝑥_𝑖(𝑡)−𝑥^∗_𝑖+^𝑛_𝑗=1𝛾_𝑖𝑗(𝑡)⁰_−∞𝑥_𝑗()(𝜃)×𝑥𝑡+𝜃𝑑𝜇_𝑖(𝑡)−𝑥^∗_𝑖𝑑𝐵_𝑖(𝑡),𝑖=1,2,…,𝑛,(2) where ℎ_𝑖(𝑡)>0, 𝑎_𝑖𝑗(𝑡) (𝑎_𝑖𝑖(𝑡)<0 if 𝑖=𝑗), 𝑏_𝑖𝑗(𝑡) and 𝑐_𝑖𝑗(𝑡) are parameter functions; noise intensities 𝜎_𝑖(𝑡)>0, 𝛼_𝑖𝑗(𝑡)≥0, 𝛽_𝑖𝑗(𝑡)≥0 and 𝛾_𝑖𝑗(𝑡)≥0 are continuous bounded function on [0,+∞); variable delay function 𝛿∶[0,∞)→[0,𝜏] and 𝛿^′(𝑡)≤0; 𝑥^∗=(𝑥^∗₁,𝑥^∗₂,…,𝑥^∗_𝑛)^𝑇 denotes an equilibrium state with respect to the deterministic part of model (2) in 𝑅^𝑛₊₊ (we always assume that such 𝑥^∗ exists in this paper). For simplicity, model (2) can be rewritten in the following form: 𝑑𝑥_𝑖(𝑡)=𝑥_𝑖(𝑡)𝑓_𝑖𝑑𝑡+𝑥_𝑖(𝑡)𝑔_𝑖𝑑𝐵_𝑖(𝑡),𝑖=1,2,…,𝑛.(3) Next, we will investigate the global positive solution and its asymptotic properties of model (2). Further, the approximation distribution of solutions to model (2) is explored. Our results will extend some classical deterministic results into the stochastic cases.

2. Preliminary

Let (Ω,ℱ,{ℱ_𝑡}_𝑡≥0,𝑃), throughout this paper unless otherwise specified, be a complete probability space with a filtration {ℱ_𝑡}_𝑡≥0 satisfying the usual conditions, and 𝐵_𝑖(𝑡) (𝑖=1,2,…,𝑛) denote the independent 1-dimensional standard Brownian motion defined on the complete probability space. We denote the nonnegative cone and positive cone by 𝑅^𝑛₊, 𝑅^𝑛₊₊, respectively; that is, 𝑅^𝑛₊={𝑥∈𝑅^𝑛∶𝑥_𝑖≥0,𝑖=1,2,…,𝑛}, 𝑅^𝑛₊₊={𝑥∈𝑅^𝑛∶𝑥_𝑖>0,𝑖=1,2,…,𝑛}. If 𝑥∈𝑅^𝑛, its norm is denoted by |𝑥|=(∑^𝑛_𝑖=1𝑥²_𝑖)^1/2; if 𝐴 is a vector or matrix, its transpose is denoted by 𝐴^𝑇; if 𝐴 is matrix, trace norm of matrix 𝐴 is denoted by |𝐴|=trace(𝐴^𝑇𝐴); if 𝐴 is negative matrix, we denote it by 𝐴<0. Suppose that 𝑦(𝑡) is a continuous bounded function on [0,+∞); we define 𝑦_𝑢=max_{𝑡∈[0,+∞)}𝑦(𝑡), 𝑦_𝑙=min_{𝑡∈[0,+∞)}𝑦(𝑡), with usual assumption inf∅=∞, where ∅ denotes the empty set.

3. The Existence and Uniqueness of Global Positive Solution

In this section, we show that model (2) has a unique global solution, and the solution will remain in 𝑅^𝑛₊₊ with probability 1.

Theorem 1. If (𝐻1) and (𝐻2) hold, then there is a unique solution 𝑥(𝑡) to model (2). Moreover, 𝑥(𝑡) remains in 𝑅^𝑛₊₊ with probability 1, where 𝑡∈𝑅.

Proof. Clearly, the coefficients of model (2) satisfy local Lipschitz continuous but do not satisfy the linear growth condition. To show that the solution 𝑥(𝑡) is global, a.s., 𝜏_𝑒=∞ is needed now, where 𝜏_𝑒 is the explosion time.
Let 𝑘₀>0 be sufficiently large such that each component of initial data 𝜉(0) is lying in the interval (1/𝑘₀,𝑘₀). For each integer 𝑘≥𝑘₀, define the stopping time 𝜏_𝑘=inf𝑡∈0,𝜏_𝑒∶𝑥_𝑖∉1𝑘,.,𝑘forsome𝑖=1,2,…,𝑛(4) Clearly, 𝜏_𝑘 is increasing as 𝑘→∞. Set 𝜏_∞=lim_𝑘→∞𝜏_𝑘; hence, 𝜏_∞≤𝜏_𝑒 a.s.. If we can prove that 𝜏_∞=∞ a.s., then 𝜏_𝑒=∞ a.s.. To prove this statement, let us define a 𝐶²-function 𝑉₁∶𝑅^𝑛₊₊→𝑅₊₊ by 𝑉₁(𝑥)=∑^𝑛_𝑖=1𝑥^𝑝_𝑖−𝑝log𝑥_𝑖, where 𝑝∈(0,1). It is easy to see that 𝑉₁(𝑥)>0 for all 𝑥∈𝑅^𝑛₊₊. Applying the Itô formula to model (2), it leads to 𝑑𝑉₁(𝑥(𝑡))=𝐿𝑉₁(𝑥(𝑡))𝑑𝑡+^𝑛_𝑖=1𝑝𝑥^𝑝_𝑖(𝑡)−1𝑔_𝑖𝑑𝐵_𝑖(𝑡),(5) where 𝐿𝑉₁(𝑥)=𝑝^𝑛_𝑖=1𝑥^𝑝_𝑖−1𝑓_𝑖+12𝑝^𝑛_𝑖=1𝑝−1𝑥^𝑝_𝑖+1𝑔²_𝑖.(6) From the elementary inequality 𝑎𝑏≤(1/2)(𝑎²+𝑏²) and Hölder's inequality, we have ^𝑛_𝑖=1𝑥^𝑝_𝑖−1𝑓_𝑖≤^𝑛_𝑖=1ℎ_𝑖𝑢𝑥^𝑝_𝑖+1+12^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢𝑥^𝑝_𝑖−1²+12^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢𝑥²_𝑗+12^𝑛𝑛_𝑖=1𝑗=1𝑏_𝑖𝑗²_𝑢𝑥^𝑝_𝑖−1²+12^𝑛𝑛_𝑖=1𝑗=1𝑥²_𝑗((𝑡+1𝑡−𝛿))2^𝑛𝑛_𝑖=1𝑗=1𝑐_𝑖𝑗²_𝑢𝑥^𝑝_𝑖−1²+12^𝑛𝑛_𝑖=1𝑗=1⁰_−∞𝑥²_𝑗(𝑡+𝜃)𝑑𝜇(𝜃).(7) By the fact that 𝐾=(‖𝜉‖_𝐶𝑟+‖𝑥‖_𝐶𝑟)²<∞, 𝑢²/𝜌_𝑗−𝑣²/(𝜌_𝑗−1)≤(𝑢+𝑣)²≤𝑢²/𝑙_𝑖+𝑣²/(1−𝑙_𝑖), 0<𝑙_𝑖<1<𝜌_𝑗, and (∑^𝑛_𝑖=1𝑎_𝑖𝑏_𝑖)²≤∑^𝑛_𝑖=1𝑎²_𝑖∑^𝑛_𝑖=1𝑏²_𝑖, it implies that 𝑔²_𝑖≤𝑥_𝑖−𝑥^∗_𝑖²1𝑙₁𝜎²_𝑖𝑢+1𝑙₂1−𝑙₁×^𝑛_𝑗=1𝛼_𝑖𝑗²_𝑢^𝑛_𝑗=1𝑥²_𝑗+1𝑙₃1−𝑙₂1−𝑙₁^𝑛_𝑗=1𝛽_𝑖𝑗²_𝑢×^𝑛_𝑗=1𝑥²_𝑗((𝑡+1𝑡−𝛿))1−𝑙₃1−𝑙₂1−𝑙₁^𝑛_𝑗=1𝛾_𝑖𝑗²_𝑢×^𝑛_𝑗=1⁰_−∞𝑥²_𝑗()(𝜃)≤𝑥𝑡+𝜃𝑑𝜇_𝑖−𝑥^∗_𝑖²1𝑙₁𝜎²_𝑖𝑢+𝐾𝑙₂1−𝑙₁^𝑛_𝑗=1𝛼_𝑖𝑗²_𝑢+𝐾𝑙₃1−𝑙₂1−𝑙₁^𝑛_𝑗=1𝛽_𝑖𝑗²_𝑢+𝐾1−𝑙₃1−𝑙₂1−𝑙₁^𝑛_𝑗=1𝛾_𝑖𝑗²_𝑢∶=𝐾₁𝑥_𝑖−𝑥^∗_𝑖²,𝑔²_𝑖≥𝑥_𝑖−𝑥^∗_𝑖²𝜎²_𝑖𝑙+^𝑛_𝑗=1𝛼_𝑖𝑗(𝑡)𝑥_𝑗²+^𝑛_𝑗=1𝛽_𝑖𝑗(𝑡)𝑥_𝑗(𝑡−𝛿(𝑡))²+^𝑛_𝑗=1𝛾_𝑖𝑗(𝑡)⁰_−∞𝑥²_𝑗(𝑡+𝜃)𝑑𝜇(𝜃)²≥𝜎²_𝑖𝑙𝑥_𝑖−𝑥^∗_𝑖²∶=𝐾₂𝑥_𝑖−𝑥^∗_𝑖²,(8) where 0<𝐾₂<𝐾₁ for any 𝑙_𝑖∈(0,1). Then (6) yields that 𝐿𝑉₁(𝑥)≤^𝑛_𝑖=1𝑝ℎ_𝑖𝑢𝑥^𝑝_𝑖+𝑝+12^𝑛_𝑗=1𝑎_𝑖𝑗𝑢+𝑏_𝑖𝑗²_𝑢+𝑐_𝑖𝑗²_𝑢𝑥^𝑝_𝑖−1²+𝑝2^𝑛_𝑗=1𝑎_𝑖𝑗𝑢𝑥²_𝑗+𝑝𝐾₁2𝑥_𝑖−𝑥^∗_𝑖²−𝑝1−𝑝𝐾₂2𝑥^𝑝_𝑖𝑥_𝑖−𝑥^∗_𝑖²+𝑝2^𝑛𝑛_𝑖=1𝑗=1𝑥²_𝑗((𝑡+𝑝𝑡−𝛿))2^𝑛𝑛_𝑖=1𝑗=1⁰_−∞𝑥²_𝑗(𝑡+𝜃)𝑑𝜇(𝜃).(9) Define 𝑉₂(𝑥(𝑡))=𝑉₁(𝑥(𝑡))+(𝑝/2)∑^𝑛_𝑖=1∑^𝑛_𝑗=1∫^𝑡_{𝑡−𝛿(𝑡)}𝑥²_𝑗(𝑠)𝑑𝑠+(𝑝/2)∑^𝑛_𝑖=1∑^𝑛_𝑗=1∫⁰_−∞∫^𝑡_𝑡+𝜃𝑥²_𝑗(𝑠)𝑑𝑠𝑑𝜇(𝜃); the proof is easily checked (details can be found at the appendix).

4. Global Asymptotic Stability

Theorem 2. If (𝐻1) and (𝐻2) hold, there exist positive numbers 𝑑₁,𝑑₂,…,𝑑_𝑛 such that 𝐷𝐴+𝐴^𝑇𝐷 is negative definite, and then lim_𝑡→∞𝑥_𝑖(𝑡)=𝑥^∗_𝑖𝑎.𝑠.,𝑖=1,2,…,𝑛;(10) that is, 𝑥^∗ is globally asymptotically stable a.s., where 𝐾=(‖𝜉‖_𝐶𝑟+‖𝑥‖_𝐶𝑟)², 𝐷=diag(𝑑₁,𝑑₂,…,𝑑_𝑛), 𝐴=(𝜂_𝑖𝑗)_𝑛×𝑛, 𝜂_𝑖𝑗=𝑎_𝑖𝑗 for 𝑖≠𝑗, and 𝜂_𝑖𝑖=𝑎_𝑖𝑖𝑢+0.5𝑥^∗_𝑖𝜎²_𝑖𝑢+𝐾^𝑛_𝑗=1𝛼_𝑖𝑗²_𝑢+𝐾^𝑛_𝑗=1𝛽_𝑖𝑗²_𝑢+𝐾^𝑛_𝑗=1𝛾_𝑖𝑗²_𝑢+0.5^𝑛_𝑗=1𝑏_𝑖𝑗𝑢+^𝑛_𝑗=1𝑐_𝑖𝑗𝑢+^𝑛_𝑗=1𝑑_𝑗𝑑_𝑖𝑏_𝑗𝑖𝑢+^𝑛_𝑗=1𝑑_𝑗𝑑_𝑖𝑐_𝑗𝑖𝑢.(11)

Proof. Define a 𝐶²-function 𝑉₃∶𝑅^𝑛₊₊→𝑅₊₊ by 𝑉₃(𝑥)=∑^𝑛_𝑖=1𝑑_𝑖(𝑥_𝑖−𝑥^∗_𝑖−𝑥^∗_𝑖log(𝑥_𝑖/𝑥^∗_𝑖)). Similar to the proof of [6] (see Theorem 2.1), we can derive 𝑑𝑉₃(𝑥(𝑡))≤𝐿𝑉₃(𝑥(𝑡))𝑑𝑡+^𝑛_𝑖=1𝑑_𝑖𝑥_𝑖(𝑡)−𝑥^∗_𝑖𝑔_𝑖𝑑𝐵_𝑖(𝑡)≤12𝑥_𝑖(𝑡)−𝑥^∗_𝑖𝐷𝐴+𝐴^𝑇𝐷𝑥_𝑖(𝑡)−𝑥^∗_𝑖^𝑇+𝑑𝑡_𝑛𝑖=1𝑑_𝑖𝑥_𝑖(𝑡)−𝑥^∗_𝑖𝑔_𝑖𝑑𝐵_𝑖(𝑡).(12) Since 𝐷𝐴+𝐴^𝑇𝐷 is negative definite, 𝐿𝑉₃(𝑥)<0 is valid along trajectories in 𝑅^𝑛₊₊ except 𝑥^∗. The proof is complete.

5. Approximation Distribution of Solution

Theorem 3. If  (𝐻1), (𝐻2), and the conditions of Theorem 2 hold, then each component 𝑥_𝑖(𝑡) of solution 𝑥(𝑡) to model (2) follows asymptotically 1-dimensional normal distribution 𝒩(𝑥^∗_𝑖,+∞), where 𝑖=1,2,…,𝑛.

Proof. By the definition of equilibrium state 𝑥^∗_𝑖≠0, then 𝑓_𝑖(𝑥^∗_𝑖)=0. Since 𝑥^∗ is stable for the deterministic part to model (2), then 𝜕_𝑥𝑖(𝑓_𝑖(𝑥)𝑥_𝑖)|_{𝑥𝑖=𝑥}^∗_𝑖∶=𝑚_𝑖<0 for 𝑖=1,2,…,𝑛.
Linearizing the 𝑖th equation of model (2) by Taylor expansion at 𝑥^∗_𝑖, then we have 𝑑𝑥_𝑖−𝑥^∗_𝑖≈𝑓_𝑖𝑥^∗_𝑖𝑥^∗_𝑖+𝜕_𝑥𝑖𝑓_𝑖(𝑥)𝑥_𝑖|_{𝑥𝑖=𝑥}^∗_𝑖×𝑥_𝑖−𝑥^∗_𝑖+𝑜𝑥_𝑖−𝑥^∗_𝑖²+𝑔𝑑𝑡_𝑖𝑥^∗_𝑖𝑥^∗_𝑖+𝜕_𝑥𝑖𝑔_𝑖(𝑥)𝑥_𝑖|_{𝑥𝑖=𝑥}^∗_𝑖×𝑥_𝑖−𝑥^∗_𝑖+𝑜𝑥_𝑖−𝑥^∗_𝑖²𝑑𝐵_𝑖(𝑡).(13) Denote 𝑥_𝑗(𝑡−𝛿(𝑡))=𝑥_𝑗𝛿(𝑡) and ∫⁰_−∞𝑥_𝑗(𝑡+𝜃)𝑑𝜇(𝜃)=𝑥_𝑗𝜃(𝑡), 𝑗=1,2,⋯,𝑛. Since 𝑔_𝑖(𝑥^∗_𝑖)=0 and 𝜕_𝑥𝑖(𝑔_𝑖(𝑥)𝑥_𝑖)|_{𝑥𝑖=𝑥}^∗_𝑖=𝜎_𝑖𝑥^∗_𝑖+𝑥^∗_𝑖∑^𝑛_𝑗=1(𝛼_𝑖𝑗𝑥_𝑗+𝛽_𝑖𝑗𝑥_𝑗𝛿+𝛾_𝑖𝑗𝑥_𝑗𝜃)|_{𝑥𝑖=𝑥}^∗_𝑖∶=𝑢_𝑖≠0, (13) can be simplified as 𝑑𝑥_𝑖−𝑥^∗_𝑖𝑑𝑡≈𝑚_𝑖𝑥_𝑖−𝑥^∗_𝑖+𝑢_𝑖𝑥_𝑖−𝑥^∗_𝑖𝑑𝐵_𝑖(𝑡)𝑑𝑡,𝑖=1,2,…,𝑛,(14) where 𝑚_𝑖<0, 𝑢_𝑖≠0. Then (14) implies that 𝑥_𝑖(𝑡)−𝑥^∗_𝑖=𝑥_𝑖(0)−𝑥^∗_𝑖𝑒^{𝑚𝑖𝑡+𝑢𝑖𝐵𝑖(𝑡)},𝑖=1,2,…,𝑛.(15) From the definition of 1-dimensional Brownian motion, we can derive that 𝑥_𝑖(𝑡)∼𝒩(𝑥^∗_𝑖+(𝑥_𝑖(0)−𝑥^∗_𝑖)𝑒^𝑚𝑖𝑡, (𝑥_𝑖(0)−𝑥^∗_𝑖)²𝑒^{𝑢2𝑖𝑡}). According to the conditions of Theorem 2, when 𝑡→+∞, one can find that 𝑢_𝑖→𝑢^∗_𝑖=𝜎_𝑖𝑥^∗_𝑖+𝑥^∗_𝑖∑^𝑛_𝑗=1(𝛼_𝑖𝑗𝑥^∗_𝑗+𝛽_𝑖𝑗𝑥^∗_𝑗𝛿+𝛾_𝑖𝑗𝑥^∗_𝑗𝜃)≠0 a.s.; thus, lim_𝑡→+∞𝑥_𝑖(𝑡)∼𝒩(𝑥^∗_𝑖+∞); in other words, 𝑥_𝑖(𝑡) asymptotically follows 1-dimensional normal distribution. When 𝑡 tends to ∞, the mean of solution 𝑥(𝑡) to model (2) is 𝑥^∗; it is just consistent with the conclusion of Theorem 2. However, the deviation between solution 𝑥_𝑖(𝑡) and mean 𝑥^∗_𝑖 may approach infinity; it is bad information for further analysis. Now, if the variance can be evaluated, the disadvantage will be improved.

6. Moment Estimation

First, let us prove one useful moment estimation.

Theorem 4. If (𝐻1) and (𝐻2) hold, then there is a positive constant 𝐺=𝐺(𝑝), such that limsup_𝑡→∞𝐸|𝑥(𝑡)|^𝑝≤𝐺, where 𝑝∈(0,1).

Proof. Define a 𝐶²-function 𝑉₄∶𝑅^𝑛₊₊→𝑅₊₊ by 𝑉₄(𝑥)=∑^𝑛_𝑖=1𝑥^𝑝_𝑖. For any given 𝜀∈(0,2𝑟), applying the Itô formula to 𝑒^𝜀𝑡𝑉₄(𝑥(𝑡)) and taking expectation, it yields that 𝑒^𝜀𝑡𝐸𝑉₄(𝑥(𝑡))=𝐸𝑉₄(𝜉(0))+𝐸^𝑡₀𝑒^𝜀𝑠𝐿𝑉₄(𝑥(𝑠))+𝜀𝑉₄(𝑥(𝑠))𝑑𝑠,(16) where 𝐿𝑉₄(𝑥)=𝑝^𝑛_𝑖=1𝑥^𝑝_𝑖𝑓_𝑖−𝑝1−𝑝2^𝑛_𝑖=1𝑥^𝑝_𝑖𝑔²_𝑖.(17) By fundamental inequality 𝑎𝑏≤(1/2)(𝑎²+𝑏²), ^𝑛_𝑖=1𝑥^𝑝_𝑖𝑓_𝑖≤^𝑛_𝑖=1ℎ_𝑖𝑢𝑥^𝑝_𝑖+12^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢𝑥_𝑖^2𝑝+12^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢𝑥²_𝑗+12^𝑛𝑛_𝑖=1𝑗=1𝑏_𝑖𝑗²_𝑢𝑥_𝑖^2𝑝+12^𝑛𝑛_𝑖=1𝑗=1𝑐_𝑖𝑗²_𝑢𝑥_𝑖^2𝑝+12^𝑛𝑛_𝑖=1𝑗=1𝑥²_𝑗((𝑠+1𝑠−𝛿))2^𝑛𝑛_𝑖=1𝑗=1⁰_−∞𝑥²_𝑗(𝑠+𝜃)𝑑𝜇_𝑛𝑖=1(𝜃),𝑥^𝑝_𝑖𝑔²_𝑖≥^𝑛_𝑖=1𝐾₂𝑥^𝑝_𝑖𝑥_𝑖−𝑥^∗_𝑖².(18) Consequently, 𝐿𝑉₄(𝑥(𝑠))≤𝑝^𝑛_𝑖=1ℎ_𝑖𝑢𝑥^𝑝_𝑖+𝑝2^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢+𝑏_𝑖𝑗²_𝑢+𝑐_𝑖𝑗²_𝑢𝑥_𝑖^2𝑝+𝑝2^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢𝑥²_𝑗+𝑝2^𝑛𝑛_𝑖=1𝑗=1𝑥²_𝑗((𝑠+1𝑠−𝛿))2^𝑛𝑛_𝑖=1𝑗=1⁰_−∞𝑥²_𝑗()(𝜃)−𝑝𝑠+𝜃𝑑𝜇1−𝑝2^𝑛_𝑖=1𝐾₂𝑥^𝑝_𝑖𝑥_𝑖−𝑥^∗_𝑖²∶=𝐻(𝑥)−𝜀𝑉₄(𝑥)+𝑛𝑝2|𝑥(𝑠−𝛿(𝑠))|²−𝑒^𝜀𝜏|𝑥|²+𝑛𝑝2⁰_−∞|𝑥(𝑠+𝜃)|²𝑑𝜇(𝜃)−𝜇_𝑟|𝑥|²,(19) where 𝐻(𝑥)=^𝑛_𝑖=1𝑝ℎ_𝑖𝑢+𝜀𝑥^𝑝_𝑖+𝑝2^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢+𝑏_𝑖𝑗²_𝑢+𝑐_𝑖𝑗²_𝑢𝑥_𝑖^2𝑝+𝑛𝑝2𝑒^𝜀𝜏+𝜇_𝑟|𝑥|²+𝑝2^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢𝑥²_𝑗−𝑝1−𝑝𝐾₂2^𝑛_𝑖=1𝑥^𝑝_𝑖𝑥_𝑖−𝑥^∗_𝑖².(20) Noting that 𝐻(𝑥) is bounded in 𝑅^𝑛₊₊; namely, 𝐾₃∶=sup_{𝑥∈𝑅𝑛++}𝐻(𝑥)<∞. Substituting this into (16), it thus follows that 𝑒^𝜀𝑡𝐸𝑉₄(𝑥(𝑡))≤𝐸𝑉₄(𝜉(0+1))𝜀𝐾₃𝑒^𝜀𝑡−1+𝑛𝑝2×𝐸^𝑡₀𝑒^𝜀𝑠|𝑥(𝑠−𝛿(𝑠))|²−𝑒^𝜀𝜏|𝑥(𝑠)|²+⁰_−∞|𝑥(𝑠+𝜃)|²𝑑𝜇(𝜃)−𝜇_𝑟|𝑥(𝑠)|²𝑑𝑠.(21) By (H1) and (H2), we have 𝐸^𝑡₀𝑒^𝜀𝑠|𝑥(𝑠−𝛿(𝑠))|²−𝑒^𝜀𝜏|𝑥(𝑠)|²𝑑𝑠≤𝐸^𝑡_−𝜏𝑒^{𝜀(𝑠+𝜏)}|𝑥(𝑠)|²𝑑𝑠−^𝑡₀𝑒^{𝜀(𝑠+𝜏)}|𝑥(𝑠)|²𝑑𝑠≤𝐸⁰_−𝜏𝑒^𝜀𝜏|𝑥(𝑠)|²𝐸𝑑𝑠,_𝑡0𝑒^𝜀𝑠0_−∞|𝑥(𝑠+𝜃)|²𝑑𝜇(𝜃)−𝜇_𝑟|𝑥(𝑠)|²𝑑𝑠=𝐸^𝑡₀𝑒^{𝜀𝑠−𝑠}_−∞𝑒^{2𝑟(𝑠+𝜃)}|𝑥(𝑠+𝜃)|²×𝑒^{−2𝑟(𝑠+𝜃)}(𝜃+𝑑𝜇)𝑑𝑠_𝑡0𝑒^𝜀𝑠0_−𝑠|𝑥(𝑠+𝜃)|²𝑑𝜇(𝜃)𝑑𝑠−𝜇_𝑟^𝑡₀𝑒^𝜀𝑠|𝑥(𝑠)|²𝜉𝑑𝑠≤𝐸²_𝐶𝑟^𝑡₀𝑒^{(𝜀−2𝑟)𝑠}𝑑𝑠⁰_−∞𝑒^−2𝑟𝜃(𝜃)+𝑑𝜇_{0−𝑡𝑡−𝜃}𝑒^𝜀𝑠|𝑥(𝑠+𝜃)|²𝑑𝑠𝑑𝜇(𝜃)−𝜇_𝑟^𝑡₀𝑒^𝜀𝑠|𝑥(𝑠)|²𝜉𝑑𝑠≤𝐸²_𝐶𝑟𝜇_𝑟^𝑡₀𝑒^{(𝜀−2𝑟)𝑠}+𝑑𝑠_0−∞𝑡0𝑒^{𝜀(𝑠−𝜃)}|𝑥(𝑠)|²𝑑𝑠𝑑𝜇(𝜃)−𝜇_𝑟^𝑡₀𝑒^𝜀𝑠|𝑥(𝑠)|²𝜉𝑑𝑠≤𝐸²_𝐶𝑟𝜇_𝑟𝑡+⁰_−∞𝑒^−𝜀𝜃(𝜃)×𝑑𝜇_𝑡0𝑒^𝜀𝑠|𝑥(𝑠)|²𝑑𝑠−𝜇_𝑟^𝑡₀𝑒^𝜀𝑠|𝑥(𝑠)|²𝜉𝑑𝑠≤𝐸²_𝐶𝑟𝜇_𝑟𝑡.(22) Hence, we derive 𝑒^𝜀𝑡𝐸𝑉₄(𝑥(𝑡))≤𝐸𝑉₄(𝜉(0+𝐾))³𝑒^𝜀𝑡𝜀+𝑛𝑝2𝐸𝜇_𝑟𝜉²_𝐶𝑟𝑡+𝐸⁰_−𝜏𝑒^𝜀𝜏|𝑥(𝑠)|²𝑑𝑠.(23) This implies that limsup_𝑡→∞𝐸𝑉₄(𝑥(𝑡))≤𝜀⁻¹𝐾₃.
Again, by the fact that |𝑥|²≤𝑛max_{1≤𝑖≤𝑛}𝑥²_𝑖, it then follows that |𝑥|^𝑝≤𝑛^𝑝/2max_{1≤𝑖≤𝑛}𝑥^𝑝_𝑖≤𝑛^𝑝/2𝑉₄(𝑥), and we have limsup_𝑡→∞𝐸|𝑥(𝑡)|^𝑝≤𝐺𝑝;(24) the assertion follows by setting 𝐺=𝐺(𝑝)=𝑛^𝑝/2𝜀⁻¹𝐾₃. The proof is complete.

Remark A. If (H1) and (H2) are valid, then limsup_𝑡→∞𝑃{|𝑥(𝑡)|≤𝜒}≥1−𝜀; that is, solution 𝑥(𝑡) to model (2) is stochastically ultimately bounded.

If we take 𝑝=1/2, then the result is valid by using of Theorem 4 and Chebyshev's inequality. The proof is omitted herewith.

Theorem 5. If (𝐻1) and (𝐻2) hold, there exists a positive constant 𝑄, and then limsup_𝑇→∞1𝑇^𝑇₀𝐷𝑥_𝑖(𝑡)𝑑𝑡≤𝑄,𝑖=1,2,…,𝑛.(25)

Proof. Rewrite (20) as 𝐻(𝑥)=𝐻₁(𝑥)−𝑛𝑝+1|𝑥|²+𝑛𝑝2𝑒^𝜀𝜏+𝜇_𝑟|𝑥|²+𝜀𝑉₄(𝑥),(26) with 𝐻₁(𝑥)=^𝑛_𝑖=1𝑝ℎ_𝑖𝑢𝑥^𝑝_𝑖+𝑝2^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢+𝑏_𝑖𝑗²_𝑢+𝑐_𝑖𝑗²_𝑢𝑥_𝑖^2𝑝+𝑝2^𝑛𝑛_𝑖=1𝑗=1𝑎_𝑖𝑗𝑢𝑥²_𝑗+𝑛𝑝+1|𝑥|²−𝑝1−𝑝2^𝑛_𝑖=1𝐾₂𝑥^𝑝_𝑖𝑥_𝑖−𝑥^∗_𝑖².(27) Clearly, 𝐻₁(𝑥) is bounded in 𝑅^𝑛₊₊; namely, 𝐾₄∶=sup_{𝑥∈𝑅𝑛++}𝐻₁(𝑥)<∞. So 𝐻(𝑥)≤𝐾₄−𝑛𝑝+1|𝑥|²+𝑛𝑝2𝑒^𝜀𝜏+𝜇_𝑟|𝑥|²+𝜀𝑉₄(𝑥).(28) By (19), we have 𝑑𝑉₄(𝑥(𝑡))≤𝐾₄−𝑛𝑝+1|𝑥(𝑡)|²+𝑛𝑝2|𝑥(𝑡−𝛿(𝑡))|²+⁰_−∞|𝑥(𝑡+𝜃)|²𝑑𝜇(𝜃)𝑑𝑡+𝑝^𝑛_𝑖=1𝑥^𝑝_𝑖(𝑡)𝑔_𝑖𝑑𝐵_𝑖(𝑡).(29) Integrating both sides of (29) from 0 to 𝑇 (𝑇>0 is arbitrary) and then taking expectations, we obtain that 𝐸𝑉₄(𝑥(𝑇))≤𝑉₄(𝜉(0))+𝐾₄𝑇−𝐸^𝑇₀|𝑥(𝑡)|²𝑑𝑡+𝑛𝑝2𝐼₁+𝑛𝑝2𝐼₂.(30) Similar to (22), by the definition of  𝛿(⋅), we can derive 𝐼₁=𝐸^𝑇₀|𝑥(𝑡−𝛿(𝑡))|²−|𝑥(𝑡)|²𝑑𝑡≤⁰_−𝜏|𝑥(𝑡)|²𝐼𝑑𝑡,₂=𝐸^𝑇₀⁰_−∞|𝑥(𝑡+𝜃)|²𝑑𝜇(𝜃)−|𝑥(𝑡)|²1𝑑𝑡≤𝜇2𝑟^𝑟𝜉²_𝐶𝑟.(31) Substituting (31) into (30), then we have limsup_𝑇→∞1𝑇^𝑇₀𝐸|𝑥(𝑡)|²𝑑𝑡≤𝐾₄∶=𝑄.(32) By the fact that 𝐷(𝑥_𝑖)=𝐸|𝑥_𝑖|²−(𝐸𝑥_𝑖)² and lim_𝑡→∞𝐸(𝑥(𝑡))=𝑥^∗, the conclusion is obviously valid. The proof is complete.

7. Example and Numerical Simulation

Now, we will simulate asymptotic behaviors of solutions to model (2). Let us consider the following 1-dimensional stochastic autonomous Logistic model: 𝑑𝑥(𝑡)=𝑥(𝑡)ℎ−𝑎𝑥(𝑡)−𝑏𝑥(𝑡−𝜏)−𝑐⁰_−∞𝑥()(𝜃(𝑡𝜎𝑥(𝑡)𝑡+𝜃𝑑𝜇)𝑑𝑡+𝑥)−𝑥^∗+𝛼𝑥(𝑡)𝑥(𝑡)−𝑥^∗+𝛽𝑥(𝑡−𝜏)𝑥(𝑡)−𝑥^∗+𝛾⁰_−∞𝑥()(𝜃)×𝑥(𝑡)𝑡+𝜃𝑑𝜇−𝑥^∗𝑑𝐵(𝑡),(33) with initial data 𝜉(𝜃)∈𝐶_𝑟, where ℎ, 𝑎, 𝑏, 𝑐, 𝜏, 𝜎, 𝛼, 𝛽, and 𝛾 are all positive numbers, and the equilibrium state is 𝑥^∗=ℎ/(𝑎+𝑏+𝑐). By Theorem 2, if 𝑎>0.5𝑥^∗(𝜎²+𝐾𝛼²+𝐾𝛽²+𝐾𝛾²), then the equilibrium state 𝑥^∗ of model (33) is globally asymptotically stable a.s.. According to Milstein method mentioned in Higham [11], the initial data of model (33) is given by 𝜉(𝜃)=0.32𝑒^0.5𝜃 and 𝜇(𝜃)=𝑒^𝜃, 𝜃∈(−∞,0]; the difference equation is followed next 𝑥^(𝑘+1)−𝑥^(𝑘)=𝑥^(𝑘)ℎ−𝑎𝑥^(𝑘)−𝑏𝑥^{(𝑘−𝜏/Δ𝑡)}−2𝑐3𝑒^{−𝑘Δ𝑡}−𝑐𝑒^{𝑘−𝑘Δ𝑡}_𝑖=1𝑥^(𝑖)𝑒^𝑖Δ𝑡Δ𝑡Δ𝑡+𝜎𝑥^(𝑘)𝑥^(𝑘)−𝑥^∗√Δ𝑡𝜉^(𝑘)+𝜎²2𝑥^(𝑘)𝑥^(𝑘)−𝑥^∗𝜉^(𝑘)2−1Δ𝑡+𝛼𝑥^(𝑘)𝑥^(𝑘)𝑥^(𝑘)−𝑥^∗√Δ𝑡𝜉^(𝑘)+𝛼²2𝑥^(𝑘)𝑥^(𝑘)𝑥^(𝑘)−𝑥^∗𝜉^(𝑘)2−1Δ𝑡+𝛽𝑥^(𝑘)𝑥^{(𝑘−𝜏/Δ𝑡)}𝑥^(𝑘)−𝑥^∗√Δ𝑡𝜉^(𝑘)+𝛽²2𝑥^(𝑘)𝑥^{(𝑘−𝜏/Δ𝑡)}𝑥^(𝑘)−𝑥^∗𝜉^(𝑘)2−1Δ𝑡+𝛾2𝑐3𝑒^{−𝑘Δ𝑡}+𝑐𝑒^{𝑘−𝑘Δ𝑡}_𝑖=1𝑥^(𝑖)𝑒^𝑖Δ𝑡Δ𝑡𝑥^(𝑘)×𝑥^(𝑘)−𝑥^∗√Δ𝑡𝜉^(𝑘)+𝛾²22𝑐3𝑒^{−𝑘Δ𝑡}+𝑐𝑒^{𝑘−𝑘Δ𝑡}_𝑖=1𝑥^(𝑖)𝑒^𝑖Δ𝑡Δ𝑡𝑥^(𝑘)×𝑥^(𝑘)−𝑥^∗𝜉^(𝑘)2−1Δ𝑡,(34) where 𝜉^(𝑘)(𝑘=1,2,…,𝑛) is the Gaussian random variable which follows the standard normal distribution 𝒩(0,1). If we fix the parameters ℎ=0.5, 𝑎=0.6, 𝑏=0.1, 𝑐=0.3, 𝜏=1, and Δ𝑡=0.01, then 𝑥^∗=ℎ/(𝑎+𝑏+𝑐)=0.5. When 𝜎=𝛼=𝛽=𝛾=0, model (33) becomes a deterministic one; then the equilibrium state 𝑥^∗ is globally asymptotically stable a.s. by Theorem 2 (see Figure 1). When 𝜎=0.1 and 𝛼=𝛽=𝛾=0, the intrinsic growth rate of model (33) is perturbed by white noise; if 𝑎>0.5𝑥^∗𝜎²+𝑏+𝑐, the equilibrium state 𝑥^∗ is globally asymptotically stable a.s. (see Figure 2). When 𝜎=0.1 and 𝛼=𝛽=𝛾=0.01, that is, each parameter of model (33) is fluctuated by white noise, if 𝑎>0.5𝑥^∗(𝜎²+𝐾𝛼²+𝐾𝛽²+𝐾𝛾²)+𝑏+𝑐 and 𝐾=(0.32𝑒^{(𝑟+0.5)𝜃}+𝑥^(𝑘))²≤(0.32+𝑥^(𝑘))² are satisfied, then the equilibrium stable 𝑥^∗ is globally asymptotically stable a.s. (see Figure 3).

Figure 1 

Model (33) becomes a deterministic model.

Figure 2 

The intrinsic growth rate of model (33) is perturbed by white noise.

Figure 3 

Each parameter of model (33) is fluctuated by white noise.

Appendix

From the definition of 𝑉₂(𝑡), we obtain 𝑑𝑉₂(𝑥(𝑡))≤𝐹(𝑥(𝑡))𝑑𝑡+^𝑛_𝑖=1𝑝𝑥^𝑝_𝑖(𝑡)−1𝑔_𝑖𝑑𝐵_𝑖(𝑡),(A.1) where 𝐹(𝑥)=^𝑛_𝑖=1𝑝ℎ_𝑖𝑢𝑥^𝑝_𝑖+𝑝+12^𝑛_𝑗=1𝑎_𝑖𝑗𝑢+𝑏_𝑖𝑗²_𝑢+𝑐_𝑖𝑗²_𝑢𝑥^𝑝_𝑖−1²+𝑝2^𝑛_𝑗=1𝑎_𝑖𝑗𝑢+2𝑥²_𝑗+𝑝𝐾₁2𝑥_𝑖−𝑥^∗_𝑖²−𝑝1−𝑝𝐾₂2𝑥^𝑝_𝑖𝑥_𝑖−𝑥^∗_𝑖².(A.2) It is easy to check that 𝐹(𝑥)≤𝐾^′<∞. Equation (A.1) becomes 𝑑𝑉₂(𝑥(𝑡))≤𝐾^′𝑑𝑡+^𝑛_𝑖=1𝑝𝑥^𝑝_𝑖(𝑡)−1𝑔_𝑖𝑑𝐵_𝑖(𝑡).(A.3) Integrating both sides of (A.3) from 0 to 𝜏_𝑘∧𝑇 (𝑇>0 is arbitrary) and then taking expectations, it yields that 𝐸𝑉₁𝑥𝜏_𝑘∧𝑇≤𝐸𝑉₂𝑥𝜏_𝑘∧𝑇≤𝑉₂(𝜉(0))+𝐾^′𝑇.(A.4) Noting that for every 𝑤∈Ω_𝑘={𝜏_𝑘≤𝑇}, by the definition of stopping time 𝜏_𝑘, 𝑥_𝑖(𝜏_𝑘,𝑤)=𝑘 or 1/𝑘 for some 𝑖=1,2,…,𝑛, 𝑉₁(𝑥(𝜏_𝑘∧𝑇))≥min{𝑘^𝑝−𝑝log𝑘,𝑘^−𝑝+𝑝log𝑘}. It then follows that 𝑃𝜏_𝑘≤𝑇𝑘^𝑝−𝑝log𝑘∧𝑘^−𝑝𝐼+𝑝log𝑘≤𝐸_{{𝜏𝑘≤𝑇}}(𝑤)𝑉₂𝑥𝜏_𝑘,𝑤≤𝑉₂(𝜉(0))+𝐾^′𝑇.(A.5)

Thus, lim_𝑘→∞𝑃{𝜏_𝑘≤𝑇}=0. Since 𝑇>0 is arbitrary, we must have 𝑃{𝜏_∞<∞}=0; then 𝑃{𝜏_∞=∞}=1 as required.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the Natural Science Foundation of China (no. 11201075), the Natural Science Foundation of Fujian Province of China (no. 2010J01005), and the Technology Innovation Platform Project of Fujian Province (no. 2009J1007).