Dynamical Behavior of the Stochastic Delay Mutualism System

Abstract

We discuss the dynamical behavior of the stochastic delay three-specie mutualism system. We develop the technique for stochastic differential equations to deal with the asymptotic property. Using it we obtain the existence of the unique positive solution, the asymptotic properties, and the nonpersistence. Finally, we give the numerical examinations to illustrate our results.

1. Introduction

The classical Lotka-Volterra model for two mutualistic species is described by the ordinary differential equation (ODE) 𝑥̇₁(𝑡)=𝑥₁(𝑡)𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡,𝑥̇₂)(𝑡)=𝑥₂(𝑡)𝑟₂+𝑎₂₁𝑥₁(𝑡)−𝑎₂₂𝑥₂(𝑡). (1) There are many extensive literatures concerned with the dynamics of this model and we do not mention them here except [1]. Goh [1] showed that if 𝑎₁₁𝑎₂₂>𝑎₁₂𝑎₂₁ holds, system (1) has a stable and globally attractive equilibrium point 𝑥^∗=(𝑥^∗₁,𝑥^∗₂) with the following property: 𝑥₁(𝑡)⟶𝑥^∗₁,𝑥₂(𝑡)⟶𝑥^∗₂,as𝑡⟶∞, (2) where 𝑥^∗₁=𝑟₁𝑎₂₂+𝑟₂𝑎₁₂𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁>0,𝑥^∗₂=𝑟₂𝑎₁₁+𝑟₁𝑎₂₁𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁>0. (3)

In fact, in many physical as well as biological systems, many studies indicate that time delay widely exists in nature, for examples, in [2–7]. When the growth rate of each specie is affected by the time delay, as a result, (1) becomes a delay differential equation (DDE) 𝑥̇₁(𝑡)=𝑥₁(𝑡)𝑟₁−𝑎₁₁𝑥₁(𝑡−𝜏)+𝑎₁₂𝑥₂(,𝑥̇𝑡−𝜏)₂(𝑡)=𝑥₂(𝑡)𝑟₂+𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₂𝑥₂(𝑡−𝜏). (4) In [4], He and Gopalsamy obtained a supercritical Hopf-bifurcation of (4) at 𝜏=𝜏^∗ (a constant) and proved that the equilibrium (𝑥^∗₁,𝑥^∗₂) is no longer asymptotically stable as the delay increases to 𝜏^∗. If only the interspecific positive feedback terms are affected by the delay, (1) becomes 𝑥̇₁(𝑡)=𝑥₁(𝑡)𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(,𝑥̇𝑡−𝜏)₂(𝑡)=𝑥₂(𝑡)𝑟₂+𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₂𝑥₂(𝑡). (5) The positive equilibrium (𝑥^∗₁,𝑥^∗₂) of (5) is globally attractive if 𝑎₁₁𝑎₂₂>𝑎₁₂𝑎₂₁ holds, which implies the delay is harmless.

Population systems are often subject to environmental noise and many authors have investigated the dynamical behaviors of stochastic population systems, for examples, in [8–26]. May [27] revealed that the parameters of the stochastic systems always fluctuate around their average values and the solution also fluctuates around its average value. If we still use 𝑟_𝑖 to denote the average growth rate, then the intrinsic growth rate becomes 𝑟_𝑖⟶𝑟_𝑖+𝜎_𝑖𝐵̇_𝑖(𝑡), (6) where 𝐵̇_𝑖(𝑡) is white noise and 𝜎_𝑖 is a positive constant representing the intensity of the noise. As a result, the mutualism system (1) becomes a stochastic differential equation (SDE) 𝑥̇₁(𝑡)=𝑥₁(𝑡)𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡)𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡,𝑥̇₂)(𝑡)=𝑥₂(𝑡)𝑟₂+𝑎₂₁𝑥₁(𝑡)−𝑎₂₂𝑥₂(𝑡)𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡). (7) Ji et al. in [14] analyzed the long-time asymptotic behavior of the system (7) and obtained the ergodic property and its stationary distribution.

Let us take a further step by considering a 3-dimensional mutualism system 𝑥̇₁(𝑡)=𝑥₁(𝑡)𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡−𝜏)+𝑎₁₃𝑥₃(,𝑥̇𝑡−𝜏)₂(𝑡)=𝑥₂(𝑡)𝑟₂+𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₂𝑥₂(𝑡)+𝑎₂₃𝑥₃(,𝑥̇𝑡−𝜏)₃(𝑡)=𝑥₃(𝑡)𝑟₃+𝑎₃₁𝑥₁(𝑡−𝜏)+𝑎₃₂𝑥₂(𝑡−𝜏)−𝑎₃₃𝑥₃(𝑡), (8) subject to the white noise. As a result, it becomes a stochastic delay differential equation (SDDE) 𝑑𝑥₁(𝑡)=𝑥₁(𝑡)𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡−𝜏)+𝑎₁₃𝑥₃(𝑡−𝜏)𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡),𝑑𝑥₂(𝑡)=𝑥₂(𝑡)𝑟₂+𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₂𝑥₂(𝑡)+𝑎₂₃𝑥₃(𝑡−𝜏)𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡),𝑑𝑥₃(𝑡)=𝑥₃(𝑡)𝑟₃+𝑎₃₁𝑥₁(𝑡−𝜏)+𝑎₃₂𝑥₂(𝑡−𝜏)−𝑎₃₃𝑥₃(𝑡)𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡). (9)

Our aim is to investigate the long-time asymptotic behavior of SDDE (9). This paper is organized as follows. In order to obtain better dynamic properties of SDDE (9), we show that there exists a unique global positive solution with any initial positive value under some assumptions in Section 2. Then, we estimate the expectation in time average of the distance between the solution of (9) and the positive equilibrium point of the deterministic model (8); namely, limsup_𝑡→∞1𝑡^𝑡₀𝐸𝑥(𝑡)−𝑥^∗≤𝐾(𝜎), (10) where 𝑥^∗ is the unique positive equilibrium point of system (8). In Section 3, we prove that system (9) is persistent in time average as the intensity of the white noise is small and yields the limit of the solution in time average. In Section 4, we obtain the nonpersistence of system (9) as the intensity of noise is big. Finally, in Section 5, we illustrate our results by some numerical examinations.

Throughout this paper, unless otherwise specified, let (Ω,{ℱ_𝑡}_𝑡≥0,𝑃) denote a complete probability space with a filtration {ℱ_𝑡}_𝑡≥0 satisfying the usual conditions (i.e., it is right continuous and ℱ₀ contains all 𝑃-null sets). Denote by 𝑅³₊ the positive cone in 𝑅³; namely, 𝑅³₊={𝑥=(𝑥₁,𝑥₂,𝑥₃)^𝑇∈𝑅³∶𝑥_𝑖>0,𝑖=1,2,3}. Let 𝜏>0 and denote by 𝐶([−𝜏,0];𝑅³₊) the family of the continuous functions from [−𝜏,0] to 𝑅³₊. For 𝑥∈𝑅³, its norm is denoted by |𝑥|=|𝑥₁|+|𝑥₂|+|𝑥₃|.

To discuss the dynamical behavior of (9), we impose the following assumption.

Assumption 1. Consider 𝑎_𝑖𝑗>0,𝑖,𝑗=1,2,3 and 𝑎₁₁>𝑎₁₂+𝑎₁₃, 𝑎₂₂>𝑎₂₁+𝑎₂₃, 𝑎₃₃>𝑎₃₁+𝑎₃₂.

Assumption 2. Consider 𝑟_𝑖>𝜎²_𝑖/2,𝑖=1,2,3.

2. Existence and Uniqueness of the Positive Solution

In population dynamics, the existence of the global positive solution is necessary. In order for a SDE to have a unique global (i.e., no explosion at any finite time) solution for any given initial value, its coefficients are generally required to satisfy the linear growth condition and local Lipschitz condition (Arnold et al. [28], Mao [23]). However, the coefficients of system (9) do not satisfy the linear growth condition, though they are locally Lipschitz continuous, so the solution of system (9) may explode at a finite time. So we prepare the following useful lemma and then yield the existence of the positive solution by using it.

Lemma 3. Under Assumption 1, then we have (1)𝑎₂₂𝑎₃₃−𝑎₂₃𝑎₃₂>0,𝑎₁₁𝑎₃₃−𝑎₁₃𝑎₃₁𝑎>0,₁₁𝑎₂₂−𝑎₁₂𝑎₂₁(2)𝑎>0,₂₂𝑎₃₃−𝑎₂₃𝑎₃₂+𝑎₃₃𝑎₂₁+𝑎₃₁𝑎₂₃+𝑎₃₂𝑎₂₁+𝑎₃₁𝑎₂₂(3)𝑎>0,₁₁𝑎₃₃−𝑎₁₃𝑎₃₁+𝑎₃₃𝑎₁₂+𝑎₃₂𝑎₁₃+𝑎₁₁𝑎₃₂+𝑎₃₁𝑎₁₂(4)𝑎>0,₁₁𝑎₂₂−𝑎₁₂𝑎₂₁+𝑎₁₂𝑎₂₃+𝑎₂₂𝑎₁₃+𝑎₁₁𝑎₂₃+𝑎₂₁𝑎₁₃(5)>0,−𝑎₁₁𝑎₁₂𝑎₁₃𝑎₂₁−𝑎₂₂𝑎₂₃𝑎₃₁𝑎₃₂−𝑎₃₃=−𝑎₁₁𝑎₂₂𝑎₃₃+𝑎₁₁𝑎₂₃𝑎₃₂+𝑎₂₂𝑎₁₃𝑎₃₁+𝑎₃₃𝑎₁₂𝑎₂₁+𝑎₁₂𝑎₂₃𝑎₃₁+𝑎₂₁𝑎₃₂𝑎₁₃<0. (11)

Proof. The assertions (1)--(5) are obviously obtained from 𝑎₁₁>𝑎₁₂+𝑎₁₃, 𝑎₂₂>𝑎₂₁+𝑎₂₃, 𝑎₃₃>𝑎₃₁+𝑎₃₂, and 𝑎_𝑖𝑗>0,𝑖,𝑗=1,2,3. Here we omit the proof.

Theorem 4. Under Assumption 1, for any given initial value 𝑥(⋅)∈𝐶([−𝜏,0];𝑅³₊), there is a unique positive solution 𝑥(𝑡) to system (9) on 𝑡≥−𝜏 and the solution will remain in 𝑅³₊ with probability  1; namely, 𝑥(𝑡)∈𝑅³₊  for all 𝑡≥−𝜏 a. s.

Proof. Firstly, define a 𝐶² function 𝑉₁∶𝑅³₊→𝑅₊ by 𝑉₁𝑥₁,𝑥₂,𝑥₃=𝑐₁𝑥₁(𝑡)−1−log𝑥₁(𝑡)+𝑐₂𝑥₂(𝑡)−1−log𝑥₂(𝑡)+𝑐₃𝑥₃(𝑡)−1−log𝑥₃(𝑡), (12) where 𝑐₁=𝑎₂₂𝑎₃₃−𝑎₂₃𝑎₃₂+𝑎₃₃𝑎₂₁+𝑎₃₁𝑎₂₃+𝑎₃₂𝑎₂₁+𝑎₃₁𝑎₂₂×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹,𝑐₂=𝑎₁₁𝑎₃₃−𝑎₁₃𝑎₃₁+𝑎₃₃𝑎₁₂+𝑎₃₂𝑎₁₃+𝑎₁₁𝑎₃₂+𝑎₃₁𝑎₁₂×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹,𝑐₃=𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁+𝑎₁₂𝑎₂₃+𝑎₂₂𝑎₁₃+𝑎₁₁𝑎₂₃+𝑎₂₁𝑎₁₃×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹. (13) From Lemma 3, 𝑐_𝑖>0,𝑖=1,2,3. So, 𝑉₁(𝑥₁,𝑥₂,𝑥₃) is nonnegative. By Itô’s formula, we obtain 𝑑𝑉₁𝑥₁,𝑥₂,𝑥₃=𝐿𝑉₁𝑥₁,𝑥₂,𝑥₃+𝑑𝑡_3𝑖=1𝑐_𝑖𝑥_𝑖(𝑡)−1𝜎_𝑖𝑑𝐵_𝑖(𝑡), (14) where 𝐿𝑉₁𝑥₁,𝑥₂,𝑥₃=𝑐₁𝑥₁(𝑡)−1𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡−𝜏)+𝑎₁₃𝑥₃(𝑡−𝜏)+𝑐₂𝑥₂(𝑡)−1𝑟₂+𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₂𝑥₂(𝑡)+𝑎₂₃𝑥₃(𝑡−𝜏)+𝑐₃𝑥₃(𝑡)−1𝑟₃+𝑎₃₁𝑥₁(𝑡−𝜏)+𝑎₃₂𝑥₂(𝑡−𝜏)−𝑎₃₃𝑥₃(𝑡+𝑐₁)𝜎²₁+𝑐₂𝜎²₂+𝑐₃𝜎²₃2. (15) By Young’s inequality, we have 𝐿𝑉₁𝑥₁,𝑥₂,𝑥₃=𝑐₁𝑟₁+𝑎₁₁𝑥₁(𝑡)−𝑎₁₂𝑥₂(𝑡−𝜏)−𝑎₁₃𝑥₃(𝑡−𝜏)−𝑎₁₁𝑥²₁(𝑡)+𝑎₁₂𝑥₁(𝑡)𝑥₂(𝑡−𝜏)+𝑎₁₃𝑥₁(𝑡)𝑥₃(𝑡−𝜏)−𝑟₁+𝜎²₁2+𝑐₂𝑟₂+𝑎₂₂𝑥₂(𝑡)−𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₃𝑥₃(𝑡−𝜏)−𝑎₂₂𝑥²₂(𝑡)+𝑎₂₁𝑥₂(𝑡)𝑥₁(𝑡−𝜏)+𝑎₂₃𝑥₂(𝑡)𝑥₃(𝑡−𝜏)−𝑟₂+𝜎²₂2+𝑐₃𝑟₃+𝑎₃₃𝑥₃(𝑡)−𝑎₃₁𝑥₁(𝑡−𝜏)−𝑎₃₂𝑥₂(𝑡−𝜏)−𝑎₃₃𝑥²₃(𝑡)+𝑎₃₁𝑥₃(𝑡)𝑥₁(𝑡−𝜏)+𝑎₃₂𝑥₃(𝑡)𝑥₂(𝑡−𝜏)−𝑟₃+𝜎²₃2≤𝑐₁𝑟₁+𝑎₁₁𝑥₁(𝑡)−𝑎₁₂𝑥₂(𝑡−𝜏)−𝑎₁₃𝑥₃(𝑡−𝜏)+−𝑎₁₁+𝑎₁₂2+𝑎₁₃2𝑥²₁(𝑡)+𝑎₁₂2𝑥²₂()+𝑎𝑡−𝜏¹³2𝑥²₃(𝑡−𝜏)−𝑟₁+𝜎²₁2+𝑐₂𝑟₂+𝑎₂₂𝑥₂(𝑡)−𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₃𝑥₃()+𝑡−𝜏−𝑎₂₂+𝑎₂₁2+𝑎₂₃2𝑥²₂(𝑡)+𝑎₂₁2𝑥²₁()+𝑎𝑡−𝜏²³2𝑥²₃(𝑡−𝜏)−𝑟₂+𝜎²₂2+𝑐₃𝑟₃+𝑎₃₃𝑥₃(𝑡)−𝑎₃₁𝑥₁(𝑡−𝜏)−𝑎₃₂𝑥₂(𝑡−𝜏)+−𝑎₃₃+𝑎₃₁2+𝑎₃₂2𝑥²₃(𝑡)+𝑎₃₁2𝑥²₁()+𝑎𝑡−𝜏³²2𝑥²₂(𝑡−𝜏)−𝑟₃+𝜎²₃2. (16) Secondly, define 𝑉₂𝑥₁,𝑥₂,𝑥₃=_𝑡^𝑡+𝜏𝑎₁₂𝑐₁2+𝑎₃₂𝑐₃2𝑥²₂()+𝑎𝑠−𝜏₁₃𝑐₁2+𝑎₂₃𝑐₂2𝑥²₃()+𝑎𝑠−𝜏₂₁𝑐₂2+𝑎₃₁𝑐₃2𝑥²₁(𝑠−𝜏)𝑑𝑠. (17) By Itô’s formula, we have 𝑑𝑉₂𝑥₁,𝑥₂,𝑥₃=𝑎₁₂𝑐₁2+𝑎₃₂𝑐₃2𝑥²₂(𝑡)+𝑎₁₃𝑐₁2+𝑎₂₃𝑐₂2𝑥²₃(𝑡)+𝑎₂₁𝑐₂2+𝑎₃₁𝑐₃2𝑥²₁(𝑡)−𝑎₁₂𝑐₁2+𝑎₃₂𝑐₃2𝑥²₂()−𝑎𝑡−𝜏₁₃𝑐₁2+𝑎₂₃𝑐₂2𝑥²₃()−𝑎𝑡−𝜏₂₁𝑐₂2+𝑎₃₁𝑐₃2𝑥²₁(𝑡−𝜏)𝑑𝑡. (18) Therefore, let 𝑉𝑥₁,𝑥₂,𝑥₃=𝑉₁𝑥₁,𝑥₂,𝑥₃+𝑉₂𝑥₁,𝑥₂,𝑥₃. (19) Equalities (16) and (18) imply that 𝑑𝑉𝑥₁,𝑥₂,𝑥₃=𝐿𝑉𝑥₁,𝑥₂,𝑥₃+³_𝑖=1𝑐_𝑖𝑥_𝑖(𝑡)−1𝜎_𝑖𝑑𝐵_𝑖(𝑡), (20) where 𝐿𝑉𝑥₁,𝑥₂,𝑥₃≤𝑐₁−𝑥²₁(𝑡)𝑎₁₁−𝑎₁₂−𝑎₁₃+𝑟₁+𝑎₁₁𝑥₁(𝑡)−𝑎₁₂𝑥₂(𝑡−𝜏)−𝑎₁₃𝑥₃(𝑡−𝜏)−𝑟₁+𝜎²₁2+𝑐₂−𝑥²₂(𝑡)𝑎₂₂−𝑎₂₁−𝑎₂₃+𝑟₂+𝑎₂₂𝑥₂(𝑡)−𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₃𝑥₃(𝑡−𝜏)−𝑟₂+𝜎²₂2+𝑐₃−𝑥²₃(𝑡)𝑎₃₃−𝑎₃₁−𝑎₃₂+𝑟₃+𝑎₃₃𝑥₃(𝑡)−𝑎₃₁𝑥₁(𝑡−𝜏)−𝑎₃₂𝑥₂(𝑡−𝜏)−𝑟₃+𝜎²₃2≤𝐾, (21) where 𝐾 is a positive constant. By the method similar to that in [24], the proof is therefore completed.

Under Assumption 1, DDE (8) has a positive equilibrium 𝑥^∗=(𝑥^∗₁,𝑥^∗₂,𝑥^∗₃) which is globally attractive while the system with the stochastic perturbation has a unique global positive solution. It is natural to ask how to estimate the distance between the solution of the deterministic system and the solution of the stochastic system. The following theorem gives us an answer.

Theorem 5. Under Assumption 1, system (9) has the following property: limsup_𝑡→∞1𝑡𝐸^𝑡₀𝑎₁₁−𝑎₁₂−𝑎₁₃𝑐₁+12𝑥₁(𝑠)−𝑥^∗₁²+𝑎₂₂−𝑎₂₁−𝑎₂₃𝑐₂+12𝑥₂(𝑠)−𝑥^∗₂²+𝑎₃₃−𝑎₃₁−𝑎₃₂𝑐₃+12𝑥₃(𝑠)−𝑥^∗₃²≤∑𝑑𝑠_3𝑖=1𝑐_𝑖𝑥^∗_𝑖𝜎²_𝑖2, (22) where 𝑥(𝑡) is a solution on 𝑡≥𝜏 to (9) with an initial value 𝑥(⋅)∈𝐶([−𝜏,0];𝑅³₊), 𝑐_𝑖(𝑖=1,2,3) is defined in the proof of Theorem 4, and 𝑥^∗=(𝑥^∗₁,𝑥^∗₂,𝑥^∗₃) is the unique positive equilibrium of system (8); namely, 𝑥^∗₁=𝑎₂₂𝑎₃₃−𝑎₂₃𝑎₃₂𝑟₁+𝑎₁₂𝑎₃₃+𝑎₁₃𝑎₃₂𝑟₂+𝑎₁₂𝑎₂₃+𝑎₁₃𝑎₂₂𝑟₃×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹,𝑥^∗₂=𝑎₂₁𝑎₃₃+𝑎₂₃𝑎₃₁𝑟₁+𝑎₁₁𝑎₃₃−𝑎₁₃𝑎₃₁𝑟₂+𝑎₁₁𝑎₂₃+𝑎₁₃𝑎₂₁𝑟₃×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹,𝑥^∗₃=𝑎₂₁𝑎₃₂+𝑎₂₂𝑎₃₁𝑟₁+𝑎₁₁𝑎₃₂+𝑎₁₂𝑎₃₁𝑟₂+𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁𝑟₃×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹. (23)

For brevity, we will give the proof in the appendix.

3. Persistence in Time Average

For convenience, we denote the unique global solution of system (9) by 𝑥(𝑡,𝜉) with an initial data 𝜉∈𝐶([−𝜏,0];𝑅⁺₃). Theorem 4 shows that the solution of system (9) will remain positive under Assumption 1. This property gives us an opportunity to investigate how the solution varies in 𝑅⁺₃. In population dynamics, one of the most attractive properties is persistence which means all species will coexist. Now we give the definition of persistence in time average.

Definition 6. System (9) is persistent in time average, if, for any initial data 𝜉∈𝐶([−𝜏,0];𝑅³₊), the solution 𝑥(𝑡,𝜉) has the property that 0<liminf_𝑡→∞1𝑡^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠≤limsup_𝑡→∞1𝑡^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠<+∞a.s.𝑖=1,2,3. (24)

To prove that system (9) is persistent in time average, we will cite a lemma. Jiang and Shi in [17] discussed a stochastic nonautonomous logistic equation 𝑑𝑁(𝑡)=𝑁(𝑡)[(𝑎(𝑡)−𝑏(𝑡)𝑁(𝑡))𝑑𝑡+𝛼(𝑡)𝑑𝐵(𝑡)], (25) where 𝐵(𝑡) is 1-dimensional standard Brownian motion, 𝑁(0)=𝑁₀, and 𝑁₀ is independent of 𝐵(𝑡). They obtained the following result.

Lemma 7 (see [17]). Assume that 𝑎(𝑡),𝑏(𝑡), and 𝛼(𝑡) are bounded continuous functions defined on [0,∞), 𝑎(𝑡)>0, and 𝑏(𝑡)>0. Then there exists a unique continuous positive solution of (25) for any initial value 𝑁(0)=𝑁₀>0, which is global and represented by 𝑁(𝑡)=exp^𝑡₀𝑎(𝑠)−𝛼²(𝑠)2(𝑠)(𝑠×1𝑑𝑠+𝛼𝑑𝐵)𝑁₀+^𝑡₀𝑏(𝑠)exp^𝑠₀𝑎(𝜏)−𝛼²(𝜏)2(𝜏)(𝜏𝑑𝜏+𝛼𝑑𝐵)𝑑𝑠⁻¹,𝑡≥0. (26) Moreover, the solution 𝑁(𝑡) has the property that lim_𝑡→∞log𝑁(𝑡)𝑡=0,𝑎.𝑠. (27)

Remark 8. Let 𝜙_𝑖(𝑡), 𝑖=1,2,3, be the solution of the following equation: 𝑑𝜙_𝑖(𝑡)=𝜙_𝑖(𝑡)𝑟_𝑖−𝑎_𝑖𝑖𝜙_𝑖(𝑡)𝑑𝑡+𝜎_𝑖𝑑𝐵_𝑖(𝑡),𝑡≥0,𝑖=1,2,3, (28) with initial value 𝜙_𝑖(0)=𝑥_𝑖(0). From Lemma 7, we have 𝜙_𝑖(𝑡)=𝑒^{(𝑟𝑖−(𝜎2𝑖/2))𝑡+𝜎𝑖𝐵𝑖(𝑡)}1/𝑥_𝑖(0)+𝑎_𝑖𝑖∫^𝑡₀𝑒^{(𝑟𝑖−(𝜎2𝑖/2))𝑠+𝜎𝑖𝐵𝑖(𝑠)}𝑑𝑠,𝑡≥0. (29) From the result in [13], we know lim_𝑡→∞1𝑡^𝑡₀𝜙_𝑖(𝑠)𝑟𝑑𝑠=^𝑖−𝜎²_𝑖/2𝑎_𝑖𝑖,lim_𝑡→∞log𝜙_𝑖(𝑡)𝑡=0,a.s.𝑖=1,2,3 (30) provided Assumption 2.

Theorem 9. Under Assumptions 1 and 2, system (9) is persistent in time average.

Proof. From Lemma 7, we know 𝑥_𝑖(𝑡)=𝑒^{∫𝑡0(𝑟𝑖+∑3𝑗=1,𝑗≠𝑖𝑎𝑖𝑗𝑥𝑗(𝑠−𝜏)−(𝜎2𝑖/2))𝑑𝑠+𝜎𝑖𝐵𝑖(𝑡)}1/𝑥_𝑖(0)+𝑎_𝑖𝑖∫^𝑡₀𝑒^{∫𝑠0(𝑟𝑖+∑3𝑗=1,𝑗≠𝑖𝑎𝑖𝑗𝑥𝑗(𝑢−𝜏)−(𝜎2𝑖/2))𝑑𝑢+𝜎𝑖𝐵𝑖(𝑠)}=(1)𝑒𝑑𝑠_{−𝜎𝑖𝐵𝑖(𝑡)}×1𝑥_𝑖(0)𝑒^{−∫𝑡0(𝑟𝑖+∑3𝑗=1,𝑗≠𝑖𝑎𝑖𝑗𝑥𝑗(𝑠−𝜏)−(𝜎2𝑖/2))𝑑𝑠}+𝑎_𝑖𝑖×^𝑡₀𝑒^{−∫𝑡𝑠(𝑟𝑖+∑3𝑗=1,𝑗≠𝑖𝑎𝑖𝑗𝑥𝑗(𝑢−𝜏)−(𝜎2𝑖/2))𝑑𝑢+𝜎𝑖𝐵𝑖(𝑠)}𝑑𝑠⁻¹. (31) From Remark 8, yields 𝑥_𝑖(𝑡)≥𝜙_𝑖(𝑡),𝑖=1,2,3. (32) Together with Lemma 7, it is easy to obtain 𝑟_𝑖−𝜎²_𝑖/2𝑎_𝑖𝑖≤liminf_𝑡→∞1𝑡^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠,liminf_𝑡→∞log𝑥_𝑖(𝑡)𝑡≥0a.s.𝑖=1,2,3. (33) The inequality limsup_𝑡→∞(1/𝑡)∫^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠<+∞ a.s. 𝑖=1,2,3 will be shown in Theorem 11.

Theorem 9 shows that system (9) is persistent in time average if the intensity of noise is small. Next we want to obtain the limit of the solution in time average of system (9). We begin from the lemma in [29].

Lemma 10. Let 𝑓∈𝐶[[0,∞)×Ω,(0,∞)], 𝐹(𝑡)∈((0,∞)×Ω,𝑅). If there exist positive constants 𝜆₀ and 𝜆 such that log𝑓(𝑡)≥𝜆𝑡−𝜆₀^𝑡₀𝑓(𝑠)𝑑𝑠+𝐹(𝑡),𝑡≥0,𝑎.𝑠., (34) and lim_𝑡→∞(𝐹(𝑡)/𝑡)=0 a.s., then liminf_𝑡→∞1𝑡^𝑡₀𝑓(𝑠)𝜆𝑑𝑠≥𝜆₀,𝑎.𝑠. (35)

Theorem 11. Under Assumptions 1 and 2, for any initial data 𝜉∈𝐶([−𝜏,0];𝑅³₊), the solution 𝑥(𝑡,𝜉) has the property that lim_𝑡→∞1𝑡^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠=𝑥^∗_𝑖𝑎.𝑠.𝑖=1,2,3, (36) where 𝑥^∗₁=𝑎₂₂𝑎₃₃−𝑎₂₃𝑎₃₂𝑟₁−𝜎²₁2+𝑎₁₂𝑎₃₃+𝑎₁₃𝑎₃₂×𝑟₂−𝜎²₂2+𝑎₁₂𝑎₂₃+𝑎₁₃𝑎₂₂𝑟₃−𝜎²₃2×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹,𝑥^∗₂=𝑎₂₁𝑎₃₃+𝑎₂₃𝑎₃₁𝑟₁−𝜎²₁2+𝑎₁₁𝑎₃₃−𝑎₁₃𝑎₃₁×𝑟₂−𝜎²₂2+𝑎₁₁𝑎₂₃+𝑎₁₃𝑎₂₁𝑟₃−𝜎²₃2×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹,𝑥^∗₃=𝑎₂₁𝑎₃₂+𝑎₂₂𝑎₃₁𝑟₁−𝜎²₁2+𝑎₁₁𝑎₃₂+𝑎₁₂𝑎₃₁×𝑟₂−𝜎²₂2+𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁𝑟₃−𝜎²₃2×𝑎₁₁𝑎₂₂𝑎₃₃−𝑎₁₁𝑎₂₃𝑎₃₂−𝑎₂₂𝑎₁₃𝑎₃₁−𝑎₃₃𝑎₁₂𝑎₂₁−𝑎₁₂𝑎₂₃𝑎₃₁−𝑎₂₁𝑎₃₂𝑎₁₃⁻¹. (37)

The mathematical derivations are lengthy; we will give the proof in the appendix.

4. Nonpersistence

In this section, we will show that the system (9) is nonpersistent if the intensity of the noise is big enough; however, it does not occur to the deterministic system. First of all, we give the definition of nonpersistence.

Definition 12. System (9) is nonpersistent, if there are positive constants 𝑏₁,𝑏₂,𝑏₃ such that lim³_{𝑡→∞𝑖=1}𝑥^𝑏𝑖_𝑖(𝑡)=0a.s. (38)

Theorem 13. Under Assumption 1, if  𝐶<0 holds, system (9) is nonpersistent, where 𝐶=𝑏₁₁(𝑟₁−(𝜎²₁/2))+𝑏₁₂(𝑟₂−(𝜎²₂/2))+𝑏₁₃(𝑟₃−(𝜎²₃/2)),𝑏₁₁>0; 𝑏₁₂,𝑏₁₃ is defined by (B.24).

Proof. It follows from (B.25) that 𝑏₁₁log𝑥₁(𝑡)+𝑏₁₂log𝑥₂(𝑡)+𝑏₁₃log𝑥₃(𝑡)≤𝑏₁₁𝑟₁+𝑏₁₂𝑟₂+𝑏₁₃𝑟₃𝑡−𝑚₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+³_𝑖=1𝑏_1𝑖𝜎_𝑖𝐵_𝑖(𝑡)+³_𝑖=1𝑏_1𝑖𝑈_𝑖≤𝑏₁₁𝑟₁+𝑏₁₂𝑟₂+𝑏₁₃𝑟₃𝑡+³_𝑖=1𝑏_1𝑖𝜎_𝑖𝐵_𝑖(𝑡)+³_𝑖=1𝑏_1𝑖𝑈_𝑖. (39) Together with lim_𝑡→∞(𝐵_𝑖(𝑡)/𝑡)=0andlim_𝑡→∞(𝑈_𝑖/𝑡)=0,𝑖=1,2,3, we have limsup_𝑡→∞1𝑡³_𝑖=1𝑏_1𝑖log𝑥_𝑖(𝑡)≤𝑏₁₁𝑟₁+𝑏₁₂𝑟₂+𝑏₁₃𝑟₃. (40) If 𝑏₁₁𝑟₁+𝑏₁₂𝑟₂+𝑏₁₃𝑟₃<0 holds, it follows lim³_{𝑡→∞𝑖=1}𝑥^𝑏1𝑖_𝑖(𝑡)=0a.s. (41) Hence, system (9) is nonpersistent. The proof is completed.

By a similar method, we can yield the following theorems.

Theorem 14. Under Assumption 1, if  𝐶<0 holds, system (9) is nonpersistent, where 𝐶=𝑏₂₁(𝑟₁−(𝜎²₁/2))+𝑏₂₂(𝑟₂−(𝜎²₂/2))+𝑏₂₃(𝑟₃−(𝜎²₃/2)), 𝑏₂₂>0, 𝑏₂₁=((−𝑎₂₁𝑎₃₃−𝑎₂₃𝑎₃₁)/(𝑎₁₃𝑎₃₁−𝑎₁₁𝑎₃₃))𝑏₂₂,𝑏₂₃=((−𝑎₁₃𝑎₂₁−𝑎₁₁𝑎₂₃)/(𝑎₁₃𝑎₃₁−𝑎₁₁𝑎₃₃))𝑏₂₂.

Theorem 15. Under Assumption 1, if  𝐶<0 holds, system (9) is nonpersistent, where 𝐶=𝑏₃₁(𝑟₁−(𝜎²₁/2))+𝑏₃₂(𝑟₂−(𝜎²₂/2))+𝑏₂₃(𝑟₃−(𝜎²₃/2)),𝑏₃₃>0,𝑏₃₁=((−𝑎₃₁𝑎₂₂−𝑎₃₂𝑎₂₁)/(𝑎₁₂𝑎₂₁−𝑎₁₁𝑎₂₂))𝑏₃₃,𝑏₃₂=((−𝑎₃₁𝑎₁₂−𝑎₃₂𝑎₁₁)/(𝑎₁₂𝑎₂₁−𝑎₁₁𝑎₂₂))𝑏₃₃.

5. Numerical Examinations

In this section, we give the numerical examinations to illustrate above results. By the method mentioned in [30], consider the discrete equation: 𝑥_1,𝑘+1=𝑥_1,𝑘+𝑥_1,𝑘𝑟₁−𝑎₁₁𝑥_1,𝑘+𝑎₁₂𝑥_{2,𝑘−𝑚}+𝑎₁₃𝑥_{3,𝑘−𝑚}Δ𝑡+𝜎₁𝜀_1,𝑘√1Δ𝑡+2𝜎²₁𝜀²_1,𝑘,𝑥Δ𝑡−Δ𝑡_2,𝑘+1=𝑥_2,𝑘+𝑥_2,𝑘𝑟₂+𝑎₂₁𝑥_{1,𝑘−𝑚}−𝑎₂₂𝑥_2,𝑘+𝑎₂₃𝑥_{3,𝑘−𝑚}Δ𝑡+𝜎₂𝜀_2,𝑘√1Δ𝑡+2𝜎²₂𝜀²_2,𝑘,𝑥Δ𝑡−Δ𝑡_3,𝑘+1=𝑥_3,𝑘+𝑥_3,𝑘𝑟₃+𝑎₃₁𝑥_{1,𝑘−𝑚}+𝑎₃₂𝑥_{2,𝑘−𝑚}−𝑎₃₃𝑥_3,𝑘Δ𝑡+𝜎₃𝜀_3,𝑘√1Δ𝑡+2𝜎²₃𝜀²_3,𝑘Δ𝑡−Δ𝑡, (42) where 𝑚 represents the integer part of 𝜏/Δ𝑡−1. Choosing suitable parameters in the system, by Matlab we get the simulation figures with initial value (𝑥₁(𝑡),𝑥₂(𝑡),𝑥₃(𝑡))≡(0.7,0.3,0.9),𝑡∈[−𝜏,0]. (For convenience we let the initial value be a constant function; otherwise we have to give 𝑚+1 values.) The time step Δ𝑡=0.005; we always choose 𝑎₁₁𝑎₁₂𝑎₁₃𝑎₂₁𝑎₂₂𝑎₂₃𝑎₃₁𝑎₃₂𝑎₃₃=,𝑟0.70.20.40.10.50.30.20.40.8₁𝑟₂𝑟₃=.0.40.20.3 (43) Then 𝑎₁₁−𝑎₁₂−𝑎₁₃=0.1,𝑎₂₂−𝑎₂₁−𝑎₂₃=0.2,𝑎₃₃−𝑎₃₁−𝑎₃₂=0.2,(𝑥^∗₁,𝑥^∗₂,𝑥^∗₃)=(2.2679,2.0268,1.9554). By choosing different intensities of the noise and time delays, we obtain the following cases to illustrate our results.

Case 1 (persistence). Choosing 𝜏=20.2𝜎₁=0.04,𝜎₂=0.08,𝜎₃=0.07, then we have 𝑟_𝑖>𝜎²_𝑖/2,𝑖=1,2,3. Hence all assumptions of Theorem 11 are satisfied. Figure 1 shows that the solution fluctuates in a small zone.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

The solution (imaginary line) of system (8) and the solution (real line) of system (9) with 𝜎1=0.04,𝜎2=0.08,𝜎3=0.07,Δ𝑡=0.005.

Case 2 (nonpersistence (we only illustrate the first situation)). (1) 𝑥₁(𝑡) Is Disturbed by a Big Noise, Which Leads to Nonpersistence. Choosing 𝜏=20.2,𝜎₁=1.53,𝜎₂=0.08,𝜎₃=0.07,𝑏₁₁=0.10,𝑏₁₂=0.1143,𝑏₁₃=0.0929, then we have 𝑟_𝑖>𝜎²_𝑖,𝑖=2,3 but 𝑏₁₁(𝑟₁−(𝜎²₁/2))+𝑏₁₂(𝑟₂−(𝜎²₂/2))+𝑏₁₃(𝑟₃−(𝜎²₃/2))<0. For convenience we shorten 𝑏_1𝑖 to 𝑏_𝑖. Hence the assumptions of Theorem 13 are satisfied, and system (9) is nonpersistent (see Figure 2(d)). In addition, since 𝑟_𝑖>𝜎²_𝑖,𝑖=2,3, 𝑥₂(𝑡),𝑥₃(𝑡) do not tend to zero in time average by Theorem 9. So lim_𝑡→∞𝑥₁(𝑡)=0 a.s. (see Figures 2(a)–2(c)).
(2) The Second Specie 𝑥₂(𝑡) Is Disturbed by the Big White Noise, Which Leads to the Nonpersistence. Choosing 𝜏=20.2,𝜎₁=0.04,𝜎₂=1.17,𝜎₃=0.07,𝑏₁=0.10,𝑏₂=0.1143,𝑏₃=0.0929, then we have 𝑟_𝑖>𝜎²_𝑖,𝑖=1,3, but 𝑏₁(𝑟₁−(𝜎²₁/2))+𝑏₂(𝑟₂−(𝜎²₂/2))+𝑏₃(𝑟₃−(𝜎²₃/2))<0. Hence the assumptions of Theorem 13 are satisfied, and system (9) is nonpersistent (see Figure 3(d)). In addition, since 𝑟_𝑖>𝜎²_𝑖,𝑖=1,3,𝑥₁(𝑡),𝑥₃(𝑡) do not tend to zero in time average by Theorem 9. So lim_𝑡→∞𝑥₂(𝑡)=0 a.s. (see Figures 3(a)–3(c)).
(3) The Third Specie 𝑥₃(𝑡) Is Disturbed by the Big White Noise, Which Leads to the Nonpersistence. Choosing 𝜏=20.2,𝜎₁=0.04,𝜎₂=0.08,𝜎₃=1.36,𝑏₁=0.10,𝑏₂=0.1143,𝑏₃=0.0929, then we have 𝑟_𝑖>𝜎²_𝑖,𝑖=1,2, but 𝑏₁(𝑟₁−(𝜎²₁/2))+𝑏₂(𝑟₂−(𝜎²₂/2))+𝑏₃(𝑟₃−(𝜎²₃/2))<0. Hence the assumptions of Theorem 13 are satisfied, and system (9) is nonpersistent (see Figure 4(d)). In addition, since 𝑟_𝑖>𝜎²_𝑖,𝑖=1,2, 𝑥₁(𝑡),𝑥₂(𝑡) do not tend to zero in time average by Theorem 9. So we have lim_𝑡→∞𝑥₃(𝑡)=0 a.s. (see Figures 4(a)–4(c)).
(4) The First Two Species 𝑥₁(𝑡),𝑥₂(𝑡) Are Disturbed by the Big White Noises, Which Leads to the Nonpersistence. Choosing 𝜏=20.2,𝜎₁=1.25,𝜎₂=1.01,𝜎₃=0.07,𝑏₁=0.10,𝑏₂=0.1143,𝑏₃=0.0929, then we have 𝑟₃>𝜎²₃, but 𝑏₁(𝑟₁−(𝜎²₁/2))+𝑏₂(𝑟₂−(𝜎²₂/2))+𝑏₃(𝑟₃−(𝜎²₃/2))<0. Hence the assumptions of Theorem 13 are satisfied; then system (9) is nonpersistent (see Figure 5(d)). In addition, since 𝑟₃>𝜎²₃, 𝑥₃(𝑡) does not tend to zero in time average by Theorem 9. So we have lim_𝑡→∞𝑥_𝑖(𝑡)=0,𝑖=1,2 a.s. (see Figures 5(a)–5(c)).
(5) The First and the Third Species 𝑥₁(𝑡), 𝑥₃(𝑡) Are Disturbed by the Big White Noises, Which Leads to the Nonpersistence. Choosing 𝜏=20.2, 𝜎₁=1.25, 𝜎₂=0.08,𝜎₃=1.12,𝑏₁=0.10, 𝑏₂=0.1143, 𝑏₃=0.0929, then we have 𝑟₂>𝜎²₂, but 𝑏₁(𝑟₁−(𝜎²₁/2))+𝑏₂(𝑟₂−(𝜎²₂/2))+𝑏₃(𝑟₃−(𝜎²₃/2))<0. Hence the assumptions of Theorem 13 are satisfied, and system (9) is nonpersistent (see Figure 6(d)). In addition, since 𝑟₂>𝜎²₂, 𝑥₂(𝑡) does not tend to zero in time average by Theorem 9. So we have lim_𝑡→∞𝑥_𝑖(𝑡)=0,𝑖=1,3 a.s. (see Figures 6(a)–6(c)).
(6) The Last Two Species 𝑥₂(𝑡), 𝑥₃(𝑡) Are Disturbed by the Big White Noises, Which Leads to Nonpersistence. Choosing 𝜏=20.2,𝜎₁=0.04,𝜎₂=1.17,𝜎₃=1.36,𝑏₁=0.10,𝑏₂=0.1143,𝑏₃=0.0929, then we have 𝑟₁>𝜎²₁, but 𝑏₁(𝑟₁−(𝜎²₁/2))+𝑏₂(𝑟₂−(𝜎²₂/2))+𝑏₃(𝑟₃−(𝜎²₃/2))<0. Hence the assumptions of Theorem 13 are satisfied, and system (9) is nonpersistent (see Figure 7(d)). In addition, since 𝑟₁>𝜎²₁, 𝑥₁(𝑡) does not tend to zero in time average by Theorem 9. So we have lim_𝑡→∞𝑥_𝑖(𝑡)=0,𝑖=2,3 a.s. (see Figures 7(a)–7(c)).
(7) All the Three Species 𝑥₁(𝑡), 𝑥₂(𝑡), 𝑥₃(𝑡) Are Disturbed by the Big White Noise, Which Leads to the Nonpersistence. Choosing 𝜏=20.2,𝜎₁=1.25,𝜎₁=1.01,𝜎₁=1.36,𝑏₁=0.10,𝑏₂=0.1143,𝑏₃=0.0929, then we have 𝑏₁(𝑟₁−(𝜎²₁/2))+𝑏₂(𝑟₂−(𝜎²₂/2))+𝑏₃(𝑟₃−(𝜎²₃/2))<0. Hence the assumptions of Theorem 13 are satisfied, and system (9) is nonpersistent (see Figure 8(d)). So we have lim_𝑡→∞𝑥_𝑖(𝑡)=0,𝑖=1,2,3 a.s. (see Figures 8(a)–8(c)).

(a)

(b)

(c)

(d)

Figure 2 

The solution (imaginary line) of system (8) and the solution (real line) of system (9) with 𝜎1=1.53,𝜎2=0.08,𝜎3=0.07,Δ𝑡=0.005.

(a)

(b)

(c)

(d)

Figure 3 

The solution (imaginary line) of system (8) and the solution (real line) of system (9) with 𝜎1=0.04,𝜎2=1.17,𝜎3=0.07,Δ𝑡=0.005.

(a)

(b)

(c)

(d)

Figure 4 

The solution (imaginary line) of system (8) and the solution (real line) of system (9) with 𝜎1=0.04,𝜎2=0.08,𝜎3=1.36,Δ𝑡=0.005.

(a)

(b)

(c)

(d)

Figure 5 

The solution (imaginary line) of system (8) and the solution (real line) of system (9) with 𝜎1=1.25,𝜎2=1.01,𝜎3=0.07,Δ𝑡=0.005.

(a)

(b)

(c)

(d)

Figure 6 

The solution (imaginary line) of system (8) and the solution (real line) of system (9) with 𝜎1=1.25,𝜎2=0.08,𝜎3=1.12,Δ𝑡=0.005.

(a)

(b)

(c)

(d)

Figure 7 

The solution (imaginary line) of system (8) and the solution (real line) of system (9) with 𝜎1=0.04,𝜎2=1.17,𝜎3=1.36,Δ𝑡=0.005.

(a)

(b)

(c)

(d)

Figure 8 

The solution (imaginary line) of system (8) and the solution (real line) of system (9) with 𝜎1=1.25,𝜎2=1.01,𝜎3=1.36,Δ𝑡=0.005.

Appendices

A. Proof of Theorem 5

Proof. Define a 𝐶² function 𝑉₁∶𝑅³₊→𝑅₊ by 𝑉₁𝑥₁,𝑥₂,𝑥₃=𝑐₁𝑥₁(𝑡)−𝑥^∗₁−𝑥^∗₁𝑥log¹(𝑡)𝑥^∗₁+𝑐₂𝑥₂(𝑡)−𝑥^∗₂−𝑥^∗₂𝑥log²(𝑡)𝑥^∗₂+𝑐₃𝑥₃(𝑡)−𝑥^∗₃−𝑥^∗₃𝑥log³(𝑡)𝑥^∗₃. (A.1) By Itô’s formula, we have 𝑑𝑉₁𝑥₁,𝑥₂,𝑥₃=𝑐₁𝑥1−^∗1𝑥₁(𝑡)𝑑𝑥₁(𝑡)+𝑐₂𝑥1−^∗2𝑥₂(𝑡)𝑑𝑥₂(𝑡)+𝑐₃𝑥1−^∗3𝑥₃(𝑡)𝑑𝑥₃(𝑡)+𝑐₁2𝑥^∗₁𝜎²₁+𝑐𝑑𝑡₂2𝑥^∗₂𝜎²₂𝑐𝑑𝑡+³2𝑥^∗₃𝜎²₃𝑑𝑡=𝑐₁𝑥₁(𝑡)−𝑥^∗₁𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡−𝜏)+𝑎₁₃𝑥₃(𝑡−𝜏)𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡)+𝑐₂𝑥₂(𝑡)−𝑥^∗₂𝑟₂+𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₂𝑥₂(𝑡)+𝑎₂₃𝑥₃(𝑡−𝜏)𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡)+𝑐₃𝑥₃(𝑡)−𝑥^∗₃𝑟₃+𝑎₃₁𝑥₁(𝑡−𝜏)+𝑎₃₂𝑥₂(𝑡−𝜏)−𝑎₃₃𝑥₃(𝑡)𝑑𝑡+𝜎₃𝑑𝐵₃(𝑡+𝑐₁)2𝑥^∗₁𝜎²₁+𝑐₂2𝑥^∗₂𝜎²₂+𝑐₃2𝑥^∗₃𝜎²₃𝑑𝑡, (A.2) where 𝐿𝑉₁𝑥₁,𝑥₂,𝑥₃=𝑐₁𝑥₁(𝑡)−𝑥^∗₁𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡−𝜏)+𝑎₁₃𝑥₃(𝑡−𝜏)+𝑐₂𝑥₂(𝑡)−𝑥^∗₂𝑟₂+𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₂𝑥₂(𝑡)+𝑎₂₃𝑥₃(𝑡−𝜏)+𝑐₃𝑥₃(𝑡)−𝑥^∗₃𝑟₃+𝑎₃₁𝑥₁(𝑡−𝜏)+𝑎₃₂𝑥₂(𝑡−𝜏)−𝑎₃₃𝑥₃(𝑡+𝑐₁)𝑥^∗₁𝜎²₁+𝑐₂𝑥^∗₂𝜎²₂+𝑐₃𝑥^∗₃𝜎²₃2. (A.3) Since 𝑥^∗=(𝑥^∗₁,𝑥^∗₂,𝑥^∗₃) is the equilibrium point of system (8), we have 𝐿𝑉₁𝑥₁,𝑥₂,𝑥₃=𝑐₁𝑥₁(𝑡)−𝑥^∗₁−𝑎₁₁𝑥₁(𝑡)−𝑥^∗₁+𝑎₁₂𝑥₂(𝑡−𝜏)−𝑥^∗₂+𝑎₁₃𝑥₃(𝑡−𝜏)−𝑥^∗₃+𝑐₂𝑥₂(𝑡)−𝑥^∗₂𝑎₂₁𝑥₁(𝑡−𝜏)−𝑥^∗₁−𝑎₂₂𝑥₂(𝑡)−𝑥^∗₂+𝑎₂₃𝑥₃(𝑡−𝜏)−𝑥^∗₃+𝑐₃𝑥₃(𝑡)−𝑥^∗₃𝑎₃₁𝑥₁(𝑡−𝜏)−𝑥^∗₁+𝑎₃₂𝑥₂(𝑡−𝜏)−𝑥^∗₂−𝑎₃₃𝑥₃(𝑡)−𝑥^∗₃+𝑐₁𝑥^∗₁𝜎²₁+𝑐₂𝑥^∗₂𝜎²₂+𝑐₃𝑥^∗₃𝜎²₃2=−𝑎₁₁𝑐₁𝑥₁(𝑡)−𝑥^∗₁²+𝑎₁₂𝑐₁𝑥₁(𝑡)−𝑥^∗₁𝑥₂(𝑡−𝜏)−𝑥^∗₂+𝑎₁₃𝑐₁𝑥₁(𝑡)−𝑥^∗₁𝑥₃(𝑡−𝜏)−𝑥^∗₃+𝑎₂₁𝑐₂𝑥₂(𝑡)−𝑥^∗₂𝑥₁(𝑡−𝜏)−𝑥^∗₁−𝑎₂₂𝑐₂𝑥₂(𝑡)−𝑥^∗₂²+𝑎₂₃𝑐₂𝑥₂(𝑡)−𝑥^∗₂𝑥₃(𝑡−𝜏)−𝑥^∗₃+𝑎₃₁𝑐₃𝑥₃(𝑡)−𝑥^∗₃𝑥₁(𝑡−𝜏)−𝑥^∗₁+𝑎₃₂𝑐₃𝑥₃(𝑡)−𝑥^∗₃𝑥₂(𝑡−𝜏)−𝑥^∗₂−𝑎₃₃𝑐₃𝑥₃(𝑡)−𝑥^∗₃²+𝑐₁𝑥^∗₁𝜎²₁+𝑐₂𝑥^∗₂𝜎²₂+𝑐₃𝑥^∗₃𝜎²₃2. (A.4) By Young inequality, we have 𝑎₁₂𝑐₁𝑥₁(𝑡)−𝑥^∗₁𝑥₂(𝑡−𝜏)−𝑥^∗₂≤𝑎₁₂𝑐₁2𝑥₁(𝑡)−𝑥^∗₁²+𝑥₂(𝑡−𝜏)−𝑥^∗₂²,𝑎₁₃𝑐₁𝑥₁(𝑡)−𝑥^∗₁𝑥₃(𝑡−𝜏)−𝑥^∗₃≤𝑎₁₃𝑐₁2𝑥₁(𝑡)−𝑥^∗₁²+𝑥₃(𝑡−𝜏)−𝑥^∗₃²,𝑎₂₁𝑐₂𝑥₂(𝑡)−𝑥^∗₂𝑥₁(𝑡−𝜏)−𝑥^∗₁≤𝑎₂₁𝑐₂2𝑥₂(𝑡)−𝑥^∗₂²+𝑥₁(𝑡−𝜏)−𝑥^∗₁²,𝑎₂₃𝑐₂𝑥₂(𝑡)−𝑥^∗₂𝑥₃(𝑡−𝜏)−𝑥^∗₃≤𝑎₂₃𝑐₂2𝑥₂(𝑡)−𝑥^∗₂²+𝑥₃(𝑡−𝜏)−𝑥^∗₃²,𝑎₃₁𝑐₃𝑥₃(𝑡)−𝑥^∗₃𝑥₁(𝑡−𝜏)−𝑥^∗₁≤𝑎₂₃𝑐₃2𝑥₃(𝑡)−𝑥^∗₃²+𝑥₁(𝑡−𝜏)−𝑥^∗₁²,𝑎₃₂𝑐₃𝑥₃(𝑡)−𝑥^∗₃𝑥₂(𝑡−𝜏)−𝑥^∗₂≤𝑎₃₂𝑐₃2𝑥₃(𝑡)−𝑥^∗₃²+𝑥₂(𝑡−𝜏)−𝑥^∗₂². (A.5) Substituting (A.5) into (A.4), we yield 𝐿𝑉₁𝑥₁,𝑥₂,𝑥₃≤−𝑎₁₁𝑐₁+𝑎₁₂𝑐₁2+𝑎₁₃𝑐₁2𝑥₁(𝑡)−𝑥^∗₁²+−𝑎₂₂𝑐₂+𝑎₂₁𝑐₂2+𝑎₂₃𝑐₂2𝑥₂(𝑡)−𝑥^∗₂²+−𝑎₃₃𝑐₃+𝑎₃₁𝑐₃2+𝑎₃₂𝑐₃2𝑥₃(𝑡)−𝑥^∗₃²+𝑎₂₁𝑐₂2+𝑎₃₁𝑐₃2𝑥₁(𝑡−𝜏)−𝑥^∗₁²+𝑎₁₂𝑐₁2+𝑎₃₂𝑐₃2𝑥₂(𝑡−𝜏)−𝑥^∗₂²+𝑎₁₃𝑐₁2+𝑎₂₃𝑐₂2𝑥₃(𝑡−𝜏)−𝑥^∗₃²+𝑐₁𝑥^∗₁𝜎²₁+𝑐₂𝑥^∗₂𝜎²₂+𝑐₃𝑥^∗₃𝜎²₃2. (A.6) From (A.6), we have 𝑑𝑉₁𝑥₁,𝑥₂,𝑥₃≤−𝑎₁₁𝑐₁+𝑎₁₂𝑐₁2+𝑎₁₃𝑐₁2𝑥₁(𝑡)−𝑥^∗₁²+−𝑎₂₂𝑐₂+𝑎₂₁𝑐₂2+𝑎₂₃𝑐₂2𝑥₂(𝑡)−𝑥^∗₂²+−𝑎₃₃𝑐₃+𝑎₃₁𝑐₃2+𝑎₃₂𝑐₃2𝑥₃(𝑡)−𝑥^∗₃²+𝑎₂₁𝑐₂2+𝑎₃₁𝑐₃2𝑥₁(𝑡−𝜏)−𝑥^∗₁²+𝑎₁₂𝑐₁2+𝑎₃₂𝑐₃2𝑥₂(𝑡−𝜏)−𝑥^∗₂²+𝑎₁₃𝑐₁2+𝑎₂₃𝑐₂2𝑥₃(𝑡−𝜏)−𝑥^∗₃²+𝑐₁𝑥^∗₁𝜎²₁+𝑐₂𝑥^∗₂𝜎²₂+𝑐₃𝑥^∗₃𝜎²₃2𝑑𝑡+³_𝑖=1𝜎_𝑖𝑑𝐵_𝑖. (A.7) Define 𝑉₂𝑥₁,𝑥₂,𝑥₃=_𝑡^𝑡+𝜏𝑎₁₂𝑐₁2+𝑎₃₂𝑐₃2𝑥₂(𝑠−𝜏)−𝑥^∗₂²+𝑎₁₃𝑐₁2+𝑎₂₃𝑐₂2𝑥₃(𝑠−𝜏)−𝑥^∗₃²+𝑎₂₁𝑐₂2+𝑎₃₁𝑐₃2𝑥₁(𝑠−𝜏)−𝑥^∗₁²𝑑𝑠. (A.8) Then by Itô’s formula, we have 𝑑𝑉₂𝑥₁,𝑥₂,𝑥₃=𝑎₁₂𝑐₁2+𝑎₃₂𝑐₃2𝑥₂(𝑡)−𝑥^∗₂²+𝑎₁₃𝑐₁2+𝑎₂₃𝑐₂2𝑥₃(𝑡)−𝑥^∗₃²+𝑎₂₁𝑐₂2+𝑎₃₁𝑐₃2+𝑥₁(𝑡)−𝑥^∗₁²−𝑎₁₂𝑐₁2+𝑎₃₂𝑐₃2𝑥₂(𝑡−𝜏)−𝑥^∗₂²−𝑎₁₃𝑐₁2+𝑎₂₃𝑐₂2𝑥₃(𝑡−𝜏)−𝑥^∗₃²−𝑎₂₁𝑐₂2+𝑎₃₁𝑐₃2𝑥₁(𝑡−𝜏)−𝑥^∗₁²𝑑𝑡. (A.9) Define 𝑉𝑥₁,𝑥₂,𝑥₃=𝑉₁𝑥₁,𝑥₂,𝑥₃+𝑉₂𝑥₁,𝑥₂,𝑥₃. (A.10) Together with (A.7) and (A.9), it implies 𝑑𝑉𝑥₁,𝑥₂,𝑥₃≤𝑎₁₂+𝑎₁₃−𝑎₁₁𝑐₁−12𝑥₁(𝑡)−𝑥^∗₁²+𝑎₂₁+𝑎₂₃−𝑎₂₂𝑐₂−12𝑥₂(𝑡)−𝑥^∗₂²+𝑎₃₁+𝑎₃₂−𝑎₃₃𝑐₃−12𝑥₃(𝑡)−𝑥^∗₃²+𝑐₁𝑥^∗₁𝜎²₁+𝑐₂𝑥^∗₂𝜎²₂+𝑐₃𝑥^∗₃𝜎²₃2𝑑𝑡+³_𝑖=1𝜎_𝑖𝑑𝐵_𝑖. (A.11) Integrating from 0 to 𝑡, taking the expectation, we have 𝐸[𝑉(𝑡)]−𝑉(0)≤−𝐸^𝑡₀𝑎₁₁−𝑎₁₂−𝑎₁₃𝑐₁+12𝑥₁(𝑠)−𝑥^∗₁²+𝑎₂₂−𝑎₂₁−𝑎₂₃𝑐₂+12𝑥₂(𝑠)−𝑥^∗₂²+𝑎₃₃−𝑎₃₁−𝑎₃₂𝑐₃+12𝑥₃(𝑠)−𝑥^∗₃²+𝑐𝑑𝑠₁𝑥^∗₁𝜎²₁+𝑐₂𝑥^∗₂𝜎²₂+𝑐₃𝑥^∗₃𝜎²₃2𝑡. (A.12) Then we yield 𝐸[𝑉(𝑡)]𝑡≤𝑉(0)𝑡−𝑎₁₁−𝑎₁₂−𝑎₁₃𝑐₁+12𝐸1𝑡^𝑡₀𝑥₁(𝑠)−𝑥^∗₁²−𝑎𝑑𝑠₂₂−𝑎₂₁−𝑎₂₃𝑐₂+12𝐸1𝑡^𝑡₀𝑥₂(𝑠)−𝑥^∗₂²−𝑎𝑑𝑠₃₃−𝑎₃₁−𝑎₃₂𝑐₃+12𝐸1𝑡^𝑡₀[𝑥₃(𝑠)−𝑥^∗₃]²+𝑐𝑑𝑠₁𝑥^∗₁𝜎²₁+𝑐₂𝑥^∗₂𝜎²₂+𝑐₃𝑥^∗₃𝜎²₃2. (A.13) Letting 𝑡→∞, therefore we have limsup_𝑡→∞1𝑡𝐸^𝑡₀𝑎₁₁−𝑎₁₂−𝑎₁₃𝑐₁+12𝑥₁(𝑠)−𝑥^∗₁²+𝑎₂₂−𝑎₂₁−𝑎₂₃𝑐₂+12𝑥₂(𝑠)−𝑥^∗₂²+𝑎₃₃−𝑎₃₁−𝑎₃₂𝑐₃+12𝑥₃(𝑠)−𝑥^∗₃²≤∑𝑑𝑠_3𝑖=1𝑐_𝑖𝑥^∗_𝑖𝜎²_𝑖2, (A.14) which is the required assertion. The proof is completed.

B. Proof of Theorem 11

Proof. It is sufficient to prove 𝑥^∗_𝑖≤liminf_𝑡→∞1𝑡^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠≤limsup_𝑡→∞1𝑡^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠≤𝑥^∗_𝑖a.s.𝑖=1,2,3. (B.1) By Itô’s formula, we have 𝑑log𝑥₁(𝑡)=𝑟₁−𝜎²₁2−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡−𝜏)+𝑎₁₃𝑥₃(𝑡−𝜏)𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡),𝑑log𝑥₂(𝑡)=𝑟₂−𝜎²₂2+𝑎₂₁𝑥₁(𝑡−𝜏)−𝑎₂₂𝑥₂(𝑡)+𝑎₂₃𝑥₃(𝑡−𝜏)𝑑𝑡+𝜎₂𝑑𝐵₂(𝑡),𝑑log𝑥₃(𝑡)=𝑟₃−𝜎²₃2+𝑎₃₁𝑥₁(𝑡−𝜏)+𝑎₃₂𝑥₂(𝑡−𝜏)−𝑎₃₃𝑥₃(𝑡)𝑑𝑡+𝜎₃𝑑𝐵₃(𝑡). (B.2) Integrating both sides of (B.2) from 0 to 𝑡, then we have log𝑥₁(𝑡)=log𝑥₁(0)+𝑟₁𝑡+𝑎₁₂⁰_−𝜏𝜉₂(𝑠)𝑑𝑠+𝑎₁₃⁰_−𝜏𝜉₃(𝑠)𝑑𝑠−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑎₁₂₀^𝑡−𝜏𝑥₂(𝑠)𝑑𝑠+𝑎₁₃₀^𝑡−𝜏𝑥₃(𝑠)𝑑𝑠+𝜎₁𝐵₁(𝑡),log𝑥₂(𝑡)=log𝑥₂(0)+𝑟₂𝑡+𝑎₂₁⁰_−𝜏𝜉₁(𝑠)𝑑𝑠+𝑎₂₃⁰_−𝜏𝜉₃(𝑠)𝑑𝑠+𝑎₂₁₀^𝑡−𝜏𝑥₁(𝑠)𝑑𝑠−𝑎₂₂^𝑡₀𝑥₂(𝑠)𝑑𝑠+𝑎₂₃₀^𝑡−𝜏𝑥₃(𝑠)𝑑𝑠+𝜎₂𝐵₂(𝑡),log𝑥₃(𝑡)=log𝑥₃(0)+𝑟₃𝑡+𝑎₃₁⁰_−𝜏𝜉₁(𝑠)𝑑𝑠+𝑎₃₂⁰_−𝜏𝜉₂(𝑠)𝑑𝑠+𝑎₃₁₀^𝑡−𝜏𝑥₁(𝑠)𝑑𝑠+𝑎₃₂₀^𝑡−𝜏𝑥₂(𝑠)𝑑𝑠−𝑎₃₃^𝑡₀𝑥₃(𝑠)𝑑𝑠+𝜎₃𝐵₃(𝑡), (B.3) where 𝑟_𝑖=𝑟_𝑖−(𝜎²_𝑖/2),𝑖=1,2,3. From (33), we know liminf_𝑡→∞1𝑡^𝑡₀𝑥₁(𝑠)𝑑𝑠≥𝑟₁𝑎₁₁∶=𝑀₁,liminf_𝑡→∞1𝑡^𝑡₀𝑥₂(𝑠)𝑑𝑠≥𝑟₂𝑎₂₂∶=𝑁₁,liminf_𝑡→∞1𝑡^𝑡₀𝑥₃(𝑠)𝑑𝑠≥𝑟₃𝑎₃₃∶=𝐿₁. (B.4) Hence, for any sufficiently small 0<𝜀₁<(1/2)min{𝑀₁,𝑁₁,𝐿₁}, there is a 𝑇₁(𝜔)>0 such that if 𝑡>𝑇₁(𝜔), (𝑡−𝜏)/𝑡≥1−𝜀₁, (1/𝑡)∫^𝑡₀𝑥₁(𝑠)𝑑𝑠≥𝑀₁−𝜀₁, (1/𝑡)∫^𝑡₀𝑥₂(𝑠)𝑑𝑠≥𝑁₁−𝜀₁, (1/𝑡)∫^𝑡₀𝑥₃(𝑠)𝑑𝑠≥𝐿₁−𝜀₁ It follows from (B.3) that for 𝑡>𝑇₁(𝜔)+𝜏, log𝑥₁(𝑡)=log𝑥₁(0)+𝑟₁𝑡+𝑎₁₂⁰_−𝜏𝑥₂(𝑠)𝑑𝑠+𝑎₁₃⁰_−𝜏𝑥₃(𝑠)𝑑𝑠−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑎₁₂₀^𝑡−𝜏𝑥₂(𝑠)𝑑𝑠+𝑎₁₃₀^𝑡−𝜏𝑥₃(𝑠)𝑑𝑠+𝜎₁𝐵₁(𝑡)=log𝑥₁(0)+𝑟₁𝑡+𝑎₁₂⁰_−𝜏𝑥₂(𝑠)𝑑𝑠+𝑎₁₃⁰_−𝜏𝑥₃(𝑠)𝑑𝑠−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑎₁₂()1𝑡−𝜏𝑡−𝜏₀^𝑡−𝜏𝑥₂(𝑠)𝑑𝑠+𝑎₁₃()1𝑡−𝜏𝑡−𝜏₀^𝑡−𝜏𝑥₃(𝑠)𝑑𝑠+𝜎₁𝐵₁(𝑡)≥log𝑥₁(0)+𝑟₁𝑡+𝑎₁₂⁰_−𝜏𝑥₂(𝑠)𝑑𝑠+𝑎₁₃⁰_−𝜏𝑥₃(𝑠)𝑑𝑠−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑎₁₂()×𝑁𝑡−𝜏₁−𝜀₁+𝑎₁₃(𝑡−𝜏)𝐿₁−𝜀₁+𝜎₁𝐵₁(𝑡)=𝑟₁+𝑎₁₂𝑁₁−𝜀₁+𝑎₁₃𝐿₁−𝜀₁𝑡−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝜎₁𝐵₁(𝑡)+𝑆₁, (B.5) where 𝑆₁=log𝑥₁(0)−𝑎₁₂𝑁₁−𝜀₁+𝑎₁₃𝐿₁−𝜀₁𝜏+𝑎₁₂⁰_−𝜏𝑥₂(𝑠)𝑑𝑠+𝑎₁₃⁰_−𝜏𝑥₃(𝑠)𝑑𝑠. (B.6) We know lim_𝑡→∞𝜎₁𝐵₁(𝑡)+𝑆₁𝑡=0. (B.7) From Lemma 10, we have liminf_𝑡→∞1𝑡^𝑡₀𝑥₁(𝑠)≥𝑑𝑠𝑟₁+𝑎₁₂𝑁₁−𝜀₁+𝑎₁₃𝐿₁−𝜀₁𝑎₁₁∶=𝑀₂. (B.8) Similarly, liminf_𝑡→∞1𝑡^𝑡₀𝑥₂(𝑠)≥𝑑𝑠𝑟₂+𝑎₂₁𝑀₁−𝜀₁+𝑎₂₃𝐿₁−𝜀₁𝑎₂₂∶=𝑁₂,liminf_𝑡→∞1𝑡^𝑡₀𝑥₃(𝑠)≥𝑑𝑠𝑟₃+𝑎₃₁𝑀₁−𝜀₁+𝑎₃₂𝑁₁−𝜀₁𝑎₃₃∶=𝐿₂. (B.9) Let 𝜀₂=(1/2)min{𝑀₂,𝑁₂,𝐿₂,𝜀₁}, continuing this process, we obtain sequences {𝑀_𝑛},{𝑁_𝑛},{𝐿_𝑛},{𝜀_𝑛}, where 𝑀_𝑛=𝑟₁+𝑎₁₂𝑁_𝑛−1−𝜀_𝑛−1+𝑎₁₃𝐿_𝑛−1−𝜀_𝑛−1𝑎₁₁,𝑁_𝑛=𝑟₂+𝑎₂₁𝑀_𝑛−1−𝜀_𝑛−1+𝑎₂₃𝐿_𝑛−1−𝜀_𝑛−1𝑎₂₂,𝐿_𝑛=𝑟₃+𝑎₃₁𝑀_𝑛−1−𝜀_𝑛−1+𝑎₃₂𝑁_𝑛−1−𝜀_𝑛−1𝑎₃₃,𝜀_𝑛=12min𝑀_𝑛,𝑁_𝑛,𝐿_𝑛,𝜀_𝑛−1. (B.10) Sequence {𝑀_𝑛},{𝑁_𝑛},{𝐿_𝑛} is nondecreasing and bounded; then we have liminf_𝑡→∞1𝑡^𝑡₀𝑥₁(𝑠)𝑑𝑠≥𝑀_𝑛,liminf_𝑡→∞1𝑡^𝑡₀𝑥₂(𝑠)𝑑𝑠≥𝑁_𝑛,liminf_𝑡→∞1𝑡^𝑡₀𝑥₃(𝑠)𝑑𝑠≥𝐿_𝑛. (B.11) Let 𝑛→∞; we derive from (B.10) that lim_𝑛→∞𝜀_𝑛=0,lim_𝑛→∞𝑀_𝑛=𝑥^∗₁,lim_𝑛→∞𝑁_𝑛=𝑥^∗₂,lim_𝑛→∞𝐿_𝑛=𝑥^∗₃, (B.12) where 𝑥^∗₁,𝑥^∗₂,𝑥^∗₃ is the unique solution of the following equation: 𝑟₁−𝑎₁₁𝑥^∗₁+𝑎₁₂𝑥^∗₂+𝑎₁₃𝑥^∗₃=0,𝑟₂+𝑎₂₁𝑥^∗₁−𝑎₂₂𝑥^∗₂+𝑎₂₃𝑥^∗₃=0,𝑟₃+𝑎₃₁𝑥^∗₁+𝑎₃₂𝑥^∗₂−𝑎₃₃𝑥^∗₃=0. (B.13) Thus we obtain the assertion liminf_𝑡→∞1𝑡^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠≥𝑥^∗_𝑖𝑖=1,2,3a.s. (B.14) Next, we will prove limsup_𝑡→∞1𝑡^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠≤𝑥^∗_𝑖𝑖=1,2,3a.s. (B.15) It follows from (B.2) that 𝑑log𝑥₁(𝑡)+𝑎_12𝑡^𝑡+𝜏𝑥₂(𝑠−𝜏)𝑑𝑠+𝑎_13𝑡^𝑡+𝜏𝑥₃()=𝑠−𝜏𝑑𝑠𝑟₁−𝑎₁₁𝑥₁(𝑡)+𝑎₁₂𝑥₂(𝑡)+𝑎₁₃𝑥₃(𝑡)𝑑𝑡+𝜎₁𝑑𝐵₁(𝑡). (B.16) Integrating both sides of (B.16) from 0 to 𝑡, we have log𝑥₁(𝑡)+𝑎_12𝑡^𝑡+𝜏𝑥₂(𝑠−𝜏)𝑑𝑠−𝑎₁₂^𝜏₀𝑥₂(𝑠−𝜏)𝑑𝑠+𝑎_13𝑡^𝑡+𝜏𝑥₃(𝑠−𝜏)𝑑𝑠−𝑎₁₃^𝜏₀𝑥₃(𝑠−𝜏)𝑑𝑠=log𝑥₁(0)+𝑟₁𝑡−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑎₁₂^𝑡₀𝑥₂(𝑠)𝑑𝑠+𝑎₁₃^𝑡₀𝑥₃(𝑠)𝑑𝑠+𝜎₁𝐵₁(𝑡). (B.17) Since 𝑥_𝑖(𝑡)>0,𝑖=1,2,3, we have log𝑥₁(𝑡)≤𝑈₁+𝑟₁𝑡−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑎₁₂^𝑡₀𝑥₂(𝑠)𝑑𝑠+𝑎₁₃^𝑡₀𝑥₃(𝑠)𝑑𝑠+𝜎₁𝐵₁(𝑡), (B.18) where 𝑈₁=𝑎₁₂^𝜏₀𝑥₂(𝑠−𝜏)𝑑𝑠+𝑎₁₃^𝜏₀𝑥₃(𝑠−𝜏)𝑑𝑠+log𝑥₁(0). (B.19) Similarly, log𝑥₂(𝑡)≤𝑈₂+𝑟₂𝑡−𝑎₂₂^𝑡₀𝑥₂(𝑠)𝑑𝑠+𝑎₂₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑎₂₃^𝑡₀𝑥₃(𝑠)𝑑𝑠+𝜎₂𝐵₂(𝑡), (B.20) where 𝑈₂=𝑎₂₁^𝜏₀𝑥₁(𝑠−𝜏)𝑑𝑠+𝑎₂₃^𝜏₀𝑥₃(𝑠−𝜏)𝑑𝑠+log𝑥₂(0).log𝑥₃(𝑡)≤𝑈₃+𝑟₃𝑡−𝑎₃₃^𝑡₀𝑥₃(𝑠)𝑑𝑠+𝑎₃₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑎₃₂^𝑡₀𝑥₂(𝑠)𝑑𝑠+𝜎₃𝐵₃(𝑡), (B.21)(B.22) where 𝑈₃=𝑎₃₁^𝜏₀𝑥₁(𝑠−𝜏)𝑑𝑠+𝑎₃₂^𝜏₀𝑥₂(𝑠−𝜏)𝑑𝑠+log𝑥₃(0). (B.23) Let 𝑏₁₁ be a positive constant 𝑏₁₂=−𝑎₁₂𝑎₃₃−𝑎₁₃𝑎₃₂𝑎₂₃𝑎₃₂−𝑎₃₃𝑎₂₂𝑏₁₁,𝑏₁₃=−𝑎₂₃𝑎₁₂−𝑎₁₃𝑎₂₂𝑎₂₃𝑎₃₂−𝑎₃₃𝑎₂₂𝑏₁₁,𝑚₁=𝑎₁₁𝑏₁₁−𝑎₂₁−𝑎₁₂𝑎₃₃−𝑎₁₃𝑎₃₂𝑎₂₃𝑎₃₂−𝑎₃₃𝑎₂₂𝑏₁₁−𝑎₃₁−𝑎₂₃𝑎₁₂−𝑎₁₃𝑎₂₂𝑎₂₃𝑎₃₂−𝑎₃₃𝑎₂₂𝑏₁₁=−𝑎₁₁𝑎₂₂𝑎₃₃+𝑎₁₁𝑎₂₃𝑎₃₂+𝑎₂₂𝑎₁₃𝑎₃₁+𝑎₃₃𝑎₁₂𝑎₂₁+𝑎₁₂𝑎₂₃𝑎₃₁+𝑎₂₁𝑎₃₂𝑎₁₃𝑎₂₃𝑎₃₂−𝑎₃₃𝑎₂₂⁻¹𝑏₁₁. (B.24) From Lemma 3, we know 𝑏₁₂>0,𝑏₁₃>0,𝑚₁>0. Since 𝑥_𝑖(𝑡)>0, from (B.18)–(B.22) we have 𝑏₁₁log𝑥₁(𝑡)+𝑏₁₂log𝑥₂(𝑡)+𝑏₁₃log𝑥₃(𝑡)=𝑏₁₁𝑟₁+𝑏₁₂𝑟₂+𝑏₁₃𝑟₃𝑡−𝑚₁^𝑡₀𝑥₁(𝑠)+𝑑𝑠_3𝑖=1𝑏_1𝑖𝜎_𝑖𝐵_𝑖(𝑡)+³_𝑖=1𝑏_1𝑖𝑈_𝑖. (B.25) From (33), we know liminf_𝑡→∞log𝑥_𝑖(𝑡)𝑡≥0𝑖=1,2,3, (B.26) which implies lim_𝑡→∞𝑏₁₁log𝑥₁(𝑡)+𝑏₁₂log𝑥₂(𝑡)+𝑏₁₃log𝑥₃(𝑡)𝑡≥0. (B.27) From lim_𝑡→∞∑³_𝑖=1𝑏_1𝑖𝜎_𝑖𝐵_𝑖(𝑡)+∑³_𝑖=1𝑏_1𝑖𝑈_𝑖𝑡=0, (B.28) we have limsup_𝑡→∞1𝑡^𝑡₀𝑥₁(𝑠)𝑏𝑑𝑠≤¹¹𝑟₁+𝑏₁₂𝑟₂+𝑏₁₃𝑟₃𝑚₁=𝑥^∗₁. (B.29) Similarly, limsup_𝑡→∞1𝑡^𝑡₀𝑥₂(𝑠)𝑑𝑠≤𝑥^∗₂,limsup_𝑡→∞1𝑡^𝑡₀𝑥₃(𝑠)𝑑𝑠≤𝑥^∗₃. (B.30) Therefore, limsup_𝑡→∞(1/𝑡)∫^𝑡₀𝑥(𝑠)𝑑𝑠≤𝑥^∗. The proof is completed.

Conflict of Interests

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

Acknowledgments

The work was supported by the NSFC of China (no. 11371085) and the Ph.D. Programs Foundation of Ministry of China (no. 200918).