A stochastic predator–prey model with two competitive preys and Ornstein–Uhlenbeck process

Abstract

In this paper, a stochastic predator–prey model with two competitive preys and Ornstein–Uhlenbeck process is formulated and analysed, which is used to obtain a better understanding of the population dynamics. At first, we validate that the stochastic system has a unique global solution with any initial value. Then we analyse the stochastic dynamics of the model in detail, including pth moment boundedness, asymptotic pathwise estimation in turn. After that, we obtain sufficient conditions for the existence of a stationary distribution of the system by adopting stochastic Lyapunov function methods. In addition, under some mild conditions, we derive the specific form of covariance matrix in the probability density near the quasi-positive equilibrium of the stochastic system. Finally, numerical illustrative examples are depicted to confirm our theoretical findings.

Keywords:

• Three-species predator–prey model
• competition between preys
• ornstein–Uhlenbeck process
• moment boundedness
• stationary distribution
• density function

1. Introduction

Predator–prey interaction with multiple prey is common in nature, for example, Pelican birds in Salton sea eat Tilapia Oreochromis mossambicus with Avian botulism [3]. In order to study the evolution law of predator–prey interaction with multiple preys, a lot of predator–prey models with multiple preys were established and analysed during the past few decades [2, 14, 16, 23, 29]. The well-known predator–prey model with competitive interaction between two preys can be expressed as follows: (1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) where $a_i$ and $b_{ij}$ $(i,j=1,2,3)$ are all positive constants, $x_1(t)$ and $x_2(t)$ represent the population densities of two competitive species and $x_3(t)$ stands for the population density of predator species. The biological significance of the parameters is as follows: $a_1$ and $a_2$ denote the intrinsic growth rates of the two competitive preys and $a_3$ represents the natural death rate of the predator population. The coefficients $b_{11}$, $b_{22}$ and $b_{33}$ denote the intra-specific competitive coefficients of two prey species $x_1$, $x_2$ and the predator species $x_3$, respectively. $b_{13}$ and $b_{23}$ represent the capturing rate of the predator at time interval, $b_{12}$ and $b_{21}$ are the inter-specific competitive coefficients of two prey species $x_1$, $x_2$ and $b_{31}$, $b_{32}$ stand for the rate of conversion of nutrients of the predator species. All the parameters in the model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ) are assumed to be positive constants.

On the other hand, in the natural world, the growth of population is always disturbed by environmental noises [4, 5], such as random changes in epidemics, droughts, floods, temperatures, etc. These changes sometimes are so strong that they can change the behaviours of population significantly [32]. May [24] pointed out that due to environmental fluctuations, the birth rates, death rates, competition coefficients, transmission coefficients and other parameters involved in the system should exhibit random fluctuation to some extent. Therefore, the parameters $a_1$, $a_2$ and $a_3$ should be regarded as random variables involved in randomness. Accordingly, many scholars have introduced white noise in the deterministic population systems to study how the noise affects the dynamics of stochastic population models [9, 10, 12, 13, 15, 19, 21, 25] and this way is common.

However, except using the common way to simulate the random disturbances, we can also assume that the parameters involved in the system satisfy the well-known Ornstein–Uhlenbeck process. The Ornstein–Uhlenbeck process is an Itô process which can be described by the following stochastic differential equation: $\mathrm{d}a(t)=\theta [\mu -a(t)]\mathrm{d}t+\sigma \mathrm{d}B(t),$where θ, μ, σ are all positive constants and θ is a mean reversion constant that denotes the strength of the process attracted by the mean [26], μ is the mean value, σ represents the degree of volatility around the mean value μ caused by shocks, $B(t)$ denotes the standard Brownian motion.

By the work of Allen [1], the stochastic mean-reverting process should incorporate the influences of environmental perturbations into the parameters of a system. So in this paper, by considering the second way of introducing randomness in model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ), we assume that $a_1(t)$, $a_2(t)$ and $a_3(t)$ are affected by three Ornstein–Uhlenbeck processes, i.e. (2) $\begin{array}{rl}\mathrm{d}{a}_1(t)& ={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\end{array}$(2) (3) $\begin{array}{rl}\mathrm{d}{a}_2(t)& ={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\end{array}$(3) (4) $\begin{array}{rl}\mathrm{d}{a}_3(t)& ={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t),\end{array}$(4) where ${\alpha }_i$, ${\overline{a}}_i$, ${\sigma }_i$ are all positive constants, and ${\alpha }_i$ denotes the speed of reversion; ${\sigma }_i$ are the intensities of the volatility of process $a_i$; ${\overline{a}}_i$ represent the mean reversion levels of $a_i(t)$, i = 1, 2, 3. $\{B_1(t){\}}_{t\ge 0}$, $\{B_2(t){\}}_{t\ge 0}$ and $\{B_3(t){\}}_{t\ge 0}$ are mutually independent standard Brownian motions defined on a complete probability space $(\mathrm{\Omega },\mathcal{F},\{{\mathcal{F}}_t{\}}_{t\ge 0},\mathbb{P})$ with a filtration $\{{\mathcal{F}}_t{\}}_{t\ge 0}$ satisfying the usual conditions [22].

By performing stochastic integral operation to (2(2) $\begin{array}{rl}\mathrm{d}{a}_1(t)& ={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\end{array}$(2) )–(4(4) $\begin{array}{rl}\mathrm{d}{a}_3(t)& ={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t),\end{array}$(4) ), respectively, we have $a_i(t)={\overline{a}}_i+(a_i^0-{\overline{a}}_i){\mathrm{e}}^{-{\alpha }_{i}t}+{\sigma }_i{\int }_0^t{\mathrm{e}}^{-{\alpha }_i(t-s)}\mathrm{d}{B}_i(s),i=1,2,3,$where $a_i(0)=a_i^0$ are the initial values of Ornstein–Uhlenbeck processes $a_i(t)$, i = 1, 2, 3. It is easy to see that the expected values and the variances of $a_i(t)$ over an interval $[0,t]$ are as follows $\mathbb{E}[a_i(t)]={\overline{a}}_i+(a_i^0-{\overline{a}}_i){\mathrm{e}}^{-{\alpha }_{i}t},\text{V}ar[a_i(t)]=\frac{{\sigma }_i^2}{2{\alpha }_i}(1-{\mathrm{e}}^{-2{\alpha }_{i}t}),i=1,2,3.$According to the property of Brownian motion, we obtain that ${\sigma }_i{\int }_0^t{\mathrm{e}}^{-{\alpha }_i(t-s)}\mathrm{d}{B}_i(s)$ follows the distribution $N(0,{\sigma }_i^2(1-{\mathrm{e}}^{-2{\alpha }_{i}t})/2{\alpha }_i)$ and $a_i$ follows the normal distribution $N({\overline{a}}_i+(a_i^0-{\overline{a}}_i){\mathrm{e}}^{-{\alpha }_{i}t},{\sigma }_i^2(1-{\mathrm{e}}^{-2{\alpha }_{i}t})/2{\alpha }_i)$. Apparently, the Ornstein–Uhlenbeck processes $a_i(t)$ follow the distribution $N(a_i^0,0)$ as t tends to zero, i = 1, 2, 3, which can be seen as a both mathematically and biologically reasonable method to simulate the random impacts of key parameters in a population system. Keeping in mind of the above facts and combining with (2(2) $\begin{array}{rl}\mathrm{d}{a}_1(t)& ={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\end{array}$(2) )–(4(4) $\begin{array}{rl}\mathrm{d}{a}_3(t)& ={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t),\end{array}$(4) ), we get the stochastic form corresponding to system (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ) which is presented as follows: (5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) Here the parameters involved in system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) are the same as those in system (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ).

The rest of the paper is arranged as follows. In Section 2, the existence and uniqueness of the global solution of the stochastic system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) is proved which is necessary to study the dynamics of a population model. In Section 3, the boundedness of the pth moment of solutions to system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) is shown. In Section 4, the asymptotic pathwise estimation of the solutions of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) is investigated. In Section 5, sufficient criteria for the existence of a stationary distribution of positive solutions to system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) are derived. In Section 6, we gain the exact form of covariance matrix in the probability density around the origin point. In Section 7, through the numerical simulations of several examples, we validate the effectiveness of our theoretical findings. Finally, we present a brief conclusion to end the paper.

Next, we present the following notations which will be used in the whole paper: $\begin{array}{rl}& {\mathbb{R}}_{+}^d=\{x=(x_1,\dots ,x_d)\in {\mathbb{R}}^d:x_i>0,1\le i\le d\}\mathrm{and}\\ & {\overline{\mathbb{R}}}_{+}^d=\{x=(x_1,\dots ,x_d)\in {\mathbb{R}}^d:x_i\ge 0,1\le i\le d\}.\\ & x\wedge y=min\{x,y\}\mathrm{for}\mathrm{any}x,y\in \mathbb{R}.\end{array}$Let $C^2({\mathbb{R}}^d;{\overline{\mathbb{R}}}_{+})$ be the set of all nonnegative functions $V(x)$ defined on ${\mathbb{R}}^d$ which means it contains all continuously twice differentiable in x. We use the notation ${\mathbf{I}}_{\mathbb{A}}$ to denote the indicator function of the set $\mathbb{A}$. If $\mathbb{G}$ is a matrix, we use the notation ${\mathbb{G}}^T$ to represent its transpose. If $\mathbb{G}$ is an invertible matrix, its inverse matrix is denoted by ${\mathbb{G}}^{-1}$. If G is a square matrix, we use the notation $|G|$ to denote its determinant. Let ${\mathbb{G}}^{(k)}$ be the kth sequential principal submatrix of the matrix $\mathbb{G}$. In addition, for any two d-dimensional symmetric matrices $\mathbb{H}$ and $\mathbb{J}$, we define $\begin{array}{rl}& \mathbb{H}\succ \mathbf{0}:\mathbb{H}\mathrm{is}\mathrm{a}\mathrm{positive}\mathrm{definite}\mathrm{matrix},\mathbb{H}⪰\mathbf{0}:\mathbb{H}\mathrm{is}\mathrm{a}\mathrm{positive}\mathrm{semi}-\mathrm{definite}\mathrm{matrix},\\ & \mathbb{H}⪰\mathbb{J}:\mathbb{H}-\mathbb{J}\mathrm{is}\mathrm{at}\mathrm{least}\mathrm{a}\mathrm{positive}\mathrm{semi}-\mathrm{definite}\mathrm{matrix}.\end{array}$In this sense, $\mathbb{H}\succ \mathbf{0}$ if $\mathbb{J}\succ \mathbf{0}$ and $\mathbb{H}⪰\mathbb{J}$.

Throughout this paper, the following hypothesis always holds: $\begin{array}{rl}\mathbf{(}\mathbf{A}\mathbf{)}& b_{12}{b}_{21}{b}_{33}+b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}-b_{11}{b}_{22}{b}_{33}-b_{11}{b}_{23}{b}_{32}-b_{13}{b}_{22}{b}_{31}<0,\\ & {\overline{a}}_2{b}_{12}{b}_{33}+{\overline{a}}_2{b}_{13}{b}_{32}+{\overline{a}}_3{b}_{12}{b}_{23}-{\overline{a}}_1{b}_{22}{b}_{33}-{\overline{a}}_1{b}_{23}{b}_{32}-{\overline{a}}_3{b}_{13}{b}_{22}<0,\\ & {\overline{a}}_1{b}_{21}{b}_{33}+{\overline{a}}_1{b}_{23}{b}_{31}+{\overline{a}}_3{b}_{13}{b}_{21}-{\overline{a}}_2{b}_{11}{b}_{33}-{\overline{a}}_2{b}_{13}{b}_{31}-{\overline{a}}_3{b}_{11}{b}_{23}<0,\\ & {\overline{a}}_1{b}_{21}{b}_{32}+{\overline{a}}_2{b}_{12}{b}_{31}+{\overline{a}}_3{b}_{11}{b}_{22}-{\overline{a}}_1{b}_{22}{b}_{31}-{\overline{a}}_2{b}_{11}{b}_{32}-{\overline{a}}_3{b}_{12}{b}_{21}<0.\end{array}$

Remark 1.1

The Assumption (A) is used to ensure the positive equilibrium $\overline{E}({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3)$ of system (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ) exists.

Now we will present basic theories on stochastic differential equations (for more details, we refer the reader to [22]).

Let $X(t)$ be a d-dimensional stochastic process on $t\ge 0$ given by the following stochastic differential equation: $\mathrm{d}X(t)=f(X(t))\mathrm{d}t+g(X(t))\mathrm{d}B(t),$where $f\in L^1({\mathbb{R}}_{+};{\mathbb{R}}^d)$, $g\in L^2({\mathbb{R}}_{+};{\mathbb{R}}^{d\times m})$ and $B(t)=\{(B_t^1,B_t^2,\dots ,B_t^d{)}^T{\}}_{t\ge 0}$ is a d-dimensional Brownian motion defined on the complete probability space $(\mathrm{\Omega },\mathcal{F},\{{\mathcal{F}}_t{\}}_{t\ge 0},\mathbb{P})$. We define the differential operator L as follows: $L=\sum _{i=1}^d{f}_i(X(t))\frac{\mathrm{\partial }}{\mathrm{\partial }{X}_i}+\frac{1}{2}\sum _{i,j=1}^d(g^T(X(t))g(X(t)){)}_{ij}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{X}_i\mathrm{\partial }{X}_j}.$If L acts on a function $V\in C^2({\mathbb{R}}^d;\mathbb{R})$, then $LV(X(t))=V_X(X(t))f(X(t))+\frac{1}{2}\mathrm{trace}[g^T(X(t))V_{XX}(X(t))g(X(t))],$where $V_X=(\mathrm{\partial }V/\mathrm{\partial }{X}_1,\dots ,\mathrm{\partial }V/\mathrm{\partial }{X}_d)$ and $V_{XX}=({\mathrm{\partial }}^{2}V/(\mathrm{\partial }{X}_i\mathrm{\partial }{X}_j){)}_{d\times d}$. If $X(t)\in {\mathbb{R}}^d$, then by means of Itô's formula, we gain $\mathrm{d}V(X(t))=LV(X(t))\mathrm{d}t+V_X(X(t))g(X(t))\mathrm{d}B(t).$

2. Existence and uniqueness of the global solution

To investigate the dynamical behaviour of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ), we should first ensure whether the solution is global or not. To this end, we present the following theorem which is related to the existence and uniqueness of the global solution to system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ).

Theorem 2.1

Given initial value $(x_1(0),x_2(0),x_3(0),a_1(0),a_2(0),a_3(0){)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$, there is a unique global solution $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T$ to system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) which will remain in ${\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$ almost surely (a.s.).

Proof.

It is noticed that all the coefficients of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) are locally Lipschitz continuous, then for any initial value $(x_1(0),x_2(0),x_3(0),a_1(0),a_2(0),a_3(0){)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$, there exists a unique local solution $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) on the interval $[0,{\tau }_e)$ a.s., where ${\tau }_e$ represents the explosion time [22]. Next, let $n_0\ge 1$ be large enough such that $x_1(0)$, $x_2(0)$, $x_3(0)$, ${\mathrm{e}}^{a_1(0)}$, ${\mathrm{e}}^{a_2(0)}$ and ${\mathrm{e}}^{a_3(0)}$ all lie within the interval $[1/n_0,n_0]$. For each integer $n\ge n_0$, define the stopping time as below [22] $\begin{array}{rl}& {\tau }_n=inf\{t\in [0,{\tau }_e):min\{x_1(t),x_2(t),x_3(t),{\mathrm{e}}^{a_1(t)},{\mathrm{e}}^{a_2(t)},{\mathrm{e}}^{a_3(t)}\}\\ & \le \frac{1}{n}\mathrm{or}max\{x_1(t),x_2(t),x_3(t),{\mathrm{e}}^{a_1(t)},{\mathrm{e}}^{a_2(t)},{\mathrm{e}}^{a_3(t)}\}\ge n\}.\end{array}$It is clear that ${\tau }_n$ is increasing as $n\to \mathrm{\infty }$. Denote by ${\tau }_{\mathrm{\infty }}=\underset{n\to \mathrm{\infty }}{lim}{\tau }_n$, then ${\tau }_{\mathrm{\infty }}\le {\tau }_e$ a.s. If ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ a.s. is true, then ${\tau }_e=\mathrm{\infty }$ a.s. and $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$ a.s. In other words, to complete the proof, ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ a.s. needs to be proved. If this assertion is not correct, then there are two constants T>0 and $ϵ\in (0,1)$ such that $\mathbb{P}\{{\tau }_{\mathrm{\infty }}\le T\}>ϵ.$Hence, there exists an integer $n_1\ge n_0$ such that (6) $\mathbb{P}\{{\tau }_n\le T\}\ge ϵ\mathrm{for}\mathrm{all}n\ge n_1.$(6) Define a $C^2$-function $V:{\mathbb{R}}_{+}^3\times {\mathbb{R}}^3\to {\overline{\mathbb{R}}}_{+}$ by $\begin{array}{rl}V(x_1,x_2,x_3,a_1,a_2,a_3)& =c_1(x_1-1-\ln x_1)+c_2(x_2-1-\ln x_2)+(x_3-1-\ln x_3)\\ & +a_1^4+a_2^4+a_3^4,\end{array}$where $c_1$ and $c_2$ are positive constants to be determined later. The nonnegativity of this function can be derived from $u-1-\ln u\ge 0,\mathrm{\forall }u>0.$Applying Itô's formula [22] to V, we obtain (7) $\begin{array}{rl}\mathrm{d}V(x_1,x_2,x_3,a_1,a_2,a_3)& =LV(x_1,x_2,x_3,a_1,a_2,a_3)\mathrm{d}t+4{\sigma }_1{a}_1^3\mathrm{d}{B}_1(t)+4{\sigma }_2{a}_2^3\mathrm{d}{B}_2(t)\\ & +4{\sigma }_3{a}_3^3\mathrm{d}{B}_3(t),\end{array}$(7) where $LV:{\mathbb{R}}_{+}^3\times {\mathbb{R}}^3\to \mathbb{R}$ is defined by $\begin{array}{rl}LV& =c_1{x}_1(a_1-b_{11}{x}_1-b_{12}{x}_2-b_{13}{x}_3)+c_2{x}_2(a_2-b_{21}{x}_1-b_{22}{x}_2-b_{23}{x}_3)\\ & +x_3(-a_3+b_{31}{x}_1+b_{32}{x}_2-b_{33}{x}_3)\\ & -c_1{a}_1+c_1{b}_{11}{x}_1+c_1{b}_{12}{x}_2+c_1{b}_{13}{x}_3-c_2{a}_2+c_2{b}_{21}{x}_1+c_2{b}_{22}{x}_2+c_2{b}_{23}{x}_3\\ & +a_3-b_{31}{x}_1-b_{32}{x}_2+b_{33}{x}_3\\ & +4{\alpha }_1{\overline{a}}_1{a}_1^3+6{\sigma }_1^2{a}_1^2+4{\alpha }_2{\overline{a}}_2{a}_2^3+6{\sigma }_2^2{a}_2^2+4{\alpha }_3{\overline{a}}_3{a}_3^3\\ & +6{\sigma }_3^2{a}_3^2-4{\alpha }_1{a}_1^4-4{\alpha }_2{a}_2^4-4{\alpha }_3{a}_3^4\\ & \le c_1{a}_1{x}_1+c_2{a}_2{x}_2-a_3{x}_3-c_1{b}_{11}{x}_1^2-c_2{b}_{22}{x}_2^2-b_{33}{x}_3^2\\ & +(b_{31}-c_1{b}_{13})x_1{x}_3+(b_{32}-c_2{b}_{23})x_2{x}_3-c_1{a}_1+a_3\\ & -c_2{a}_2+(c_1{b}_{11}+c_2{b}_{21})x_1+(c_1{b}_{12}+c_2{b}_{22})x_2+(c_1{b}_{13}+c_2{b}_{23}+b_{33})x_3\\ & +4{\alpha }_1{\overline{a}}_1{a}_1^3+6{\sigma }_1^2{a}_1^2+4{\alpha }_2{\overline{a}}_2{a}_2^3\\ & +6{\sigma }_2^2{a}_2^2+4{\alpha }_3{\overline{a}}_3{a}_3^3+6{\sigma }_3^2{a}_3^2-4{\alpha }_1{a}_1^4-4{\alpha }_2{a}_2^4-4{\alpha }_3{a}_3^4\\ & \le c_1(\frac{2}{3}{x}_1^{\frac{3}{2}}+\frac{1}{3}|a_1{|}^3)+c_2(\frac{2}{3}{x}_2^{\frac{3}{2}}+\frac{1}{3}|a_2{|}^3)+(\frac{2}{3}{x}_3^{\frac{3}{2}}+\frac{1}{3}|a_3{|}^3)\\ & -c_1{b}_{11}{x}_1^2-c_2{b}_{22}{x}_2^2-b_{33}{x}_3^2+(c_1{b}_{11}+c_2{b}_{21})x_1\\ & +(c_1{b}_{12}+c_2{b}_{22})x_2+(c_1{b}_{13}+c_2{b}_{23}+b_{33})x_3+c_1|a_1|+c_2|a_2|+|a_3|\\ & +4{\alpha }_1{\overline{a}}_1{a}_1^3+6{\sigma }_1^2{a}_1^2+4{\alpha }_2{\overline{a}}_2{a}_2^3+6{\sigma }_2^2{a}_2^2\\ & +4{\alpha }_3{\overline{a}}_3{a}_3^3+6{\sigma }_3^2{a}_3^2-4{\alpha }_1{a}_1^4-4{\alpha }_2{a}_2^4-4{\alpha }_3{a}_3^4+(b_{31}-c_1{b}_{13})x_1{x}_3\\ & +(b_{32}-c_2{b}_{23})x_2{x}_3\\ & \le \underset{x_1\in (0,\mathrm{\infty })}{sup}\{-c_1{b}_{11}{x}_1^2+\frac{2c_1}{3}{x}_1^{\frac{3}{2}}+(c_1{b}_{11}+c_2{b}_{21})x_1\}\\ & +\underset{x_2\in (0,\mathrm{\infty })}{sup}\{-c_2{b}_{22}{x}_2^2+\frac{2c_2}{3}{x}_2^{\frac{3}{2}}+(c_1{b}_{12}+c_2{b}_{22})x_2\}\\ & +\underset{x_3\in (0,\mathrm{\infty })}{sup}\{-b_{33}{x}_3^2+\frac{2}{3}{x}_3^{\frac{3}{2}}+(c_1{b}_{13}+c_2{b}_{23}+b_{33})x_3\}\\ & +\underset{a_1\in (-\mathrm{\infty },\mathrm{\infty })}{sup}\{-4{\alpha }_1{a}_1^4+4{\alpha }_1{\overline{a}}_1{a}_1^3+\frac{c_1}{3}|a_1{|}^3+6{\sigma }_1^2{a}_1^2+\phantom{\frac{c_1}{3}}{c}_1|a_1|\}\\ & +\underset{a_2\in (-\mathrm{\infty },\mathrm{\infty })}{sup}\{-4{\alpha }_2{a}_2^4+4{\alpha }_2{\overline{a}}_2{a}_2^3+\frac{c_2}{3}|a_2{|}^3+6{\sigma }_2^2{a}_2^2+c_2|a_2|\}\\ & +\underset{a_3\in (-\mathrm{\infty },\mathrm{\infty })}{sup}\{-4{\alpha }_3{a}_3^4+4{\alpha }_3{\overline{a}}_3{a}_3^3+\frac{1}{3}|a_3{|}^3+6{\sigma }_3^2{a}_3^2+|a_3|\}\\ & +(b_{31}-c_1{b}_{13})x_1{x}_3+(b_{32}-c_2{b}_{23})x_2{x}_3,\end{array}$where in the second inequality, we have used the Young inequality $a_1{x}_1\le \frac{2}{3}{x}_1^{\frac{3}{2}}+\frac{1}{3}|a_1{|}^3,a_2{x}_2\le \frac{2}{3}{x}_2^{\frac{3}{2}}+\frac{1}{3}|a_2{|}^3\mathrm{and}-a_3{x}_3\le \frac{2}{3}{x}_3^{\frac{3}{2}}+\frac{1}{3}|a_3{|}^3.$Choose $c_1=b_{31}/b_{13}$, $c_2=b_{32}/b_{23}$ such that $b_{31}-c_1{b}_{13}=0$ and $b_{32}-c_2{b}_{23}=0$, then we have $\begin{array}{rl}LV& \le \underset{x_1\in (0,\mathrm{\infty })}{sup}\{-c_1{b}_{11}{x}_1^2+\frac{2c_1}{3}{x}_1^{\frac{3}{2}}+(c_1{b}_{11}+c_2{b}_{21})x_1\}\\ & +\underset{x_2\in (0,\mathrm{\infty })}{sup}\{-c_2{b}_{22}{x}_2^2+\frac{2c_2}{3}{x}_2^{\frac{3}{2}}+(c_1{b}_{12}+c_2{b}_{22})x_2\}\\ & +\underset{x_3\in (0,\mathrm{\infty })}{sup}\{-b_{33}{x}_3^2+\frac{2}{3}{x}_3^{\frac{3}{2}}+(c_1{b}_{13}+c_2{b}_{23}+b_{33})x_3\}\\ & +\underset{a_1\in (-\mathrm{\infty },\mathrm{\infty })}{sup}\{-4{\alpha }_1{a}_1^4+4{\alpha }_1{\overline{a}}_1{a}_1^3+\frac{c_1}{3}|a_1{|}^3+6{\sigma }_1^2{a}_1^2+\phantom{\frac{c_1}{3}}{c}_1|a_1|\}\\ & +\underset{a_2\in (-\mathrm{\infty },\mathrm{\infty })}{sup}\{-4{\alpha }_2{a}_2^4+4{\alpha }_2{\overline{a}}_2{a}_2^3+\frac{c_2}{3}|a_2{|}^3+6{\sigma }_2^2{a}_2^2+c_2|a_2|\}\\ & +\underset{a_3\in (-\mathrm{\infty },\mathrm{\infty })}{sup}\{-4{\alpha }_3{a}_3^4+4{\alpha }_3{\overline{a}}_3{a}_3^3+\frac{1}{3}|a_3{|}^3+6{\sigma }_3^2{a}_3^2+|a_3|\}\\ & :=K_1.\end{array}$Here $K_1$ is a positive constant which is independent of the variables $x_1$, $x_2$, $x_3$, $a_1$, $a_2$ and $a_3$. The rest of the proof is similar to that of Zhang and Yuan [31] and so it is omitted here. The proof is confirmed.

3. The pth moment boundedness

In this section, the boundedness of the pth moment of solutions to system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) is shown by establishing some Lyapunov functions.

Theorem 3.1

Given initial value $(x_1(0),x_2(0),x_3(0),a_1(0),a_2(0),a_3(0){)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$ and any integer $p\ge 3$, there exists a positive constant $M(p)$ such that the solution $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) has the property $\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}[\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}+\frac{a_1^{2p}(t)}{2p}+\frac{a_2^{2p}(t)}{2p}+\frac{a_3^{2p}(t)}{2p}]\le M(p).$

Proof.

For any integer $p\ge 3$, define $\begin{array}{rl}V(x,a)& =\frac{{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^p}{p}+\frac{a_1^{2p}}{2p}+\frac{a_2^{2p}}{2p}+\frac{a_3^{2p}}{2p}\\ & =V_1(x)+V_2(a),x=(x_1,x_2,x_3{)}^T\in {\mathbb{R}}_{+}^3,a=(a_1,a_2,a_3{)}^T\in {\mathbb{R}}^3,\end{array}$where $V_1(x):=\frac{{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^p}{p},V_2(a):=\frac{a_1^{2p}}{2p}+\frac{a_2^{2p}}{2p}+\frac{a_3^{2p}}{2p}.$Applying Itô's formula [22] to $V_1(x)$, we have $\begin{array}{rl}L{V}_1& ={(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}\{\frac{b_{31}}{b_{13}}{x}_1(a_1-b_{11}{x}_1-b_{12}{x}_2-b_{13}{x}_3)\\ & +\frac{b_{32}}{b_{23}}{x}_2(a_2-b_{21}{x}_1-b_{22}{x}_2-b_{23}{x}_3)\\ & +\phantom{\frac{b_{31}}{b_{13}}}{x}_3(-a_3+b_{31}{x}_1+b_{32}{x}_2-b_{33}{x}_3)\}\\ & ={(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(-\frac{b_{11}{b}_{31}}{b_{13}}{x}_1^2-\frac{b_{22}{b}_{32}}{b_{23}}{x}_2^2-b_{33}{x}_3^2)\\ & +{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2-a_3{x}_3).\end{array}$Let $A=\left(\begin{array}{ccc}{\displaystyle \frac{b_{11}{b}_{31}}{b_{13}}}& 0& 0\\ 0& {\displaystyle \frac{b_{22}{b}_{32}}{b_{23}}}& 0\\ 0& 0& b_{33}\end{array}\right),$then A is positive definite. Let ${\lambda }_0$ be the minimum eigenvalue of the matrix A, i.e.${\lambda }_0=min\{b_{11}{b}_{31}/b_{13},b_{22}{b}_{32}/b_{23},b_{33}\}$, then we get $\left(x_1{x}_2{x}_3\right)A\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right)\ge {\lambda }_0(x_1^2+x_2^2+x_3^2).$In addition, we have ${(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^2\le (\frac{b_{31}^2}{b_{13}^2}+\frac{b_{32}^2}{b_{23}^2}+1)(x_1^2+x_2^2+x_3^2).$Hence (8) $\begin{array}{rl}L{V}_1& \le -{\lambda }_0{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(x_1^2+x_2^2+x_3^2)\\ & +{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2-a_3{x}_3)\\ & \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p+1}\\ & +{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2\phantom{\frac{b_{32}}{b_{23}}}-a_3{x}_3).\end{array}$(8) Similarly, applying Itô's formula [22] to $V_2(r)$, we obtain (9) $\begin{array}{rl}L{V}_2& ={\alpha }_1{a}_1^{2p-1}({\overline{a}}_1-a_1)+{\alpha }_2{a}_2^{2p-1}({\overline{a}}_2-a_2)+{\alpha }_3{a}_3^{2p-1}({\overline{a}}_3-a_3)\\ & +\frac{2p-1}{2}({\sigma }_1^2{a}_1^{2p-2}+{\sigma }_2^2{a}_2^{2p-2}+{\sigma }_3^2{a}_3^{2p-2})\\ & \le -\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2,\end{array}$(9) where $\begin{array}{rl}{K}_2& :=\underset{(a_1,a_2,a_3{)}^T\in {\mathbb{R}}^3}{sup}\{-\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+{\alpha }_1{\overline{a}}_1{a}_1^{2p-1}+{\alpha }_2{\overline{a}}_2{a}_2^{2p-1}\\ & +{\alpha }_3{\overline{a}}_3{a}_3^{2p-1}+\frac{2p-1}{2}{\sigma }_1^2{a}_1^{2p-2}+\phantom{\frac{2p-1}{2}}\frac{2p-1}{2}{\sigma }_2^2{a}_2^{2p-2}+\frac{2p-1}{2}{\sigma }_3^2{a}_3^{2p-2}\}<\mathrm{\infty }.\end{array}$From (8(8) $\begin{array}{rl}L{V}_1& \le -{\lambda }_0{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(x_1^2+x_2^2+x_3^2)\\ & +{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2-a_3{x}_3)\\ & \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p+1}\\ & +{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2\phantom{\frac{b_{32}}{b_{23}}}-a_3{x}_3).\end{array}$(8) ) and (9(9) $\begin{array}{rl}L{V}_2& ={\alpha }_1{a}_1^{2p-1}({\overline{a}}_1-a_1)+{\alpha }_2{a}_2^{2p-1}({\overline{a}}_2-a_2)+{\alpha }_3{a}_3^{2p-1}({\overline{a}}_3-a_3)\\ & +\frac{2p-1}{2}({\sigma }_1^2{a}_1^{2p-2}+{\sigma }_2^2{a}_2^{2p-2}+{\sigma }_3^2{a}_3^{2p-2})\\ & \le -\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2,\end{array}$(9) ), it follows that (10) $\begin{array}{rl}LV& \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p+1}\\ & +{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2-a_3{x}_3)\\ & -\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2\\ & =-\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{V}_1^{\frac{p+1}{p}}+p{V}_1\frac{\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2-a_3{x}_3}{\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3}\\ & -\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2\\ & \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{V}_1^{\frac{p+1}{p}}+p(|a_1|+|a_2|+|a_3|)V_1\\ & -\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2\\ & \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{V}_1^{\frac{p+1}{p}}+p(\frac{|a_1{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}})\\ & +p(\frac{|a_2{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}})\\ & +p(\frac{|a_3{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}})-\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2\\ & \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}-\frac{{\alpha }_1}{4}{a}_1^{2p}-\frac{{\alpha }_2}{4}{a}_2^{2p}-\frac{{\alpha }_3}{4}{a}_3^{2p}+K_3,\end{array}$(10) where in the second inequality, we have used the Young inequality $\begin{array}{rl}& |a_1|V_1\le \frac{|a_1{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}},|a_2|V_1\le \frac{|a_2{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}}\mathrm{and}\\ & |a_3|V_1\le \frac{|a_3{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}}\end{array}$and $\begin{array}{rl}{K}_3& :=\underset{(x_1,x_2,x_3,a_1,a_2,a_3{)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3}{sup}\{-\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}-\frac{{\alpha }_1}{4}{a}_1^{2p}-\frac{{\alpha }_2}{4}{a}_2^{2p}\\ & -\frac{{\alpha }_3}{4}{a}_3^{2p}+\frac{3p(p+\frac{1}{2})}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}}\\ & +\phantom{\frac{2}{2p+1}}\frac{p}{p+\frac{3}{2}}|a_1{|}^{p+\frac{3}{2}}+\frac{p}{p+\frac{3}{2}}|a_2{|}^{p+\frac{3}{2}}+\frac{p}{p+\frac{3}{2}}|a_3{|}^{p+\frac{3}{2}}+K_2\}<\mathrm{\infty }.\end{array}$Hence, for any $0<k<(p/2)min\{{\alpha }_1,{\alpha }_2,{\alpha }_3\}$, by means of (10(10) $\begin{array}{rl}LV& \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p+1}\\ & +{(\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3)}^{p-1}(\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2-a_3{x}_3)\\ & -\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2\\ & =-\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{V}_1^{\frac{p+1}{p}}+p{V}_1\frac{\frac{b_{31}}{b_{13}}{a}_1{x}_1+\frac{b_{32}}{b_{23}}{a}_2{x}_2-a_3{x}_3}{\frac{b_{31}}{b_{13}}{x}_1+\frac{b_{32}}{b_{23}}{x}_2+x_3}\\ & -\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2\\ & \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{V}_1^{\frac{p+1}{p}}+p(|a_1|+|a_2|+|a_3|)V_1\\ & -\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2\\ & \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2}{V}_1^{\frac{p+1}{p}}+p(\frac{|a_1{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}})\\ & +p(\frac{|a_2{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}})\\ & +p(\frac{|a_3{|}^{p+\frac{3}{2}}}{p+\frac{3}{2}}+\frac{p+\frac{1}{2}}{p+\frac{3}{2}}{V}_1^{1+\frac{2}{2p+1}})-\frac{{\alpha }_1}{2}{a}_1^{2p}-\frac{{\alpha }_2}{2}{a}_2^{2p}-\frac{{\alpha }_3}{2}{a}_3^{2p}+K_2\\ & \le -\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}-\frac{{\alpha }_1}{4}{a}_1^{2p}-\frac{{\alpha }_2}{4}{a}_2^{2p}-\frac{{\alpha }_3}{4}{a}_3^{2p}+K_3,\end{array}$(10) ), we obtain (11) $\begin{array}{rl}L({\mathrm{e}}^{kt}V(x(t),a(t)))& ={\mathrm{e}}^{kt}(kV(x(t),a(t))+LV(x(t),a(t)))\\ & \le {\mathrm{e}}^{kt}(-\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}-\frac{{\alpha }_1}{4}{a}_1^{2p}-\frac{{\alpha }_2}{4}{a}_2^{2p}\\ & -\frac{{\alpha }_3}{4}{a}_3^{2p}+k{V}_1+\frac{k}{2p}{a}_1^{2p}+\frac{k}{2p}{a}_2^{2p}\frac{k}{2p}{a}_3^{2p}+K_3\phantom{\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}})\\ & \le H(p){\mathrm{e}}^{kt},\end{array}$(11) where $\begin{array}{rl}H(p)& :=\underset{(x_1,x_2,x_3,a_1,a_2,a_3{)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3}{sup}\{-\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}-\frac{{\alpha }_1}{4}{a}_1^{2p}\\ & -\frac{{\alpha }_2}{4}{a}_2^{2p}-\frac{{\alpha }_3}{4}{a}_3^{2p}+k{V}_1+\frac{k}{2p}{a}_1^{2p}+\phantom{\frac{k}{2p}}\frac{k}{2p}{a}_2^{2p}+\frac{k}{2p}{a}_3^{2p}+K_3\phantom{\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}}\}<\mathrm{\infty }.\end{array}$Furthermore, according to (11(11) $\begin{array}{rl}L({\mathrm{e}}^{kt}V(x(t),a(t)))& ={\mathrm{e}}^{kt}(kV(x(t),a(t))+LV(x(t),a(t)))\\ & \le {\mathrm{e}}^{kt}(-\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}-\frac{{\alpha }_1}{4}{a}_1^{2p}-\frac{{\alpha }_2}{4}{a}_2^{2p}\\ & -\frac{{\alpha }_3}{4}{a}_3^{2p}+k{V}_1+\frac{k}{2p}{a}_1^{2p}+\frac{k}{2p}{a}_2^{2p}\frac{k}{2p}{a}_3^{2p}+K_3\phantom{\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}})\\ & \le H(p){\mathrm{e}}^{kt},\end{array}$(11) ), we have (12) $\begin{array}{rl}d[{\mathrm{e}}^{kt}V(x,a)]& ={\mathrm{e}}^{kt}kV(x,a)\mathrm{d}t+{\mathrm{e}}^{kt}dV(x,a)\\ & \le H(p){\mathrm{e}}^{kt}\mathrm{d}t+{\sigma }_1{\mathrm{e}}^{kt}{a}_1^{2p-1}\mathrm{d}{B}_1(t)+{\sigma }_2{\mathrm{e}}^{kt}{a}_2^{2p-1}\mathrm{d}{B}_2(t)+{\sigma }_3{\mathrm{e}}^{kt}{a}_3^{2p-1}\mathrm{d}{B}_3(t).\end{array}$(12) By integrating (12(12) $\begin{array}{rl}d[{\mathrm{e}}^{kt}V(x,a)]& ={\mathrm{e}}^{kt}kV(x,a)\mathrm{d}t+{\mathrm{e}}^{kt}dV(x,a)\\ & \le H(p){\mathrm{e}}^{kt}\mathrm{d}t+{\sigma }_1{\mathrm{e}}^{kt}{a}_1^{2p-1}\mathrm{d}{B}_1(t)+{\sigma }_2{\mathrm{e}}^{kt}{a}_2^{2p-1}\mathrm{d}{B}_2(t)+{\sigma }_3{\mathrm{e}}^{kt}{a}_3^{2p-1}\mathrm{d}{B}_3(t).\end{array}$(12) ) from 0 to t and then taking expectations on both sides, one has $\mathbb{E}V(x(t),a(t))\le {\mathrm{e}}^{-kt}V(x(0),a(0))+\frac{H(p)}{k}(1-{\mathrm{e}}^{-kt}).$Letting $t\to \mathrm{\infty }$ leads to that $\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}V(x(t),a(t))\le \frac{H(p)}{k}:=M(p).$The proof is confirmed.

Remark 3.1

In view of Theorem 3.1, it is easy to obtain that there exists a positive constant $S(p)$ such that $\begin{array}{rl}& \underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}(x_1^p(t))\le S(p),\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}(x_2^p(t))\le S(p),\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}(x_3^p(t))\le S(p),\\ & \underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}(a_1^{2p}(t))\le S(p),\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}(a_2^{2p}(t))\le S(p),\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}(a_3^{2p}(t))\le S(p).\end{array}$

4. Asymptotic pathwise estimation

In this section, we will verify pathwise properties of the solution $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ).

Theorem 4.1

Given initial value $(x_1(0),x_2(0),x_3(0),a_1(0),a_2(0),a_3(0){)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$, the solution $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) has the property $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln (\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t))}{\ln t}\le \frac{1}{p}\mathrm{a}.\mathrm{s}.,\mathrm{for}\mathrm{any}\mathrm{integer}p\ge 3.$In addition $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_1(t)}{t}\le 0,\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_2(t)}{t}\le 0\mathrm{and}\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_3(t)}{t}\le 0\mathrm{a}.\mathrm{s}.$

Proof.

According to Theorem 3.1, one can see that there is a positive constant $M_1(p)$ such that $\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}\left(\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}\right)\le M_1(p)\mathrm{for}\mathrm{any}\mathrm{integer}p\ge 3,$which together with the continuity of $x_1(t)$, $x_2(t)$ and $x_3(t)$ implies that there is a positive constant $M_2(p)$ such that (13) $\mathbb{E}\left(\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}\right)\le M_2(p)\mathrm{for}\mathrm{all}t\in [0,\mathrm{\infty }).$(13) For sufficiently small $\delta >0$, $k=1,2,\dots$, in view of (13(13) $\mathbb{E}\left(\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}\right)\le M_2(p)\mathrm{for}\mathrm{all}t\in [0,\mathrm{\infty }).$(13) ), we get $\begin{array}{rl}& \mathbb{E}\left[\underset{k\delta \le t\le (k+1)\delta }{sup}\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}\right]\\ & \le \mathbb{E}\left[\frac{(\frac{b_{31}}{b_{13}}{x}_1(k\delta )+\frac{b_{32}}{b_{23}}{x}_2(k\delta )+x_3(k\delta ){)}^p}{p}\right]+I\le M_2(p)+I,\end{array}$where $\begin{array}{rl}I& :=\mathbb{E}[\underset{k\delta \le t\le (k+1)\delta }{sup}|{\int }_{k\delta }^t(-\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}(s)-\frac{{\alpha }_1}{4}{a}_1^{2p}(s)-\frac{{\alpha }_2}{4}{a}_2^{2p}(s)\\ & -\frac{{\alpha }_3}{4}{a}_3^{2p}(s)+K_3\phantom{\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}})\mathrm{d}s|]\\ & \le \mathbb{E}[\underset{k\delta \le t\le (k+1)\delta }{sup}({\int }_{k\delta }^t\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{V}_1^{\frac{p+1}{p}}(s)\mathrm{d}s\\ & +\frac{{\alpha }_1}{4}{\int }_{k\delta }^t|a_1^{2p}(s)|\mathrm{d}s+\frac{{\alpha }_2}{4}{\int }_{k\delta }^t|a_2^{2p}(s)|\mathrm{d}s+{\int }_{k\delta }^t|K_3|\mathrm{d}s+\frac{{\alpha }_3}{4}{\int }_{k\delta }^t|a_3^{2p}(s)|\mathrm{d}s\phantom{\frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}}]\\ & \le \frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}{\int }_{k\delta }^{(k+1)\delta }\mathbb{E}{V}_1^{\frac{p+1}{p}}(s)\mathrm{d}s+\frac{{\alpha }_1}{4}{\int }_{k\delta }^{(k+1)\delta }\mathbb{E}|a_1^{2p}(s)|\mathrm{d}s\\ & +\frac{{\alpha }_2}{4}{\int }_{k\delta }^{(k+1)\delta }\mathbb{E}|a_2^{2p}(s)|\mathrm{d}s+\frac{{\alpha }_3}{4}{\int }_{k\delta }^{(k+1)\delta }\mathbb{E}|a_3^{2p}(s)|\mathrm{d}s+{\int }_{k\delta }^{(k+1)\delta }|K_3|\mathrm{d}s\\ & \le \frac{{\lambda }_0{b}_{13}^2{b}_{23}^2{p}^{\frac{p+1}{p}}}{2(b_{13}^2{b}_{23}^2+b_{13}^2{b}_{32}^2+b_{23}^2{b}_{31}^2)}M(p)+\frac{{\alpha }_1\delta }{4}S(p)+\frac{{\alpha }_2\delta }{4}S(p)+\frac{{\alpha }_3\delta }{4}S(p)+|K_3|\delta \\ & :=M_3(p).\end{array}$So $\begin{array}{rl}& \mathbb{E}\left[\underset{k\delta \le t\le (k+1)\delta }{sup}\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}\right]\\ & \le \mathbb{E}\left[\frac{(\frac{b_{31}}{b_{13}}{x}_1(k\delta )+\frac{b_{32}}{b_{23}}{x}_2(k\delta )+x_3(k\delta ){)}^p}{p}\right]+M_3(p)\\ \\ & \le M_2(p)+M_3(p)\\ \\ & :=M_4(p).\end{array}$Let ${ϵ}_x>0$ be arbitrary. By Chebyshev's inequality [22], we obtain $\begin{array}{rl}& \mathbb{P}\{\underset{k\delta \le t\le (k+1)\delta }{sup}\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}>(k\delta {)}^{1+{ϵ}_x}\}\\ & \le \frac{\mathbb{E}\left[\underset{k\delta \le t\le (k+1)\delta }{sup}\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}\right]}{(k\delta {)}^{1+{ϵ}_x}}\\ & \le \frac{M_4(p)}{(k\delta {)}^{1+{ϵ}_x}},k=1,2,\dots \end{array}$By Borel–Cantelli Lemma [22], one obtains that for almost all $\omega \in \mathrm{\Omega }$, (14) $\underset{k\delta \le t\le (k+1)\delta }{sup}\frac{{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t))}^p}{p}\le (k\delta {)}^{1+{ϵ}_x},$(14) holds for all but finitely many k. As a result, there exists a $k_0(\omega )$ for almost all $\omega \in \mathrm{\Omega }$, for which (14(14) $\underset{k\delta \le t\le (k+1)\delta }{sup}\frac{{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t))}^p}{p}\le (k\delta {)}^{1+{ϵ}_x},$(14) ) holds whenever $k\ge k_0$. Thus, for almost all $\omega \in \mathrm{\Omega }$, if $k\ge k_0$ and $k\delta \le t\le (k+1)\delta$, $\frac{\ln \left(\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}\right)}{\ln t}\le \frac{(1+{ϵ}_x)\ln (k\delta )}{\ln (k\delta )}=1+{ϵ}_x\mathrm{a}.\mathrm{s}.$Therefore $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln \left(\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}\right)}{\ln t}\le 1+{ϵ}_x\mathrm{a}.\mathrm{s}.$Letting ${ϵ}_x\to 0^{+}$ gives that $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln {(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t))}^p}{\ln t}\le 1\mathrm{a}.\mathrm{s}.$and hence $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln (\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t))}{\ln t}\le \frac{1}{p}\mathrm{a}.\mathrm{s}.$Since the integer p is arbitrary and greater than 0, it can be seen that $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_1(t)}{t}\le 0,\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_2(t)}{t}\le 0\mathrm{and}\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_3(t)}{t}\le 0\mathrm{a}.\mathrm{s}.$This completes the proof.

5. Existence of a stationary distribution

In this section, sufficient conditions which are used to ensure the existence of a stationary distribution of the stochastic system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) are presented.

Let $X(t)$ be a homogeneous Markov process defined in ${\mathbb{R}}^d$ which is described by the following stochastic differential equation: (15) $\mathrm{d}X(t)=b(X(t))\mathrm{d}t+\sum _{r=1}^k{\sigma }_r(X(t))\mathrm{d}{B}_r(t).$(15)

Lemma 5.1

[11]

Assume that the vectors $b(X)$, ${\sigma }_1(X),\dots ,{\sigma }_k(X)$ $(t\ge t_0,X\in {\mathbb{R}}^d)$ are continuous functions with respect to X and satisfy the following conditions:

$(A_1)$ There is a constant B such that (16) $\begin{array}{rl}& |b(X_1)-b(X_2)|+\sum _{r=1}^k|{\sigma }_r(X_1)-{\sigma }_r(X_2)|\le B|X_1-X_2|,\\ & |b(X)|+\sum _{r=1}^k|{\sigma }_r(X)|\le B(1+|X|).\end{array}$(16) $(A_2)$ There exists a nonnegative $C^2$-function $V(X)$ in ${\mathbb{R}}^d$ such that $LV(X)\le -1$ outside some compact set. Then system (15(15) $\mathrm{d}X(t)=b(X(t))\mathrm{d}t+\sum _{r=1}^k{\sigma }_r(X(t))\mathrm{d}{B}_r(t).$(15) ) admits at least one stationary distribution on ${\mathbb{R}}^d$.

Remark 5.1

We can use the global existence of solutions of (15(15) $\mathrm{d}X(t)=b(X(t))\mathrm{d}t+\sum _{r=1}^k{\sigma }_r(X(t))\mathrm{d}{B}_r(t).$(15) ) to replace the conditions (16(16) $\begin{array}{rl}& |b(X_1)-b(X_2)|+\sum _{r=1}^k|{\sigma }_r(X_1)-{\sigma }_r(X_2)|\le B|X_1-X_2|,\\ & |b(X)|+\sum _{r=1}^k|{\sigma }_r(X)|\le B(1+|X|).\end{array}$(16) ) in Lemma 5.1 and the conclusion can be found in Remark 5 of Xu et al. [30].

Theorem 5.1

Let Assumption (A) hold. If the following conditions hold $\begin{array}{rl}& b_{11}>\frac{{\sigma }_1}{2},b_{22}>\frac{{\sigma }_2}{2},b_{33}>\frac{{\sigma }_3}{2},\\ & 4b_{13}{b}_{23}{b}_{31}{b}_{32}(b_{11}-\frac{{\sigma }_1}{2})(b_{22}-\frac{{\sigma }_2}{2})>(b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}{)}^2,\\ & \omega <min\{\frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})-\frac{ϵ}{2}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}}){\overline{x}}_1^2,\\ & +\frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})-\frac{1}{2ϵ}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}}){\overline{x}}_2^2,(b_{33}-\frac{{\sigma }_3}{2}){\overline{x}}_3^2\},\end{array}$and $ϵ>0$ is a sufficiently small constant satisfying $\frac{b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}}{2b_{13}{b}_{32}(b_{22}-\frac{{\sigma }_2}{2})}<ϵ<\frac{2b_{23}{b}_{31}(b_{11}-\frac{{\sigma }_1}{2})}{b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}},$then system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) has at least one stationary distribution $\pi (\cdot )$ on ${\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$, where $\omega =b_{31}{\sigma }_1/(2b_{13}{\alpha }_1)+b_{32}{\sigma }_2/(2b_{23}{\alpha }_2)+{\sigma }_3/(2{\alpha }_3)$ and $({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3)$ is the positive equilibrium of system (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ).

Proof.

By Theorem 2.1, we have proved that there is a global solution to system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ). Hence, in order to complete the proof of Theorem 5.1, we only need to prove the condition $(A_2)$. That is to say, we need to construct a nonnegative $C^2$-function $V(x_1,x_2,x_3,a_1,a_2,a_3)$ and a compact set $D\subset {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$ such that $LV\le -1$ for all $(x_1,x_2,x_3,a_1,a_2,a_3)\in ({\mathbb{R}}_{+}^3\times {\mathbb{R}}^3)\setminus D$.

Define a nonnegative $C^2$-function $V:{\mathbb{R}}_{+}^3\times {\mathbb{R}}^3\to \mathbb{R}$ as follows: $\begin{array}{rl}V(x_1,x_2,x_3,a_1,a_2,a_3)& =\frac{b_{31}}{b_{13}}(x_1-{\overline{x}}_1-{\overline{x}}_1\ln \frac{x_1}{{\overline{x}}_1})+\frac{b_{32}}{b_{23}}(x_2-{\overline{x}}_2-{\overline{x}}_2\ln \frac{x_2}{{\overline{x}}_2})\\ & +(x_3-{\overline{x}}_3-{\overline{x}}_3\ln \frac{x_3}{{\overline{x}}_3})+\frac{b_{31}}{2b_{13}{\alpha }_1{\sigma }_1}(a_1-{\overline{a}}_1{)}^2\\ & +\frac{b_{32}}{2b_{23}{\alpha }_2{\sigma }_2}(a_2-{\overline{a}}_2{)}^2+\frac{1}{2{\alpha }_3{\sigma }_3}(a_3-{\overline{a}}_3{)}^2,\end{array}$where $({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3)$ is the positive equilibrium of system (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ) which satisfies the following algebraic equations (17) $\{\begin{array}{l}{\overline{a}}_1-b_{11}{\overline{x}}_1-b_{12}{\overline{x}}_2-b_{13}{\overline{x}}_3=0,\\ {\overline{a}}_2-b_{21}{\overline{x}}_1-b_{22}{\overline{x}}_2-b_{23}{\overline{x}}_3=0,\\ -{\overline{a}}_3+b_{31}{\overline{x}}_1+b_{32}{\overline{x}}_2-b_{33}{\overline{x}}_3=0.\end{array}$(17) For convenience, denote by $\begin{array}{rl}{V}_1(x_1)& =x_1-{\overline{x}}_1-{\overline{x}}_1\ln \frac{x_1}{{\overline{x}}_1},V_2(x_2)=x_2-{\overline{x}}_2-{\overline{x}}_2\ln \frac{x_2}{{\overline{x}}_2},\\ V_3(x_3)& =x_3-{\overline{x}}_3-{\overline{x}}_3\ln \frac{x_3}{{\overline{x}}_3},\\ V_4(a_1)& =\frac{b_{31}}{2b_{13}{\alpha }_1{\sigma }_1}(a_1-{\overline{a}}_1{)}^2,V_5(a_2)=\frac{b_{32}}{2b_{23}{\alpha }_2{\sigma }_2}(a_2-{\overline{a}}_2{)}^2,\\ V_6(a_3)& =\frac{1}{2{\alpha }_3{\sigma }_3}(a_3-{\overline{a}}_3{)}^2.\end{array}$Applying Itô's formula [22] to $V_1$ and using (17(17) $\{\begin{array}{l}{\overline{a}}_1-b_{11}{\overline{x}}_1-b_{12}{\overline{x}}_2-b_{13}{\overline{x}}_3=0,\\ {\overline{a}}_2-b_{21}{\overline{x}}_1-b_{22}{\overline{x}}_2-b_{23}{\overline{x}}_3=0,\\ -{\overline{a}}_3+b_{31}{\overline{x}}_1+b_{32}{\overline{x}}_2-b_{33}{\overline{x}}_3=0.\end{array}$(17) ), we have (18) $\begin{array}{rl}L{V}_1& =(x_1-{\overline{x}}_1)({\overline{a}}_1-b_{11}{x}_1-b_{12}{x}_2-b_{13}{x}_3)+(x_1-{\overline{x}}_1)(a_1-{\overline{a}}_1)\\ & =(x_1-{\overline{x}}_1)(b_{11}{\overline{x}}_1+b_{12}{\overline{x}}_2+b_{13}{\overline{x}}_3-b_{11}{x}_1-b_{12}{x}_2-b_{13}{x}_3)+(x_1-{\overline{x}}_1)(a_1-{\overline{a}}_1)\\ & =-b_{11}(x_1-{\overline{x}}_1{)}^2-b_{12}(x_1-{\overline{x}}_1)(x_2-{\overline{x}}_2)-b_{13}(x_1-{\overline{x}}_1)(x_3-{\overline{x}}_3)\\ & +(x_1-{\overline{x}}_1)(a_1-{\overline{a}}_1)\\ & \le -(b_{11}-\frac{{\sigma }_1}{2})(x_1-{\overline{x}}_1{)}^2-b_{12}(x_1-{\overline{x}}_1)(x_2-{\overline{x}}_2)\\ & -b_{13}(x_1-{\overline{x}}_1)(x_3-{\overline{x}}_3)+\frac{(a_1-{\overline{a}}_1{)}^2}{2{\sigma }_1},\end{array}$(18) where in the above inequality, we have used the Young inequality $(x_1-{\overline{x}}_1)(a_1-{\overline{a}}_1)\le \frac{{\sigma }_1}{2}(x_1-{\overline{x}}_1{)}^2+\frac{1}{2{\sigma }_1}(a_1-{\overline{a}}_1{)}^2$and ${\sigma }_1\in (0,2b_{11})$ is a sufficiently small constant to ensure the coefficient in front of $(x_1-{\overline{x}}_1{)}^2$ is negative.

Similarly, we obtain (19) $\begin{array}{rl}L{V}_2& =(x_2-{\overline{x}}_2)({\overline{a}}_2-b_{21}{x}_1-b_{22}{x}_2-b_{23}{x}_3)+(x_2-{\overline{x}}_2)(a_2-{\overline{a}}_2)\\ & =(x_2-{\overline{x}}_2)(b_{21}{\overline{x}}_1+b_{22}{\overline{x}}_2+b_{23}{\overline{x}}_3-b_{21}{x}_1-b_{22}{x}_2-b_{23}{x}_3)\\ & +(x_2-{\overline{x}}_2)(a_2-{\overline{a}}_2)\\ & =-b_{22}(x_2-{\overline{x}}_2{)}^2-b_{21}(x_1-{\overline{x}}_1)(x_2-{\overline{x}}_2)-b_{23}(x_2-{\overline{x}}_2)(x_3-{\overline{x}}_3)\\ & +(x_2-{\overline{x}}_2)(a_2-{\overline{a}}_2)\\ & \le -(b_{22}-\frac{{\sigma }_2}{2})(x_2-{\overline{x}}_2{)}^2-b_{21}(x_1-{\overline{x}}_1)(x_2-{\overline{x}}_2)\\ & -b_{23}(x_2-{\overline{x}}_2)(x_3-{\overline{x}}_3)+\frac{(a_2-{\overline{a}}_2{)}^2}{2{\sigma }_2},\end{array}$(19) (20) $\begin{array}{rl}L{V}_3& =(x_3-{\overline{x}}_3)(-{\overline{a}}_3+b_{31}{x}_1+b_{32}{x}_2-b_{33}{x}_3)+(x_3-{\overline{x}}_3)(-a_3+{\overline{a}}_3)\\ & =(x_3-{\overline{x}}_3)(-b_{31}{\overline{x}}_1-b_{32}{\overline{x}}_2+b_{33}{\overline{x}}_3+b_{31}{x}_1+b_{32}{x}_2-b_{33}{x}_3)\\ & +(x_3-{\overline{x}}_3)(-a_3+{\overline{a}}_3)\\ & =-b_{33}(x_3-{\overline{x}}_3{)}^2+b_{31}(x_1-{\overline{x}}_1)(x_3-{\overline{x}}_3)+b_{32}(x_2-{\overline{x}}_2)(x_3-{\overline{x}}_3)\\ & +(x_3-{\overline{x}}_3)(-a_3+{\overline{a}}_3)\\ & \le -(b_{33}-\frac{{\sigma }_3}{2})(x_3-{\overline{x}}_3{)}^2+b_{31}(x_1-{\overline{x}}_1)(x_3-{\overline{x}}_3)+b_{32}(x_2-{\overline{x}}_2)(x_3-{\overline{x}}_3)\\ & +\frac{(a_3-{\overline{a}}_3{)}^2}{2{\sigma }_3},\end{array}$(20) where ${\sigma }_2\in (0,2b_{22})$, ${\sigma }_3\in (0,2b_{33})$ are sufficiently small constants to ensure the coefficients in front of $(x_2-{\overline{x}}_2{)}^2$, $(x_3-{\overline{x}}_3{)}^2$ are negative, respectively.

Moreover, we have (21) $\begin{array}{rl}L{V}_4& =\frac{b_{31}}{b_{13}{\sigma }_1}(a_1-{\overline{a}}_1)({\overline{a}}_1-a_1)+\frac{b_{31}{\sigma }_1}{2b_{13}{\alpha }_1}\\ & =-\frac{b_{31}}{b_{13}{\sigma }_1}(a_1-{\overline{a}}_1{)}^2+\frac{b_{31}{\sigma }_1}{2b_{13}{\alpha }_1},\end{array}$(21) (22) $\begin{array}{rl}L{V}_5& =\frac{b_{32}}{b_{23}{\sigma }_2}(a_2-{\overline{a}}_2)({\overline{a}}_2-a_2)+\frac{b_{32}{\sigma }_2}{2b_{23}{\alpha }_2}\\ & =-\frac{b_{32}}{b_{23}{\sigma }_2}(a_2-{\overline{a}}_2{)}^2+\frac{b_{32}{\sigma }_2}{2b_{23}{\alpha }_2},\end{array}$(22) and (23) $\begin{array}{rl}L{V}_6& =\frac{1}{{\sigma }_3}(a_3-{\overline{a}}_3)({\overline{a}}_3-a_3)+\frac{{\sigma }_3}{2{\alpha }_3}\\ & =-\frac{1}{{\sigma }_3}(a_3-{\overline{a}}_3{)}^2+\frac{{\sigma }_3}{2{\alpha }_3}.\end{array}$(23) Then according to (18(18) $\begin{array}{rl}L{V}_1& =(x_1-{\overline{x}}_1)({\overline{a}}_1-b_{11}{x}_1-b_{12}{x}_2-b_{13}{x}_3)+(x_1-{\overline{x}}_1)(a_1-{\overline{a}}_1)\\ & =(x_1-{\overline{x}}_1)(b_{11}{\overline{x}}_1+b_{12}{\overline{x}}_2+b_{13}{\overline{x}}_3-b_{11}{x}_1-b_{12}{x}_2-b_{13}{x}_3)+(x_1-{\overline{x}}_1)(a_1-{\overline{a}}_1)\\ & =-b_{11}(x_1-{\overline{x}}_1{)}^2-b_{12}(x_1-{\overline{x}}_1)(x_2-{\overline{x}}_2)-b_{13}(x_1-{\overline{x}}_1)(x_3-{\overline{x}}_3)\\ & +(x_1-{\overline{x}}_1)(a_1-{\overline{a}}_1)\\ & \le -(b_{11}-\frac{{\sigma }_1}{2})(x_1-{\overline{x}}_1{)}^2-b_{12}(x_1-{\overline{x}}_1)(x_2-{\overline{x}}_2)\\ & -b_{13}(x_1-{\overline{x}}_1)(x_3-{\overline{x}}_3)+\frac{(a_1-{\overline{a}}_1{)}^2}{2{\sigma }_1},\end{array}$(18) )–(23(23) $\begin{array}{rl}L{V}_6& =\frac{1}{{\sigma }_3}(a_3-{\overline{a}}_3)({\overline{a}}_3-a_3)+\frac{{\sigma }_3}{2{\alpha }_3}\\ & =-\frac{1}{{\sigma }_3}(a_3-{\overline{a}}_3{)}^2+\frac{{\sigma }_3}{2{\alpha }_3}.\end{array}$(23) ), we have $\begin{array}{rl}LV& \le -(b_{33}-\frac{{\sigma }_3}{2})(x_3-{\overline{x}}_3{)}^2-\frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})(x_1-{\overline{x}}_1{)}^2\\ & -\frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})(x_2-{\overline{x}}_2{)}^2-\frac{b_{31}}{2b_{13}{\sigma }_1}(a_1-{\overline{a}}_1{)}^2\\ & -\frac{b_{32}}{2b_{23}{\sigma }_2}(a_2-{\overline{a}}_2{)}^2-\frac{1}{2{\sigma }_3}(a_3-{\overline{a}}_3{)}^2+\frac{b_{31}{\sigma }_1}{2b_{13}{\alpha }_1}+\frac{b_{32}{\sigma }_2}{2b_{23}{\alpha }_2}+\frac{{\sigma }_3}{2{\alpha }_3}\\ & -(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})(x_1-{\overline{x}}_1)(x_2-{\overline{x}}_2)\\ & \le -(b_{33}-\frac{{\sigma }_3}{2})(x_3-{\overline{x}}_3{)}^2-\frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})(x_1-{\overline{x}}_1{)}^2\\ & -\frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})(x_2-{\overline{x}}_2{)}^2-\frac{b_{31}}{2b_{13}{\sigma }_1}(a_1-{\overline{a}}_1{)}^2\\ & -\frac{b_{32}}{2b_{23}{\sigma }_2}(a_2-{\overline{a}}_2{)}^2-\frac{1}{2{\sigma }_3}(a_3-{\overline{a}}_3{)}^2+\frac{b_{31}{\sigma }_1}{2b_{13}{\alpha }_1}\\ & +\frac{b_{32}{\sigma }_2}{2b_{23}{\alpha }_2}+\frac{{\sigma }_3}{2{\alpha }_3}+\frac{ϵ}{2}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})(x_1-{\overline{x}}_1{)}^2\\ & +\frac{1}{2ϵ}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})(x_2-{\overline{x}}_2{)}^2\\ & =-[\frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})-\frac{ϵ}{2}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})](x_1\\ & -{\overline{x}}_1{)}^2-[\frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})-\frac{1}{2ϵ}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})](x_2-{\overline{x}}_2{)}^2\\ & -(b_{33}-\frac{{\sigma }_3}{2})(x_3-{\overline{x}}_3{)}^2-\frac{b_{31}}{2b_{13}{\sigma }_1}(a_1-{\overline{a}}_1{)}^2-\frac{b_{32}}{2b_{23}{\sigma }_2}(a_2-{\overline{a}}_2{)}^2\\ & -\frac{1}{2{\sigma }_3}(a_3-{\overline{a}}_3{)}^2+\frac{b_{31}{\sigma }_1}{2b_{13}{\alpha }_1}+\frac{b_{32}{\sigma }_2}{2b_{23}{\alpha }_2}+\frac{{\sigma }_3}{2{\alpha }_3},\end{array}$where in the second inequality we have used the Young inequality $-(x_1-{\overline{x}}_1)(x_2-{\overline{x}}_2)\le \frac{ϵ}{2}(x_1-{\overline{x}}_1{)}^2+\frac{1}{2ϵ}(x_2-{\overline{x}}_2{)}^2$and $\frac{b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}}{2b_{13}{b}_{32}(b_{22}-\frac{{\sigma }_2}{2})}<ϵ<\frac{2b_{23}{b}_{31}(b_{11}-\frac{{\sigma }_1}{2})}{b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}}$such that $\begin{array}{rl}& \frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})-\frac{ϵ}{2}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})>0\mathrm{and}\\ & \frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})-\frac{1}{2ϵ}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})>0.\end{array}$Let $\omega =\frac{b_{31}{\sigma }_1}{2b_{13}{\alpha }_1}+\frac{b_{32}{\sigma }_2}{2b_{23}{\alpha }_2}+\frac{{\sigma }_3}{2{\alpha }_3},$and note that if the condition $\begin{array}{rl}\omega & <min\{[\frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})-\frac{ϵ}{2}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})]\\ & {\overline{x}}_1^2,[\frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})-\frac{1}{2ϵ}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})]{\overline{x}}_2^2,(b_{33}-\frac{{\sigma }_3}{2}){\overline{x}}_3^2\}\end{array}$holds, then the ellipsoid $\begin{array}{rl}& -[\frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})-\frac{ϵ}{2}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})](x_1-{\overline{x}}_1{)}^2\\ & -[\frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})-\frac{1}{2ϵ}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}})](x_2-{\overline{x}}_2{)}^2\\ & -(b_{33}-\frac{{\sigma }_3}{2})(x_3-{\overline{x}}_3{)}^2-\frac{b_{31}}{2b_{13}{\sigma }_1}(a_1-{\overline{a}}_1{)}^2\\ & -\frac{b_{32}}{2b_{23}{\sigma }_2}(a_2-{\overline{a}}_2{)}^2-\frac{1}{2{\sigma }_3}(a_3-{\overline{a}}_3{)}^2+\omega =0\end{array}$lies entirely in ${\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$. Thus, we can take $D_{ϵ}$ to be a neighbourhood of the ellipsoid with ${\overline{D}}_{ϵ}\subset {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$, where $D_{ϵ}:=(ϵ,1/ϵ)\times (ϵ,1/ϵ)\times (ϵ,1/ϵ)\times (-1/ϵ,1/ϵ)\times (-1/ϵ,1/ϵ)\times (-1/ϵ,1/ϵ)$ and ${\overline{D}}_{ϵ}$ denotes the compact closure of $D_{ϵ}$. So for $(x_1,x_2,x_3,a_1,a_2,a_3{)}^T\in ({\mathbb{R}}_{+}^3\times {\mathbb{R}}^3)\setminus D_{ϵ}$, $LV\le -1$, which implies that the condition $(A_2)$ in Lemma 5.1 holds. According to Lemma 5.1, we obtain that there is a solution $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) which is a stationary Markov process. This completes the proof.

6. Local density function for system (5)

In this section, we will get the specific form of covariance matrix in the probability density around the quasi-positive equilibrium of the stochastic system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ).

6.1. Equivalent transformation of system (5)

Firstly, we define a quasi-positive equilibrium $E^{\ast }=(x_1^{\ast },x_2^{\ast },x_3^{\ast },a_1^{\ast },a_2^{\ast },a_3^{\ast }{)}^T$ involved in stochasticity by the equations (24) $\{\begin{array}{l}{x}_1^{\ast }(a_1^{\ast }-b_{11}{x}_1^{\ast }-b_{12}{x}_2^{\ast }-b_{13}{x}_3^{\ast })=0,\\ x_2^{\ast }(a_2^{\ast }-b_{21}{x}_1^{\ast }-b_{22}{x}_2^{\ast }-b_{23}{x}_3^{\ast })=0,\\ x_3^{\ast }(-a_3^{\ast }+b_{31}{x}_1^{\ast }+b_{32}{x}_2^{\ast }-b_{33}{x}_3^{\ast })=0,\\ {\alpha }_1({\overline{a}}_1-a_1^{\ast })=0,\\ {\alpha }_2({\overline{a}}_2-a_2^{\ast })=0,\\ {\alpha }_3({\overline{a}}_3-a_3^{\ast })=0.\end{array}$(24) By solving Equation (24(24) $\{\begin{array}{l}{x}_1^{\ast }(a_1^{\ast }-b_{11}{x}_1^{\ast }-b_{12}{x}_2^{\ast }-b_{13}{x}_3^{\ast })=0,\\ x_2^{\ast }(a_2^{\ast }-b_{21}{x}_1^{\ast }-b_{22}{x}_2^{\ast }-b_{23}{x}_3^{\ast })=0,\\ x_3^{\ast }(-a_3^{\ast }+b_{31}{x}_1^{\ast }+b_{32}{x}_2^{\ast }-b_{33}{x}_3^{\ast })=0,\\ {\alpha }_1({\overline{a}}_1-a_1^{\ast })=0,\\ {\alpha }_2({\overline{a}}_2-a_2^{\ast })=0,\\ {\alpha }_3({\overline{a}}_3-a_3^{\ast })=0.\end{array}$(24) ), we obtain that if Assumption (A) holds, then $x_1^{\ast }={\overline{x}}_1>0,x_2^{\ast }={\overline{x}}_2>0,x_3^{\ast }={\overline{x}}_3>0,a_1^{\ast }={\overline{a}}_1,a_2^{\ast }={\overline{a}}_2,a_3^{\ast }={\overline{a}}_3,$where $({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3)$ is the positive equilibrium of system (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ).

Let $(u_1,u_2,u_3,u_4,u_5,u_6{)}^T=(x_1-x_1^{\ast },x_2-x_2^{\ast },x_3-x_3^{\ast },a_1-a_1^{\ast },a_2-a_2^{\ast },a_3-a_3^{\ast }{)}^T$. In view of Itô's integral and system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ), the corresponding linearised system of (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) around $E^{\ast }$ can be expressed as follows (25) $\{\begin{array}{l}\mathrm{d}{u}_1=(-q_{11}{u}_1-q_{12}{u}_2-q_{13}{u}_3+q_{14}{u}_4)\mathrm{d}t,\\ \mathrm{d}{u}_2=(-q_{21}{u}_1-q_{22}{u}_2-q_{23}{u}_3+q_{25}{u}_5)\mathrm{d}t,\\ \mathrm{d}{u}_3=(q_{31}{u}_1+q_{32}{u}_2-q_{33}{u}_3-q_{36}{u}_6)\mathrm{d}t,\\ \mathrm{d}{u}_4=-{\alpha }_1{u}_4\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{u}_5=-{\alpha }_2{u}_5\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{u}_6=-{\alpha }_3{u}_6\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t),\end{array}$(25) where $q_{11}=2b_{11}{x}_1^{\ast }+b_{12}{x}_2^{\ast }+b_{13}{x}_3^{\ast }-a_1^{\ast }=b_{11}{x}_1^{\ast }>0$, $q_{12}=b_{12}{x}_1^{\ast }>0$, $q_{13}=b_{13}{x}_1^{\ast }>0$, $q_{14}=x_1^{\ast }>0$, $q_{21}=b_{21}{x}_2^{\ast }>0$, $q_{22}=b_{21}{x}_1^{\ast }+2b_{22}{x}_2^{\ast }+b_{23}{x}_3^{\ast }-{\overline{a}}_2=b_{22}{x}_2^{\ast }>0$, $q_{23}=b_{23}{x}_2^{\ast }>0$, $q_{25}=x_2^{\ast }>0$, $q_{31}=b_{31}{x}_3^{\ast }>0$, $q_{32}=b_{32}{x}_3^{\ast }>0$, $q_{33}=2b_{33}{x}_3^{\ast }+{\overline{a}}_3-b_{31}{x}_1^{\ast }-b_{32}{x}_2^{\ast }=b_{33}{x}_3^{\ast }>0$, $q_{36}=x_3^{\ast }>0$.

6.2. Local density function for system (5)

In this subsection, we will obtain the specific form of covariance matrix in the probability density around the quasi-positive equilibrium of the stochastic system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ). To this end, we need to introduce a definition and a lemma.

Definition 6.1

[20]

The characteristic polynomial of the square matrix $A_n$ is defined as ${\phi }_{A_n}(\lambda )={\lambda }^n+a_1{\lambda }^{n-1}+\dots +a_{n-1}\lambda +a_n$, then $A_n$ is called a Hurwitz matrix if and only if $A_n$ has all negative real-part eigenvalues, i.e. $H_p=\left|\begin{array}{ccccc}{a}_1& a_3& a_5& \dots & a_{2p-1}\\ 1& a_2& a_4& \dots & a_{2p-2}\\ 0& a_1& a_3& \dots & a_{2p-3}\\ 0& 1& a_2& \dots & a_{2p-4}\\ \dots & \dots & \dots & \dots & \dots \\ 0& 0& 0& \dots & a_p\end{array}\right|>0,p=1,2,\dots ,n,$where the complementary definition $a_j=0$, j>n.

Lemma 6.1

[7]

For the three-dimensional real algebraic equation $H_0^2+A_0{\mathrm{\Sigma }}_0+{\mathrm{\Sigma }}_0{A}_0^T=0$, where $H_0=\mathrm{diag}(1,0,0)$, and $A_0=\left(\begin{array}{ccc}-{\tau }_1& -{\tau }_2& -{\tau }_3\\ 1& 0& 0\\ 0& 1& 0\end{array}\right).$If ${\mathrm{\Sigma }}_0$ is a three-dimensional symmetric matrix, ${\tau }_1>0,{\tau }_3>0\mathrm{and}{\tau }_1{\tau }_2-{\tau }_3>0,$then ${\mathrm{\Sigma }}_0$ is a positive definite matrix, where $\begin{array}{r}{\mathrm{\Sigma }}_0=\left(\begin{array}{ccc}{\displaystyle \frac{{\tau }_2}{2({\tau }_1{\tau }_2-{\tau }_3)}}& 0& -{\displaystyle \frac{1}{2({\tau }_1{\tau }_2-{\tau }_3)}}\\ 0& {\displaystyle \frac{1}{2({\tau }_1{\tau }_2-{\tau }_3)}}& 0\\ -{\displaystyle \frac{1}{2({\tau }_1{\tau }_2-{\tau }_3)}}& 0& {\displaystyle \frac{{\tau }_1}{2{\tau }_3({\tau }_1{\tau }_2-{\tau }_3)}}\end{array}\right).\end{array}$Here $A_0$ is called the standard $R_1$ matrix.

Theorem 6.1

Let $(u_1(t),u_2(t),u_3(t),u_4(t),u_5(t),u_6(t){)}^T$ be a solution of system (25(25) $\{\begin{array}{l}\mathrm{d}{u}_1=(-q_{11}{u}_1-q_{12}{u}_2-q_{13}{u}_3+q_{14}{u}_4)\mathrm{d}t,\\ \mathrm{d}{u}_2=(-q_{21}{u}_1-q_{22}{u}_2-q_{23}{u}_3+q_{25}{u}_5)\mathrm{d}t,\\ \mathrm{d}{u}_3=(q_{31}{u}_1+q_{32}{u}_2-q_{33}{u}_3-q_{36}{u}_6)\mathrm{d}t,\\ \mathrm{d}{u}_4=-{\alpha }_1{u}_4\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{u}_5=-{\alpha }_2{u}_5\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{u}_6=-{\alpha }_3{u}_6\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t),\end{array}$(25) ) with the initial value $(u_1(0),u_2(0),u_3(0),u_4(0),u_5(0),u_6(0){)}^T\in {\mathbb{R}}^6$. If Assumption (A) holds, then there is a local density function of a multivariate normal distribution $\mathrm{\Phi }(u_1,u_2,u_3,u_4,u_5,u_6)$ around the origin point $(0,0,0,0,0,0{)}^T$, which takes the form $\mathrm{\Phi }(u_1,u_2,u_3,u_4,u_5,u_6)=(2\pi {)}^{-3}|\mathrm{\Sigma }{|}^{-\frac{1}{2}}{\mathrm{e}}^{-\frac{1}{2}(u_1,u_2,u_3,u_4,u_5,u_6){\mathrm{\Sigma }}^{-1}(u_1,u_2,u_3,u_4,u_5,u_6{)}^T},$where $\mathrm{\Sigma }=({\sigma }_{ij}{)}_{6\times 6}$ denotes the covariance matrix of $(u_1,u_2,u_3,u_4,u_5,u_6{)}^T$ which is positive definite and takes the following form

1. If ${\omega }_1=0$, ${\omega }_2=0$, ${\omega }_3=0$, then $\begin{array}{rl}\mathrm{\Sigma }& ={\tilde{\rho }}_1^2(J_6{J}_5{J}_2{J}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{01}[(J_6{J}_5{J}_2{J}_1{)}^{-1}{]}^T+{\tilde{\rho }}_2^2(P_6{P}_5{P}_2{P}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{02}[(P_6{P}_5{P}_2{P}_1{)}^{-1}{]}^T\\ & +{\tilde{\rho }}_3^2(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{03}[(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{]}^T.\end{array}$

2. If ${\omega }_1=0$, ${\omega }_2\ne 0$, ${\omega }_3\ne 0$, then $\begin{array}{rl}\mathrm{\Sigma }& ={\tilde{\rho }}_1^2(J_6{J}_5{J}_2{J}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{01}[(J_6{J}_5{J}_2{J}_1{)}^{-1}{]}^T+{\rho }_2^2(P_4{P}_3{P}_2{P}_1{)}^{-1}{\mathrm{\Sigma }}_{02}[(P_4{P}_3{P}_2{P}_1{)}^{-1}{]}^T\\ & +{\rho }_3^2(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{\mathrm{\Sigma }}_{03}[(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{]}^T.\end{array}$

3. If ${\omega }_1\ne 0$, ${\omega }_2=0$, ${\omega }_3\ne 0$, then $\begin{array}{rl}\mathrm{\Sigma }& ={\rho }_1^2(J_4{J}_3{J}_2{J}_1{)}^{-1}{\mathrm{\Sigma }}_{01}[(J_4{J}_3{J}_2{J}_1{)}^{-1}{]}^T+{\tilde{\rho }}_2^2(P_6{P}_5{P}_2{P}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{02}[(P_6{P}_5{P}_2{P}_1{)}^{-1}{]}^T\\ & +{\rho }_3^2(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{\mathrm{\Sigma }}_{03}[(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{]}^T.\end{array}$

4. If ${\omega }_1\ne 0$, ${\omega }_2\ne 0$, ${\omega }_3=0$, then $\begin{array}{rl}\mathrm{\Sigma }& ={\rho }_1^2(J_4{J}_3{J}_2{J}_1{)}^{-1}{\mathrm{\Sigma }}_{01}[(J_4{J}_3{J}_2{J}_1{)}^{-1}{]}^T+{\rho }_2^2(P_4{P}_3{P}_2{P}_1{)}^{-1}{\mathrm{\Sigma }}_{02}[(P_4{P}_3{P}_2{P}_1{)}^{-1}{]}^T\\ & +{\tilde{\rho }}_3^2(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{03}[(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{]}^T.\end{array}$

5. If ${\omega }_1=0$, ${\omega }_2=0$, ${\omega }_3\ne 0$, then $\begin{array}{rl}\mathrm{\Sigma }& ={\tilde{\rho }}_1^2(J_6{J}_5{J}_2{J}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{01}[(J_6{J}_5{J}_2{J}_1{)}^{-1}{]}^T+{\tilde{\rho }}_2^2(P_6{P}_5{P}_2{P}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{02}[(P_6{P}_5{P}_2{P}_1{)}^{-1}{]}^T\\ & +{\rho }_3^2(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{\mathrm{\Sigma }}_{03}[(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{]}^T.\end{array}$

6. If ${\omega }_1=0$, ${\omega }_2\ne 0$, ${\omega }_3=0$, then $\begin{array}{rl}\mathrm{\Sigma }& ={\tilde{\rho }}_1^2(J_6{J}_5{J}_2{J}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{01}[(J_6{J}_5{J}_2{J}_1{)}^{-1}{]}^T+{\rho }_2^2(P_4{P}_3{P}_2{P}_1{)}^{-1}{\mathrm{\Sigma }}_{02}[(P_4{P}_3{P}_2{P}_1{)}^{-1}{]}^T\\ & +{\tilde{\rho }}_3^2(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{03}[(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{]}^T.\end{array}$

7. If ${\omega }_1\ne 0$, ${\omega }_2=0$, ${\omega }_3=0$, then $\begin{array}{rl}\mathrm{\Sigma }& ={\rho }_1^2(J_4{J}_3{J}_2{J}_1{)}^{-1}{\mathrm{\Sigma }}_{01}[(J_4{J}_3{J}_2{J}_1{)}^{-1}{]}^T+{\tilde{\rho }}_2^2(P_6{P}_5{P}_2{P}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{02}[(P_6{P}_5{P}_2{P}_1{)}^{-1}{]}^T\\ & +{\tilde{\rho }}_3^2(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{03}[(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{]}^T.\end{array}$

8. If ${\omega }_1\ne 0$, ${\omega }_2\ne 0$, ${\omega }_3\ne 0$, then $\begin{array}{rl}\mathrm{\Sigma }& ={\rho }_1^2(J_4{J}_3{J}_2{J}_1{)}^{-1}{\mathrm{\Sigma }}_{01}[(J_4{J}_3{J}_2{J}_1{)}^{-1}{]}^T+{\rho }_2^2(P_4{P}_3{P}_2{P}_1{)}^{-1}{\mathrm{\Sigma }}_{02}[(P_4{P}_3{P}_2{P}_1{)}^{-1}{]}^T\\ & +{\rho }_3^2(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{\mathrm{\Sigma }}_{03}[(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{]}^T,\end{array}$

here $\begin{array}{rl}{\omega }_1& =\frac{q_{21}{q}_{22}{q}_{31}-q_{21}^2{q}_{32}-q_{21}{q}_{31}{q}_{33}-q_{23}{q}_{31}^2}{q_{31}^2},\\ {\omega }_2& =\frac{q_{11}{q}_{12}{q}_{32}-q_{12}^2{q}_{31}-q_{13}{q}_{32}^2-q_{12}{q}_{32}{q}_{33}}{q_{32}^2},\\ {\omega }_3& =\frac{q_{13}{q}_{22}{q}_{23}+q_{13}^2{q}_{21}-q_{12}{q}_{23}^2-q_{11}{q}_{13}{q}_{23}}{q_{23}^2},{\rho }_1={\sigma }_1{q}_{14}{q}_{31}{\omega }_1,\\ {\tilde{\rho }}_1& ={\sigma }_1{q}_{14}{q}_{31},{\rho }_2={\sigma }_2{q}_{25}{q}_{32}{\omega }_2,\\ {\tilde{\rho }}_2& ={\sigma }_2{q}_{25}{q}_{32},{\rho }_3={\sigma }_3{q}_{23}{q}_{36}{\omega }_3,{\tilde{\rho }}_3={\sigma }_3{q}_{23}{q}_{36},\\ {\tau }_1& =q_{11}+q_{22}+q_{33},{\tau }_2=q_{11}{q}_{22}+q_{11}{q}_{33}+q_{13}{q}_{31}\\ & +q_{22}{q}_{33}+q_{23}{q}_{32}-q_{12}{q}_{21},\\ {\tau }_3& =q_{11}{q}_{22}{q}_{33}+q_{11}{q}_{23}{q}_{32}+q_{13}{q}_{22}{q}_{31}-q_{13}{q}_{21}{q}_{32}\\ & -q_{12}{q}_{21}{q}_{33}-q_{12}{q}_{23}{q}_{31},\\ p_1& =q_{11}+q_{33}+\frac{q_{21}{q}_{32}}{q_{31}},p_2=q_{11}{q}_{33}+q_{12}{q}_{21}+q_{13}{q}_{31}+\frac{q_{11}{q}_{21}{q}_{32}}{q_{31}},\\ p_3& =({\alpha }_1+\frac{q_{21}{q}_{32}}{q_{31}}-q_{22})[q_{32}(q_{11}-q_{22}+\frac{q_{21}{q}_{32}}{q_{31}})-q_{12}{q}_{31}],\\ r_1& =q_{22}+q_{33}+\frac{q_{12}{q}_{31}}{q_{32}},r_2=q_{22}(q_{33}+\frac{q_{12}{q}_{31}}{q_{32}})+q_{23}{q}_{32}-q_{12}{q}_{21},\\ r_3& =({\alpha }_2+\frac{q_{12}{q}_{31}}{q_{32}}-q_{11})[q_{31}(q_{22}-q_{11}+\frac{q_{12}{q}_{31}}{q_{32}})-q_{21}{q}_{32}],\\ l_1& =q_{22}+q_{33}+\frac{q_{13}{q}_{21}}{q_{23}},l_2=q_{13}{q}_{31}+q_{22}{q}_{33}+q_{23}{q}_{32}+\frac{q_{13}{q}_{21}{q}_{33}}{q_{23}},\\ l_3& =({\alpha }_3+\frac{q_{13}{q}_{21}}{q_{23}}-q_{11})[q_{21}(q_{11}-q_{33}-\frac{q_{13}{q}_{21}}{q_{23}})-q_{23}{q}_{31}],\\ a_{3(25)}& =q_{25}(q_{32}{\omega }_1+\frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2}),\\ a_{3(26)}& =q_{36}{\omega }_1(q_{22}+q_{33})-\frac{q_{21}{q}_{36}(q_{32}{\omega }_1+\frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2})}{q_{31}},\\ a_{5(23)}& =q_{12}{q}_{21}-q_{11}(q_{33}+\frac{q_{21}{q}_{32}}{q_{31}})-q_{13}{q}_{31},\\ a_{5(24)}& =q_{11}{q}_{32}+\frac{q_{21}{q}_{32}^2}{q_{31}}-q_{22}{q}_{32}-q_{12}{q}_{31},\\ c_{3(25)}& =q_{14}(q_{31}{\omega }_2+\frac{(q_{12}{q}_{31}-q_{11}{q}_{32}{)}^2}{q_{32}^2}),\\ c_{3(26)}& =q_{36}{\omega }_2(q_{11}+q_{33})-\frac{q_{12}{q}_{36}}{q_{32}}(q_{31}{\omega }_2+\frac{(q_{12}{q}_{31}-q_{11}{q}_{32}{)}^2}{q_{32}^2}),\\ c_{5(23)}& =q_{12}{q}_{21}-q_{22}(q_{33}+\frac{q_{12}{q}_{31}}{q_{32}})-q_{23}{q}_{32},\\ c_{5(24)}& =q_{31}(q_{22}-q_{11}+\frac{q_{12}{q}_{31}}{q_{32}})-q_{21}{q}_{32},\\ d_{3(25)}& =q_{14}(\frac{(q_{11}{q}_{23}-q_{13}{q}_{21}{)}^2}{q_{23}^2}-{\omega }_3{q}_{21}),\\ d_{3(26)}& =-q_{25}[\frac{q_{13}(q_{11}{q}_{23}-q_{13}{q}_{21}{)}^2}{q_{23}^3}+{\omega }_3(q_{11}+q_{22}-\frac{q_{13}{q}_{21}}{q_{23}})],\\ d_{5(23)}& =-(q_{13}{q}_{31}+q_{22}{q}_{33}+q_{23}{q}_{32}+\frac{q_{13}{q}_{21}{q}_{33}}{q_{23}}),\\ d_{5(24)}& =q_{21}(q_{11}-q_{33}-\frac{q_{13}{q}_{21}}{q_{23}})-q_{23}{q}_{31},\\ {\mathrm{\Sigma }}_{01}^{\ast }& =2({\alpha }_1^3+{\alpha }_1^2{\tau }_1+{\alpha }_1{\tau }_2+{\tau }_3)({\tau }_1{\tau }_2-{\tau }_3),\\ {\mathrm{\Sigma }}_{02}^{\ast }& =2({\alpha }_2^3+{\alpha }_2^2{\tau }_1+{\alpha }_2{\tau }_2+{\tau }_3)({\tau }_1{\tau }_2-{\tau }_3),\\ {\mathrm{\Sigma }}_{03}^{\ast }& =2({\alpha }_3^3+{\alpha }_3^2{\tau }_1+{\alpha }_3{\tau }_2+{\tau }_3)({\tau }_1{\tau }_2-{\tau }_3),\end{array}$and $\begin{array}{rl}{J}_1& =\left(\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),J_2=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& {\displaystyle \frac{q_{21}}{q_{31}}}& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ J_3& =\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& {\omega }_1{q}_{31}& -{\omega }_1(q_{22}+q_{33})& {\omega }_1{q}_{32}+{\displaystyle \frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2}}& 0& 0\\ 0& 0& {\omega }_1& {\displaystyle \frac{q_{21}{q}_{32}-q_{22}{q}_{31}}{q_{31}}}& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ J_4& =\left(\begin{array}{cccccc}{q}_{14}{q}_{31}{\omega }_1& -{\tau }_1& -{\tau }_2& -{\tau }_3& a_{3(25)}& a_{3(26)}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ J_5& =\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& q_{31}& -q_{33}-{\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}}& q_{32}& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ J_6& =\left(\begin{array}{cccccc}{q}_{14}{q}_{31}& -(q_{11}+q_{33}+{\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}})& a_{5(23)}& a_{5(24)}& q_{25}{q}_{32}& q_{33}{q}_{36}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ P_1& =\left(\begin{array}{cccccc}0& 0& 0& 0& 1& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),P_2=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& {\displaystyle \frac{q_{12}}{q_{32}}}& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ P_3& =\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& {\omega }_2{q}_{32}& -{\omega }_2(q_{11}+q_{33})& {\omega }_2{q}_{31}+{\displaystyle \frac{(q_{12}{q}_{31}-q_{11}{q}_{32}{)}^2}{q_{32}^2}}& 0& 0\\ 0& 0& {\omega }_2& {\displaystyle \frac{q_{12}{q}_{31}-q_{11}{q}_{32}}{q_{32}}}& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ P_4& =\left(\begin{array}{cccccc}{q}_{25}{q}_{32}{\omega }_2& -{\tau }_1& -{\tau }_2& -{\tau }_3& c_{3(25)}& c_{3(26)}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ P_5& =\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& q_{32}& -q_{33}-{\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}}& q_{31}& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ P_6& =\left(\begin{array}{cccccc}{q}_{25}{q}_{32}& -(q_{22}+q_{33}+{\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}})& c_{5(23)}& c_{5(24)}& q_{14}{q}_{31}& q_{33}{q}_{36}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ Q_1& =\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 1\\ 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\end{array}\right),Q_2=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& -{\displaystyle \frac{q_{13}}{q_{23}}}& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ Q_3& =\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& -{\omega }_3{q}_{23}& -{\omega }_3(q_{11}+q_{22})& {\displaystyle \frac{(q_{13}{q}_{21}-q_{11}{q}_{23}{)}^2}{q_{23}^2}-{\omega }_3{q}_{21}}& 0& 0\\ 0& 0& {\omega }_3& {\displaystyle \frac{q_{13}{q}_{21}-q_{11}{q}_{23}}{q_{23}}}& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ Q_4& =\left(\begin{array}{cccccc}{q}_{23}{q}_{36}{\omega }_3& -{\tau }_1& -{\tau }_2& -{\tau }_3& d_{3(25)}& d_{3(26)}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ Q_5& =\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& -q_{23}& -q_{22}-{\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}}& -q_{21}& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)\\ Q_6& =\left(\begin{array}{cccccc}{q}_{23}{q}_{36}& -(q_{22}+q_{33}+{\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}})& d_{5(23)}& d_{5(24)}& -q_{14}{q}_{21}& -q_{22}{q}_{25}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),\\ {\mathrm{\Sigma }}_{01}& =(\begin{array}{ccc}{\displaystyle \frac{{\alpha }_1^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_1{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0& -{\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}\\ 0& {\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0\\ -{\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0& {\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}\\ 0& -{\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ -{\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0& 0\\ 0& 0& 0\\ {\displaystyle \frac{{\alpha }_1{\tau }_1({\alpha }_1+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_1{\tau }_3{\mathrm{\Sigma }}_{01}^{\ast }}}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}),\\ {\tilde{\mathrm{\Sigma }}}_{01}& =(\begin{array}{ccc}{\displaystyle \frac{{\alpha }_1{p}_1+p_2}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0& -{\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}\\ 0& {\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0\\ -{\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0& {\displaystyle \frac{{\alpha }_1+p_1}{2{\alpha }_1{p}_1{p}_2({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}),\\ {\mathrm{\Sigma }}_{02}& =(\begin{array}{ccc}{\displaystyle \frac{{\alpha }_2^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_2{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0& -{\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}\\ 0& {\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0\\ -{\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0& {\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}\\ 0& -{\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ -{\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0& 0\\ 0& 0& 0\\ {\displaystyle \frac{{\alpha }_2{\tau }_1({\alpha }_2+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_2{\tau }_3{\mathrm{\Sigma }}_{02}^{\ast }}}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}),\\ {\tilde{\mathrm{\Sigma }}}_{02}& =(\begin{array}{cc}{\displaystyle \frac{{\alpha }_2{r}_1+r_2}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0\\ 0& {\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}\\ -{\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\\ & \begin{array}{cccc}-{\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0& 0& 0\\ 0& 0& 0& 0\\ {\displaystyle \frac{{\alpha }_2+r_1}{2{\alpha }_2{r}_1{r}_2({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}),\\ {\mathrm{\Sigma }}_{03}& =(\begin{array}{ccc}{\displaystyle \frac{{\alpha }_3^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_3{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0& -{\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}\\ 0& {\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0\\ -{\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0& {\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}\\ 0& -{\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ -{\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0& 0\\ 0& 0& 0\\ {\displaystyle \frac{{\alpha }_3{\tau }_1({\alpha }_3+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_3{\tau }_3{\mathrm{\Sigma }}_{03}^{\ast }}}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array})\\ ,\\ {\tilde{\mathrm{\Sigma }}}_{03}& =(\begin{array}{cc}{\displaystyle \frac{{\alpha }_3{l}_1+l_2}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0\\ 0& {\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}\\ -{\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\\ & \begin{array}{cccc}-{\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0& 0& 0\\ 0& 0& 0& 0\\ {\displaystyle \frac{{\alpha }_3+l_1}{2{\alpha }_3{l}_1{l}_2({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}).\end{array}$

Proof.

Let $X=(u_1,u_2,u_3,u_4,u_5,u_6{)}^T$, $B(t)=(0,0,0,B_1(t),B_2(t),B_3(t){)}^T$, $\begin{array}{rl}A& =\left(\begin{array}{cccccc}-q_{11}& -q_{12}& -q_{13}& q_{14}& 0& 0\\ -q_{21}& -q_{22}& -q_{23}& 0& q_{25}& 0\\ q_{31}& q_{32}& -q_{33}& 0& 0& -q_{36}\\ 0& 0& 0& -{\alpha }_1& 0& 0\\ 0& 0& 0& 0& -{\alpha }_2& 0\\ 0& 0& 0& 0& 0& -{\alpha }_3\end{array}\right),\\ H& =\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\sigma }_1& 0& 0\\ 0& 0& 0& 0& {\sigma }_2& 0\\ 0& 0& 0& 0& 0& {\sigma }_3\end{array}\right).\end{array}$With these notations, system (25(25) $\{\begin{array}{l}\mathrm{d}{u}_1=(-q_{11}{u}_1-q_{12}{u}_2-q_{13}{u}_3+q_{14}{u}_4)\mathrm{d}t,\\ \mathrm{d}{u}_2=(-q_{21}{u}_1-q_{22}{u}_2-q_{23}{u}_3+q_{25}{u}_5)\mathrm{d}t,\\ \mathrm{d}{u}_3=(q_{31}{u}_1+q_{32}{u}_2-q_{33}{u}_3-q_{36}{u}_6)\mathrm{d}t,\\ \mathrm{d}{u}_4=-{\alpha }_1{u}_4\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{u}_5=-{\alpha }_2{u}_5\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{u}_6=-{\alpha }_3{u}_6\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t),\end{array}$(25) ) can be written in the following equivalent form: (26) $\mathrm{d}X(t)=AX(t)\mathrm{d}t+H\mathrm{d}B(t).$(26) On the basis of Gardiner [6], system (26(26) $\mathrm{d}X(t)=AX(t)\mathrm{d}t+H\mathrm{d}B(t).$(26) ) has a unique density function $\mathrm{\Phi }(u_1,u_2,u_3,u_4,u_5,u_6,t)$, which can be determined by the following six-dimensional Fokker–Planck equation (27) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}\mathrm{\Phi }(X(t),t)+\frac{\mathrm{\partial }}{\mathrm{\partial }X(t)}[AX(t)\mathrm{\Phi }(X(t),t)]-\frac{{\sigma }_1^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_4^2}\mathrm{\Phi }(X(t),t)-\frac{{\sigma }_2^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_5^2}\mathrm{\Phi }(X(t),t)\\ & -\frac{{\sigma }_3^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_6^2}\mathrm{\Phi }(X(t),t)=0.\end{array}$(27) Next, we will give the accurate expression of the density function by solving Equation (27(27) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}\mathrm{\Phi }(X(t),t)+\frac{\mathrm{\partial }}{\mathrm{\partial }X(t)}[AX(t)\mathrm{\Phi }(X(t),t)]-\frac{{\sigma }_1^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_4^2}\mathrm{\Phi }(X(t),t)-\frac{{\sigma }_2^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_5^2}\mathrm{\Phi }(X(t),t)\\ & -\frac{{\sigma }_3^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_6^2}\mathrm{\Phi }(X(t),t)=0.\end{array}$(27) ). Note that $\mathrm{\partial }\mathrm{\Phi }(X(t),t)/\mathrm{\partial }t=0$ under a stationary case, then (27(27) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}\mathrm{\Phi }(X(t),t)+\frac{\mathrm{\partial }}{\mathrm{\partial }X(t)}[AX(t)\mathrm{\Phi }(X(t),t)]-\frac{{\sigma }_1^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_4^2}\mathrm{\Phi }(X(t),t)-\frac{{\sigma }_2^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_5^2}\mathrm{\Phi }(X(t),t)\\ & -\frac{{\sigma }_3^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_6^2}\mathrm{\Phi }(X(t),t)=0.\end{array}$(27) ) becomes $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }{u}_1}[(-q_{11}{u}_1-q_{12}{u}_2-q_{13}{u}_3+q_{14}{u}_4)\mathrm{\Phi }]+\frac{\mathrm{\partial }}{\mathrm{\partial }{u}_2}[(-q_{21}{u}_1-q_{22}{u}_2-q_{23}{u}_3+q_{25}{u}_5)\mathrm{\Phi }]\\ & +\frac{\mathrm{\partial }}{\mathrm{\partial }{u}_3}[(q_{31}{u}_1+q_{32}{u}_2-q_{33}{u}_3-q_{36}{u}_6)\mathrm{\Phi }]\\ & +\frac{\mathrm{\partial }}{\mathrm{\partial }{u}_4}(-{\alpha }_1{u}_4\mathrm{\Phi })+\frac{\mathrm{\partial }}{\mathrm{\partial }{u}_5}(-{\alpha }_2{u}_5\mathrm{\Phi })+\frac{\mathrm{\partial }}{\mathrm{\partial }{u}_6}(-{\alpha }_3{u}_6\mathrm{\Phi })-\frac{{\sigma }_1^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_4^2}\mathrm{\Phi }\\ & -\frac{{\sigma }_2^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_5^2}\mathrm{\Phi }-\frac{{\sigma }_3^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{u}_6^2}\mathrm{\Phi }=0.\end{array}$Additionally, because the diffusion matrix H of system (26(26) $\mathrm{d}X(t)=AX(t)\mathrm{d}t+H\mathrm{d}B(t).$(26) ) is a constant matrix, then by the work of Gardiner [6], we get (28) $H^2+A\mathrm{\Sigma }+\mathrm{\Sigma }{A}^T=0.$(28) Because of the mutual independence of Brownian motions $B_1(t)$, $B_2(t)$ and $B_3(t)$, in view of the finite independent superposition principle, system (28(28) $H^2+A\mathrm{\Sigma }+\mathrm{\Sigma }{A}^T=0.$(28) ) can be equivalent to the sum of solutions of the following three sub-equations $H_i^2+A{\mathrm{\Sigma }}_i+{\mathrm{\Sigma }}_i{A}^T=0,i=1,2,3,$where $H_1^2=\mathrm{diag}(0,0,0,{\sigma }_1^2,0,0)$, $H_2^2=\mathrm{diag}(0,0,0,0,{\sigma }_2^2,0)$, $H_3^2=\mathrm{diag}(0,0,0,0,0,{\sigma }_3^2)$, ${\mathrm{\Sigma }}_i$ $(i=1,2,3)$ are their corresponding solutions. One can easily get that $\mathrm{\Sigma }={\mathrm{\Sigma }}_1+{\mathrm{\Sigma }}_2+{\mathrm{\Sigma }}_3$ and $H^2=H_1^2+H_2^2+H_3^2$.

Let $A^{(3)}:=\left(\begin{array}{ccc}-q_{11}& -q_{12}& -q_{13}\\ -q_{21}& -q_{22}& -q_{23}\\ q_{31}& q_{32}& -q_{33}\end{array}\right).$In order to prove the matrix $A^{(3)}$ is a Hurwitz matrix, in view of Definition 6.1, we only need to show that all the eigenvalues of $A^{(3)}$ have negative real parts. To this end, define the characteristic equation of $A^{(3)}$ by ${\phi }_{A^{(3)}}(\lambda )={\lambda }^3+{\tau }_1{\lambda }^2+{\tau }_2\lambda +{\tau }_3=0,$where $\begin{array}{rl}{\tau }_1& =q_{11}+q_{22}+q_{33}>0,{\tau }_2=q_{11}{q}_{22}+q_{11}{q}_{33}+q_{13}{q}_{31}+q_{22}{q}_{33}\\ & +q_{23}{q}_{32}-q_{12}{q}_{21}>0,\\ {\tau }_3& =q_{11}{q}_{22}{q}_{33}+q_{11}{q}_{23}{q}_{32}+q_{13}{q}_{22}{q}_{31}-q_{13}{q}_{21}{q}_{32}-q_{12}{q}_{21}{q}_{33}-q_{12}{q}_{23}{q}_{31}>0.\end{array}$By direct computation, one has $\begin{array}{rl}{\tau }_1{\tau }_2-{\tau }_3& =(q_{11}+q_{22}+q_{33})(q_{11}{q}_{22}+q_{11}{q}_{33}+q_{13}{q}_{31}+q_{22}{q}_{33}+q_{23}{q}_{32}-q_{12}{q}_{21})\\ & -(q_{11}{q}_{22}{q}_{33}+q_{11}{q}_{23}{q}_{32}\\ & +q_{13}{q}_{22}{q}_{31}-q_{13}{q}_{21}{q}_{32}-q_{12}{q}_{21}{q}_{33}-q_{12}{q}_{23}{q}_{31}{)}^2\\ & =x_1^{\ast }{x}_2^{\ast }{x}_3^{\ast }(b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}+2b_{11}{b}_{22}{b}_{33})+(x_1^{\ast }{)}^2{x}_2^{\ast }(b_{11}^2{b}_{22}-b_{11}{b}_{12}{b}_{21})\\ & +x_1^{\ast }(x_2^{\ast }{)}^2(b_{11}{b}_{22}^2-b_{12}{b}_{21}{b}_{22})\\ & +(x_1^{\ast }{)}^2{x}_3^{\ast }(b_{11}^2{b}_{33}+b_{11}{b}_{13}{b}_{31})+x_1^{\ast }(x_3^{\ast }{)}^2(b_{13}{b}_{31}{b}_{33}+b_{11}{b}_{33}^2)\\ & +(x_2^{\ast }{)}^2{x}_3^{\ast }(b_{22}^2{b}_{33}+b_{22}{b}_{23}{b}_{32})\\ & +x_2^{\ast }(x_3^{\ast }{)}^2(b_{22}{b}_{33}^2+b_{23}{b}_{32}{b}_{33})\\ & >0.\end{array}$By means of Lemma 6.1, we obtain that the matrix $A^{(3)}$ has all negative real-part eigenvalues.

Next, we will accomplish the proof of Theorem 6.1 in three steps.

Step 1. Consider the algebraic equation (29) $H_1^2+A{\mathrm{\Sigma }}_1+{\mathrm{\Sigma }}_1{A}^T=0,$(29) where $H_1^2=\mathrm{diag}(0,0,0,{\sigma }_1^2,0,0)$.

Define $A_1=J_{1}A{J}_1^{-1}$, where the ordering matrix $J_1$ is given by $J_1=\left(\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right).$Direct calculation leads to that $A_1=\left(\begin{array}{cccccc}-{\alpha }_1& 0& 0& 0& 0& 0\\ q_{14}& -q_{11}& -q_{13}& -q_{12}& 0& 0\\ 0& q_{31}& -q_{33}& q_{32}& 0& -q_{36}\\ 0& -q_{21}& -q_{23}& -q_{22}& q_{25}& 0\\ 0& 0& 0& 0& -{\alpha }_2& 0\\ 0& 0& 0& 0& 0& -{\alpha }_3\end{array}\right).$Define $A_2=J_2{A}_1{J}_2^{-1}$, where the elimination matrix $J_2$ takes the form $J_2=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& {\displaystyle \frac{q_{21}}{q_{31}}}& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$then we obtain $A_2=\left(\begin{array}{cccccc}-{\alpha }_1& 0& 0& 0& 0& 0\\ q_{14}& -q_{11}& {\displaystyle \frac{q_{12}{q}_{21}}{q_{31}}-q_{13}}& -q_{12}& 0& 0\\ 0& q_{31}& -q_{33}-{\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}}& q_{32}& 0& -q_{36}\\ 0& 0& {\omega }_1& {\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}-q_{22}}& q_{25}& -{\displaystyle \frac{q_{21}{q}_{36}}{q_{31}}}\\ 0& 0& 0& 0& -{\alpha }_2& 0\\ 0& 0& 0& 0& 0& -{\alpha }_3\end{array}\right),$where ${\omega }_1=(q_{21}{q}_{22}{q}_{31}-q_{21}^2{q}_{32}-q_{21}{q}_{31}{q}_{33}-q_{23}{q}_{31}^2)/q_{31}^2$. In view of the value of ${\omega }_1$, we will divide the following discussion into two parts.

Case 1. If ${\omega }_1\ne 0$, consider the standard transformation matrix $J_3=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& {\omega }_1{q}_{31}& -{\omega }_1(q_{22}+q_{33})& {\omega }_1{q}_{32}+{\displaystyle \frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2}}& 0& 0\\ 0& 0& {\omega }_1& {\displaystyle \frac{q_{21}{q}_{32}-q_{22}{q}_{31}}{q_{31}}}& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{A}_3& =J_3{A}_2{J}_3^{-1}\\ & =(\begin{array}{ccc}-{\alpha }_1& 0& 0\\ q_{14}{q}_{31}{\omega }_1& -(q_{11}+q_{22}+q_{33})& a_{3(23)}\\ 0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ a_{3(24)}& a_{3(25)}& a_{3(26)}\\ 0& {\displaystyle \frac{q_{25}(q_{21}{q}_{32}-q_{22}{q}_{31})}{q_{31}}}& a_{3(36)}\\ 0& q_{25}& -{\displaystyle \frac{q_{21}{q}_{36}}{q_{31}}}\\ 0& -{\alpha }_2& 0\\ 0& 0& -{\alpha }_3\end{array}),\end{array}$where $\begin{array}{rl}{a}_{3(23)}& =q_{12}{q}_{21}-q_{11}{q}_{22}-q_{11}{q}_{33}-q_{13}{q}_{31}-q_{22}{q}_{33}-q_{23}{q}_{32},\\ a_{3(24)}& =q_{13}{q}_{21}{q}_{32}+q_{12}{q}_{21}{q}_{33}+q_{12}{q}_{23}{q}_{31}-q_{11}{q}_{22}{q}_{33}-q_{11}{q}_{23}{q}_{32}-q_{13}{q}_{22}{q}_{31},\\ a_{3(25)}& =q_{25}(q_{32}{\omega }_1+\frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2}),\\ a_{3(26)}& =q_{36}{\omega }_1(q_{22}+q_{33})-\frac{q_{21}{q}_{36}(q_{32}{\omega }_1+\frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2})}{q_{31}},\\ a_{3(36)}& =-q_{36}{\omega }_1-\frac{q_{21}{q}_{36}(q_{21}{q}_{32}-q_{22}{q}_{31})}{q_{31}^2}.\end{array}$Moreover, define $A_4=J_4{A}_3{J}_4^{-1}$, where the standard transformation matrix $J_4$ takes the form $J_4=\left(\begin{array}{cccccc}{q}_{14}{q}_{31}{\omega }_1& -{\tau }_1& -{\tau }_2& -{\tau }_3& a_{3(25)}& a_{3(26)}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right).$Direct calculation gives that $\begin{array}{rl}{A}_4& =(\begin{array}{ccc}-({\tau }_1+{\alpha }_1)& -({\tau }_1{\alpha }_1+{\tau }_2)& -({\tau }_2{\alpha }_1+{\tau }_3)\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}-{\tau }_3{\alpha }_1& a_{4(15)}& a_{4(16)}\\ 0& 0& 0\\ 0& {\displaystyle \frac{q_{25}(q_{21}{q}_{32}-q_{22}{q}_{31})}{q_{31}}}& a_{4(36)}\\ 0& q_{25}& -{\displaystyle \frac{q_{21}{q}_{36}}{q_{31}}}\\ 0& -{\alpha }_2& 0\\ 0& 0& -{\alpha }_3\end{array}),\end{array}$where $\begin{array}{rl}{a}_{4(15)}& =q_{25}\{({\alpha }_1-{\alpha }_2)(q_{32}{\omega }_1+\frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2})-{\tau }_3-\frac{{\tau }_2(q_{21}{q}_{32}-q_{22}{q}_{31})}{q_{31}}\},\\ a_{4(16)}& =\frac{q_{21}{q}_{36}{\tau }_3}{q_{31}}+({\alpha }_1-{\alpha }_3-{\tau }_2)\\ & \times [q_{36}{\omega }_1(q_{22}+q_{33})-\frac{q_{21}{q}_{36}(q_{32}{\omega }_1+\frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2})}{q_{31}}],\\ a_{4(36)}& =q_{36}{\omega }_1(q_{22}+q_{33})-\frac{q_{21}{q}_{36}(q_{32}{\omega }_1+\frac{(q_{21}{q}_{32}-q_{22}{q}_{31}{)}^2}{q_{31}^2})}{q_{31}}.\end{array}$In addition, an equivalent algebraic equation of (6.6) can be described by $(J_4{J}_3{J}_2{J}_1)H_1^2(J_4{J}_3{J}_2{J}_1{)}^T+A_4(J_4{J}_3{J}_2{J}_1){\mathrm{\Sigma }}_1(J_4{J}_3{J}_2{J}_1{)}^T+(J_4{J}_3{J}_2{J}_1){\mathrm{\Sigma }}_1(J_4{J}_3{J}_2{J}_1{)}^T{A}_4^T=0.$Let ${\mathrm{\Sigma }}_{01}=\frac{1}{{\rho }_1^2}(J_4{J}_3{J}_2{J}_1){\mathrm{\Sigma }}_1(J_4{J}_3{J}_2{J}_1{)}^T,G_0=\mathrm{diag}(1,0,0,0,0,0),$where ${\rho }_1={\sigma }_1{q}_{14}{q}_{31}{\omega }_1$, then we have (30) $G_0^2+A_4{\mathrm{\Sigma }}_{01}+{\mathrm{\Sigma }}_{01}{A}_4^T=0.$(30) By solving Equation (30(30) $G_0^2+A_4{\mathrm{\Sigma }}_{01}+{\mathrm{\Sigma }}_{01}{A}_4^T=0.$(30) ), we obtain $\begin{array}{rl}{\mathrm{\Sigma }}_{01}& =(\begin{array}{ccc}{\displaystyle \frac{{\alpha }_1^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_1{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0& -{\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}\\ 0& {\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0\\ -{\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0& {\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}\\ 0& -{\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ -{\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0& 0\\ 0& 0& 0\\ {\displaystyle \frac{{\alpha }_1{\tau }_1({\alpha }_1+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_1{\tau }_3{\mathrm{\Sigma }}_{01}^{\ast }}}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}),\end{array}$where ${\mathrm{\Sigma }}_{01}^{\ast }=2({\alpha }_1^3+{\alpha }_1^2{\tau }_1+{\alpha }_1{\tau }_2+{\tau }_3)({\tau }_1{\tau }_2-{\tau }_3)$.

Noticeably, the matrix ${\mathrm{\Sigma }}_{01}$ is positive semi-definite and its submatrix $\begin{array}{rl}{\mathrm{\Sigma }}_{01}^{(4)}& =(\begin{array}{cc}{\displaystyle \frac{{\alpha }_1^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_1{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0\\ 0& {\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}\\ -{\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0\\ 0& -{\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}\end{array}\\ & \begin{array}{cc}-{\displaystyle \frac{{\alpha }_1{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0\\ 0& -{\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}\\ {\displaystyle \frac{{\alpha }_1+{\tau }_1}{{\mathrm{\Sigma }}_{01}^{\ast }}}& 0\\ 0& {\displaystyle \frac{{\alpha }_1{\tau }_1({\alpha }_1+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_1{\tau }_3{\mathrm{\Sigma }}_{01}^{\ast }}}\end{array})\end{array}$is positive definite. Thus, the matrix ${\mathrm{\Sigma }}_1={\rho }_1^2(J_4{J}_3{J}_2{J}_1{)}^{-1}{\mathrm{\Sigma }}_{01}[(J_4{J}_3{J}_2{J}_1{)}^{-1}{]}^T$ is also positive semi-definite and there is a positive constant ${\eta }_1$ such that ${\mathrm{\Sigma }}_1⪰{\eta }_1\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right).$Case 2. If ${\omega }_1=0$, consider the standardized transformation matrix $J_5=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& q_{31}& -q_{33}-{\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}}& q_{32}& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{A}_5& =J_5{A}_2{J}_5^{-1}=(\begin{array}{ccc}-{\alpha }_1& 0& 0\\ q_{14}{q}_{31}& -(q_{11}+q_{33}+{\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}})& a_{5(23)}\\ 0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ a_{5(24)}& q_{25}{q}_{32}& q_{33}{q}_{36}\\ 0& 0& -q_{36}\\ {\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}-q_{22}}& q_{25}& -{\displaystyle \frac{q_{21}{q}_{36}}{q_{31}}}\\ 0& -{\alpha }_2& 0\\ 0& 0& -{\alpha }_3\end{array}),\end{array}$where $\begin{array}{rl}{a}_{5(23)}& =q_{12}{q}_{21}-q_{11}(q_{33}+\frac{q_{21}{q}_{32}}{q_{31}})-q_{13}{q}_{31},\\ a_{5(24)}& =q_{11}{q}_{32}+\frac{q_{21}{q}_{32}^2}{q_{31}}-q_{22}{q}_{32}-q_{12}{q}_{31}.\end{array}$Then continue to do a matrix similarity transformation to $A_5$, choose the matrix $J_6=\left(\begin{array}{cccccc}{q}_{14}{q}_{31}& -(q_{11}+q_{33}+{\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}})& a_{5(23)}& a_{5(24)}& q_{25}{q}_{32}& q_{33}{q}_{36}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{A}_6& =J_6{A}_5{J}_6^{-1}=(\begin{array}{ccc}-(p_1+{\alpha }_1)& -(p_1{\alpha }_1+p_2)& -p_2{\alpha }_1\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}{p}_3& a_{6(15)}& a_{6(16)}\\ 0& 0& 0\\ 0& 0& -q_{36}\\ {\displaystyle \frac{q_{21}{q}_{32}}{q_{31}}-q_{22}}& q_{25}& -{\displaystyle \frac{q_{21}{q}_{36}}{q_{31}}}\\ 0& -{\alpha }_2& 0\\ 0& 0& -{\alpha }_3\end{array}),\end{array}$where $\begin{array}{rl}{p}_1& =q_{11}+q_{33}+\frac{q_{21}{q}_{32}}{q_{31}},p_2=-a_{5(23)},p_3=({\alpha }_1+\frac{q_{21}{q}_{32}}{q_{31}}-q_{22})a_{5(24)},\\ a_{6(15)}& =q_{25}(q_{11}{q}_{32}-q_{12}{q}_{31}-q_{22}{q}_{32}+\frac{q_{21}{q}_{32}^2}{q_{31}})+q_{25}{q}_{32}({\alpha }_1-{\alpha }_2),\\ a_{6(16)}& =q_{13}{q}_{31}{q}_{36}+q_{33}{q}_{36}({\alpha }_1-{\alpha }_3+q_{11})+\frac{q_{21}{q}_{32}{q}_{36}}{q_{31}}(q_{22}-\frac{q_{21}{q}_{32}}{q_{31}}).\end{array}$Furthermore, an equivalent algebraic equation of (6.6) can be described by $(J_6{J}_5{J}_2{J}_1)H_1^2(J_6{J}_5{J}_2{J}_1{)}^T+A_6(J_6{J}_5{J}_2{J}_1){\mathrm{\Sigma }}_1(J_6{J}_5{J}_2{J}_1{)}^T+(J_6{J}_5{J}_2{J}_1){\mathrm{\Sigma }}_1(J_6{J}_5{J}_2{J}_1{)}^T{A}_6^T=0.$Denote by ${\tilde{\mathrm{\Sigma }}}_{01}=\frac{1}{{\tilde{\rho }}_1^2}(J_6{J}_5{J}_2{J}_1){\mathrm{\Sigma }}_1(J_6{J}_5{J}_2{J}_1{)}^T,{\tilde{\rho }}_1={\sigma }_1{q}_{14}{q}_{31},$then (31) $G_0^2+A_6{\tilde{\mathrm{\Sigma }}}_{01}+{\tilde{\mathrm{\Sigma }}}_{01}{A}_6^T=0.$(31) By solving Equation (31(31) $G_0^2+A_6{\tilde{\mathrm{\Sigma }}}_{01}+{\tilde{\mathrm{\Sigma }}}_{01}{A}_6^T=0.$(31) ), we obtain $\begin{array}{rl}{\tilde{\mathrm{\Sigma }}}_{01}& =(\begin{array}{ccc}{\displaystyle \frac{{\alpha }_1{p}_1+p_2}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0& -{\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}\\ 0& {\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0\\ -{\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0& {\displaystyle \frac{{\alpha }_1+p_1}{2{\alpha }_1{p}_1{p}_2({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}).\end{array}$It is obvious that the matrix ${\tilde{\mathrm{\Sigma }}}_{01}$ is positive semi-definite and its submatrix $\begin{array}{rl}{\tilde{\mathrm{\Sigma }}}_{01}^{(4)}& =(\begin{array}{cc}{\displaystyle \frac{{\alpha }_1{p}_1+p_2}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0\\ 0& {\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}\\ -{\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0\\ 0& 0\end{array}\\ & \begin{array}{cc}-{\displaystyle \frac{1}{2p_1({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0\\ 0& 0\\ {\displaystyle \frac{{\alpha }_1+p_1}{2{\alpha }_1{p}_1{p}_2({\alpha }_1{p}_1+p_2+{\alpha }_1^2)}}& 0\\ 0& 0\end{array})\end{array}$is also positive semi-definite. So the matrix ${\mathrm{\Sigma }}_1={\tilde{\rho }}_1^2(J_6{J}_5{J}_2{J}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{01}[(J_6{J}_5{J}_2{J}_1{)}^{-1}{]}^T$ is positive semi-definite and there is a positive constant ${\tilde{\eta }}_1$ such that ${\mathrm{\Sigma }}_1⪰{\tilde{\eta }}_1\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right).$Step 2. Consider the algebraic equation (32) $H_2^2+A{\mathrm{\Sigma }}_2+{\mathrm{\Sigma }}_2{A}^T=0,$(32) where $H_2^2=\mathrm{diag}(0,0,0,0,{\sigma }_2^2,0)$.

Let $C_1=P_{1}A{P}_1^{-1}$, where the ordering matrix $P_1$ takes the form $P_1=\left(\begin{array}{cccccc}0& 0& 0& 0& 1& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right).$Then $C_1=\left(\begin{array}{cccccc}-{\alpha }_2& 0& 0& 0& 0& 0\\ q_{25}& -q_{22}& -q_{23}& -q_{21}& 0& 0\\ 0& q_{32}& -q_{33}& q_{31}& 0& -q_{36}\\ 0& -q_{12}& -q_{13}& -q_{11}& q_{14}& 0\\ 0& 0& 0& 0& -{\alpha }_1& 0\\ 0& 0& 0& 0& 0& -{\alpha }_3\end{array}\right).$Let $C_2=P_2{C}_1{P}_2^{-1}$, where the elimination matrix $P_2$ takes the form $P_2=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& {\displaystyle \frac{q_{12}}{q_{32}}}& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$then we get $C_2=\left(\begin{array}{cccccc}-{\alpha }_2& 0& 0& 0& 0& 0\\ q_{25}& -q_{22}& {\displaystyle \frac{q_{12}{q}_{21}}{q_{32}}-q_{23}}& -q_{21}& 0& 0\\ 0& q_{32}& -q_{33}-{\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}}& q_{31}& 0& -q_{36}\\ 0& 0& {\omega }_2& {\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}-q_{11}}& q_{14}& -{\displaystyle \frac{q_{12}{q}_{36}}{q_{32}}}\\ 0& 0& 0& 0& -{\alpha }_1& 0\\ 0& 0& 0& 0& 0& -{\alpha }_3\end{array}\right),$where ${\omega }_2=(q_{11}{q}_{12}{q}_{32}-q_{12}^2{q}_{31}-q_{13}{q}_{32}^2-q_{12}{q}_{32}{q}_{33})/q_{32}^2$. By means of the value of ${\omega }_2$, we will divide the following discussion into two parts.

Case 1. If ${\omega }_2\ne 0$, consider the standardized transformation matrix $P_3=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& {\omega }_2{q}_{32}& -{\omega }_2(q_{11}+q_{33})& {\omega }_2{q}_{31}+{\displaystyle \frac{(q_{12}{q}_{31}-q_{11}{q}_{32}{)}^2}{q_{32}^2}}& 0& 0\\ 0& 0& {\omega }_2& {\displaystyle \frac{q_{12}{q}_{31}-q_{11}{q}_{32}}{q_{32}}}& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{C}_3& =P_3{C}_2{P}_3^{-1}\\ & =(\begin{array}{cccc}-{\alpha }_2& 0& 0& 0\\ {\omega }_2{q}_{25}{q}_{32}& -(q_{11}+q_{22}+q_{33})& c_{3(23)}& c_{3(24)}\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\\ & \begin{array}{cc}0& 0\\ c_{3(25)}& c_{3(26)}\\ {\displaystyle \frac{q_{14}(q_{12}{q}_{31}-q_{11}{q}_{32})}{q_{32}}}& c_{3(36)}\\ q_{14}& -{\displaystyle \frac{q_{12}{q}_{36}}{q_{32}}}\\ -{\alpha }_1& 0\\ 0& -{\alpha }_3\end{array}),\end{array}$where $\begin{array}{rl}{c}_{3(23)}& =q_{12}{q}_{21}-q_{11}{q}_{22}-q_{11}{q}_{33}-q_{22}{q}_{33}-q_{23}{q}_{32}+q_{31}{\omega }_2+\frac{q_{12}{q}_{31}}{q_{32}}(q_{33}-q_{11})\\ & +\frac{q_{12}^2{q}_{31}^2}{q_{32}^2},\\ c_{3(24)}& ={\omega }_2(q_{22}{q}_{31}-q_{21}{q}_{32})+\frac{q_{12}{q}_{31}-q_{11}{q}_{32}}{q_{32}^2}\\ & \times (q_{12}{q}_{22}{q}_{31}+q_{22}{q}_{32}{q}_{33}+q_{23}{q}_{32}^2-q_{12}{q}_{21}{q}_{32}),\\ c_{3(25)}& =q_{14}(q_{31}{\omega }_2+\frac{(q_{12}{q}_{31}-q_{11}{q}_{32}{)}^2}{q_{32}^2}),\\ c_{3(26)}& =q_{36}{\omega }_2(q_{11}+q_{33})-\frac{q_{12}{q}_{36}}{q_{32}}(q_{31}{\omega }_2+\frac{(q_{12}{q}_{31}-q_{11}{q}_{32}{)}^2}{q_{32}^2}),\\ c_{3(36)}& =-q_{36}{\omega }_2-\frac{q_{12}{q}_{36}(q_{12}{q}_{31}-q_{11}{q}_{32})}{q_{32}^2}.\end{array}$Then continue to do a matrix similarity transformation to $C_3$, choose the matrix $P_4=\left(\begin{array}{cccccc}{q}_{25}{q}_{32}{\omega }_2& -{\tau }_1& -{\tau }_2& -{\tau }_3& c_{3(25)}& c_{3(26)}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{C}_4& =P_4{C}_3{P}_4^{-1}\\ & =(\begin{array}{ccc}-({\tau }_1+{\alpha }_2)& -({\tau }_1{\alpha }_2+{\tau }_2)& -({\tau }_2{\alpha }_2+{\tau }_3)\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}-{\tau }_3{\alpha }_2& c_{4(15)}& c_{4(16)}\\ 0& 0& 0\\ 0& {\displaystyle \frac{q_{14}(q_{12}{q}_{31}-q_{11}{q}_{32})}{q_{32}}}& c_{4(36)}\\ 0& q_{14}& -{\displaystyle \frac{q_{12}{q}_{36}}{q_{32}}}\\ 0& -{\alpha }_1& 0\\ 0& 0& -{\alpha }_3\end{array}),\end{array}$where $\begin{array}{rl}{c}_{4(15)}& =q_{14}\{({\alpha }_2-{\alpha }_1)[q_{11}^2-\frac{q_{12}{q}_{31}(q_{11}-q_{33})}{q_{32}}]-{\tau }_3+\frac{{\tau }_2(q_{11}{q}_{32}-q_{12}{q}_{31})}{q_{32}}\},\\ c_{4(16)}& =q_{36}\{({\alpha }_3-{\alpha }_2)[\frac{q_{12}}{q_{32}}(q_{33}^2-q_{13}{q}_{31})+q_{11}{q}_{13}+q_{13}{q}_{33}]\\ & -{\tau }_2(q_{13}+\frac{q_{12}{q}_{33}}{q_{32}})+\frac{q_{12}{\tau }_3}{q_{32}}\},\\ c_{4(36)}& =\frac{q_{36}(q_{12}{q}_{33}+q_{13}{q}_{32})}{q_{32}}.\end{array}$Moreover, an equivalent algebraic equation of (32(32) $H_2^2+A{\mathrm{\Sigma }}_2+{\mathrm{\Sigma }}_2{A}^T=0,$(32) ) can be described by $\begin{array}{rl}& (P_4{P}_3{P}_2{P}_1)H_2^2(P_4{P}_3{P}_2{P}_1{)}^T+C_4(P_4{P}_3{P}_2{P}_1){\mathrm{\Sigma }}_2(P_4{P}_3{P}_2{P}_1{)}^T\\ & +(P_4{P}_3{P}_2{P}_1){\mathrm{\Sigma }}_2(P_4{P}_3{P}_2{P}_1{)}^T{C}_4^T=0.\end{array}$Let ${\mathrm{\Sigma }}_{02}=\frac{1}{{\rho }_2^2}(P_4{P}_3{P}_2{P}_1){\mathrm{\Sigma }}_2(P_4{P}_3{P}_2{P}_1{)}^T,$where ${\rho }_2={\sigma }_2{q}_{25}{q}_{32}{\omega }_2$, then we have (33) $G_0^2+C_4{\mathrm{\Sigma }}_{02}+{\mathrm{\Sigma }}_{02}{C}_4^T=0.$(33) By solving Equation (33(33) $G_0^2+C_4{\mathrm{\Sigma }}_{02}+{\mathrm{\Sigma }}_{02}{C}_4^T=0.$(33) ), we have $\begin{array}{rl}{\mathrm{\Sigma }}_{02}& =(\begin{array}{ccc}{\displaystyle \frac{{\alpha }_2^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_2{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0& -{\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}\\ 0& {\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0\\ -{\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0& {\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}\\ 0& -{\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ -{\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0& 0\\ 0& 0& 0\\ {\displaystyle \frac{{\alpha }_2{\tau }_1({\alpha }_2+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_2{\tau }_3{\mathrm{\Sigma }}_{02}^{\ast }}}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}),\end{array}$where ${\mathrm{\Sigma }}_{02}^{\ast }=2({\alpha }_2^3+{\alpha }_2^2{\tau }_1+{\alpha }_2{\tau }_2+{\tau }_3)({\tau }_1{\tau }_2-{\tau }_3)$.

It is easy to see that the matrix ${\mathrm{\Sigma }}_{02}$ is positive semi-definite and its submatrix $\begin{array}{rl}{\mathrm{\Sigma }}_{02}^{(4)}& =(\begin{array}{cc}{\displaystyle \frac{{\alpha }_2^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_2{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0\\ 0& {\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}\\ -{\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0\\ 0& -{\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}\end{array}\\ & \begin{array}{cc}-{\displaystyle \frac{{\alpha }_2{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0\\ 0& -{\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}\\ {\displaystyle \frac{{\alpha }_2+{\tau }_1}{{\mathrm{\Sigma }}_{02}^{\ast }}}& 0\\ 0& {\displaystyle \frac{{\alpha }_2{\tau }_1({\alpha }_2+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_2{\tau }_3{\mathrm{\Sigma }}_{02}^{\ast }}}\end{array})\end{array}$is positive definite. Accordingly, the matrix ${\mathrm{\Sigma }}_2={\rho }_2^2(P_4{P}_3{P}_2{P}_1{)}^{-1}{\mathrm{\Sigma }}_{02}[(P_4{P}_3{P}_2{P}_1{)}^{-1}{]}^T$ is also positive semi-definite and there is a positive constant ${\eta }_2$ such that ${\mathrm{\Sigma }}_2⪰{\eta }_2\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right).$Case 2. If ${\omega }_2=0$, consider the standardized transformation matrix $P_5=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& q_{32}& -q_{33}-{\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}}& q_{31}& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{C}_5& =P_5{C}_2{P}_5^{-1}\\ & =(\begin{array}{ccc}-{\alpha }_2& 0& 0\\ q_{25}{q}_{32}& -(q_{22}+q_{33}+{\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}})& c_{5(23)}\\ 0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ c_{5(24)}& q_{14}{q}_{31}& q_{33}{q}_{36}\\ 0& 0& -q_{36}\\ {\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}-q_{11}}& q_{14}& -{\displaystyle \frac{q_{12}{q}_{36}}{q_{32}}}\\ 0& -{\alpha }_1& 0\\ 0& 0& -{\alpha }_3\end{array}),\end{array}$where $\begin{array}{rl}{c}_{5(23)}& =q_{12}{q}_{21}-q_{22}(q_{33}+\frac{q_{12}{q}_{31}}{q_{32}})-q_{23}{q}_{32},\\ c_{5(24)}& =q_{31}(q_{22}-q_{11}+\frac{q_{12}{q}_{31}}{q_{32}})-q_{21}{q}_{32}.\end{array}$Then continue to do a matrix similarity transformation to $C_5$, choose the matrix $P_6=\left(\begin{array}{cccccc}{q}_{25}{q}_{32}& -(q_{22}+q_{33}+{\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}})& c_{5(23)}& c_{5(24)}& q_{14}{q}_{31}& q_{33}{q}_{36}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{C}_6& =P_6{C}_5{P}_6^{-1}\\ & =\left(\begin{array}{cccccc}-(r_1+{\alpha }_2)& -(r_1{\alpha }_2+r_2)& -r_2{\alpha }_2& r_3& c_{6(15)}& c_{6(16)}\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& -q_{36}\\ 0& 0& 0& {\displaystyle \frac{q_{12}{q}_{31}}{q_{32}}-q_{11}}& q_{14}& -{\displaystyle \frac{q_{12}{q}_{36}}{q_{32}}}\\ 0& 0& 0& 0& -{\alpha }_1& 0\\ 0& 0& 0& 0& 0& -{\alpha }_3\end{array}\right),\end{array}$where $\begin{array}{rl}{r}_1& =q_{22}+q_{33}+\frac{q_{12}{q}_{31}}{q_{32}},r_2=-c_{5(23)},r_3=({\alpha }_2+\frac{q_{12}{q}_{31}}{q_{32}}-q_{11})c_{5(24)},\\ c_{6(15)}& =q_{14}{q}_{31}(q_{22}-q_{11}+\frac{q_{12}{q}_{31}}{q_{32}})-q_{14}{q}_{21}{q}_{32}+q_{14}{q}_{31}({\alpha }_2-{\alpha }_1),\\ c_{6(16)}& =q_{33}{q}_{36}({\alpha }_2-{\alpha }_3)+q_{22}{q}_{33}{q}_{36}+q_{23}{q}_{32}{q}_{36}+\frac{q_{11}{q}_{12}{q}_{31}{q}_{36}}{q_{32}}-\frac{q_{12}^2{q}_{31}^2{q}_{36}}{q_{32}^2}.\end{array}$Furthermore, an equivalent algebraic equation of (6.6) can be described by $\begin{array}{rl}& (P_6{P}_5{P}_2{P}_1)H_2^2(P_6{P}_5{P}_2{P}_1{)}^T+C_6(P_6{P}_5{P}_2{P}_1){\mathrm{\Sigma }}_2(P_6{P}_5{P}_2{P}_1{)}^T\\ & +(P_6{P}_5{P}_2{P}_1){\mathrm{\Sigma }}_2(P_6{P}_5{P}_2{P}_1{)}^T{C}_6^T=0.\end{array}$Let ${\tilde{\mathrm{\Sigma }}}_{02}=\frac{1}{{\tilde{\rho }}_2^2}(P_6{P}_5{P}_2{P}_1){\mathrm{\Sigma }}_2(P_6{P}_5{P}_2{P}_1{)}^T,{\tilde{\rho }}_2={\sigma }_2{q}_{25}{q}_{32},$then (34) $G_0^2+C_6{\tilde{\mathrm{\Sigma }}}_{02}+{\tilde{\mathrm{\Sigma }}}_{02}{C}_6^T=0.$(34) By solving Equation (34(34) $G_0^2+C_6{\tilde{\mathrm{\Sigma }}}_{02}+{\tilde{\mathrm{\Sigma }}}_{02}{C}_6^T=0.$(34) ), we obtain $\begin{array}{rl}{\tilde{\mathrm{\Sigma }}}_{02}& =(\begin{array}{cc}{\displaystyle \frac{{\alpha }_2{r}_1+r_2}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0\\ 0& {\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}\\ -{\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\\ & \begin{array}{cccc}-{\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0& 0& 0\\ 0& 0& 0& 0\\ {\displaystyle \frac{{\alpha }_2+r_1}{2{\alpha }_2{r}_1{r}_2({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}).\end{array}$Obviously, the matrix ${\tilde{\mathrm{\Sigma }}}_{02}$ is positive semi-definite and its submatrix ${\tilde{\mathrm{\Sigma }}}_{02}^{(4)}=\left(\begin{array}{cccc}{\displaystyle \frac{{\alpha }_2{r}_1+r_2}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0& -{\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0\\ 0& {\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0& 0\\ -{\displaystyle \frac{1}{2r_1({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0& {\displaystyle \frac{{\alpha }_2+r_1}{2{\alpha }_2{r}_1{r}_2({\alpha }_2{r}_1+r_2+{\alpha }_2^2)}}& 0\\ 0& 0& 0& 0\end{array}\right)$is also positive semi-definite. So the matrix ${\mathrm{\Sigma }}_2={\tilde{\rho }}_2^2(P_6{P}_5{P}_2{P}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{02}[(P_6{P}_5{P}_2{P}_1{)}^{-1}{]}^T$ is positive semi-definite and there is a positive constant ${\tilde{\eta }}_2$ such that ${\mathrm{\Sigma }}_2⪰{\tilde{\eta }}_2\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right).$Step 3. Consider the algebraic equation (29) $H_3^2+A{\mathrm{\Sigma }}_3+{\mathrm{\Sigma }}_3{A}^T=0,$(29) where $H_3^2=\mathrm{diag}(0,0,0,0,0,{\sigma }_3^2)$.

Let $D_1=Q_{1}A{Q}_1^{-1}$, where the ordering matrix $Q_1$ takes the form $Q_1=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 1\\ 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\end{array}\right).$By a simple computation, we have $D_1=\left(\begin{array}{cccccc}-{\alpha }_3& 0& 0& 0& 0& 0\\ -q_{36}& -q_{33}& q_{32}& q_{31}& 0& 0\\ 0& -q_{23}& -q_{22}& -q_{21}& 0& q_{25}\\ 0& -q_{13}& -q_{12}& -q_{11}& q_{14}& 0\\ 0& 0& 0& 0& -{\alpha }_1& 0\\ 0& 0& 0& 0& 0& -{\alpha }_2\end{array}\right).$Choose the elimination matrix $Q_2=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& -{\displaystyle \frac{q_{13}}{q_{23}}}& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{D}_2& =Q_2{D}_1{Q}_2^{-1}\\ & =\left(\begin{array}{cccccc}-{\alpha }_3& 0& 0& 0& 0& 0\\ -q_{36}& -q_{33}& q_{32}+{\displaystyle \frac{q_{13}{q}_{31}}{q_{23}}}& q_{31}& 0& 0\\ 0& -q_{23}& -q_{22}-{\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}}& -q_{21}& 0& q_{25}\\ 0& 0& {\omega }_3& {\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}-q_{11}}& q_{14}& -{\displaystyle \frac{q_{13}{q}_{25}}{q_{23}}}\\ 0& 0& 0& 0& -{\alpha }_1& 0\\ 0& 0& 0& 0& 0& -{\alpha }_2\end{array}\right),\end{array}$where ${\omega }_3=(q_{13}{q}_{22}{q}_{23}+q_{13}^2{q}_{21}-q_{12}{q}_{23}^2-q_{11}{q}_{13}{q}_{23})/q_{23}^2$. According to the value of ${\omega }_3$, we will divide the following discussion into two parts.

Case 1. If ${\omega }_3\ne 0$, consider the standardized transformation matrix $Q_3=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& -{\omega }_3{q}_{23}& -{\omega }_3(q_{11}+q_{22})& {\displaystyle \frac{(q_{13}{q}_{21}-q_{11}{q}_{23}{)}^2}{q_{23}^2}-{\omega }_3{q}_{21}}& 0& 0\\ 0& 0& {\omega }_3& {\displaystyle \frac{q_{13}{q}_{21}-q_{11}{q}_{23}}{q_{23}}}& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{D}_3& =Q_3{D}_2{Q}_3^{-1}\\ & =(\begin{array}{cccc}-{\alpha }_3& 0& 0& 0\\ q_{23}{q}_{36}{\omega }_3& -(q_{11}+q_{22}+q_{33})& d_{3(23)}& d_{3(24)}\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\\ & \begin{array}{cc}0& 0\\ d_{3(25)}& d_{3(26)}\\ {\displaystyle \frac{q_{14}(q_{13}{q}_{21}-q_{11}{q}_{23})}{q_{23}}}& d_{3(36)}\\ q_{14}& -{\displaystyle \frac{q_{13}{q}_{25}}{q_{23}}}\\ -{\alpha }_1& 0\\ 0& -{\alpha }_2\end{array}),\end{array}$where $\begin{array}{rl}{d}_{3(23)}& =\frac{q_{13}^2{q}_{21}^2}{q_{23}^2}+\frac{q_{13}{q}_{21}(q_{22}-q_{21})}{q_{23}}-q_{11}{q}_{22}-q_{11}{q}_{33}-q_{13}{q}_{31}\\ & -q_{21}{\omega }_3-q_{22}{q}_{33}-q_{23}{q}_{32},\\ d_{3(24)}& =\frac{q_{13}{q}_{21}{q}_{33}}{q_{23}}(q_{22}+q_{13}-q_{11})+q_{13}{q}_{21}{q}_{32}-{\omega }_3(q_{21}{q}_{33}+q_{23}{q}_{31})\\ & -q_{11}{q}_{23}{q}_{32}-q_{11}{q}_{13}{q}_{31}-q_{11}{q}_{22}{q}_{33},\\ d_{3(25)}& =q_{14}(\frac{(q_{11}{q}_{23}-q_{13}{q}_{21}{)}^2}{q_{23}^2}-{\omega }_3{q}_{21}),\\ d_{3(26)}& =-q_{25}[\frac{q_{13}(q_{11}{q}_{23}-q_{13}{q}_{21}{)}^2}{q_{23}^3}+{\omega }_3(q_{11}+q_{22}-\frac{q_{13}{q}_{21}}{q_{23}})],\\ d_{3(36)}& =q_{25}[\frac{q_{13}(q_{11}{q}_{23}-q_{13}{q}_{21})}{q_{23}^2}+{\omega }_3].\end{array}$Then continue to do a matrix similarity transformation to $D_3$, choose the matrix $Q_4=\left(\begin{array}{cccccc}{q}_{23}{q}_{36}{\omega }_3& -{\tau }_1& -{\tau }_2& -{\tau }_3& d_{3(25)}& d_{3(26)}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{D}_4& =Q_4{D}_3{Q}_4^{-1}\\ & =(\begin{array}{cccc}-({\tau }_1+{\alpha }_3)& -({\tau }_1{\alpha }_3+{\tau }_2)& -({\tau }_2{\alpha }_3+{\tau }_3)& -{\tau }_3{\alpha }_3\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\\ & \begin{array}{cc}{d}_{4(15)}& d_{4(16)}\\ 0& 0\\ {\displaystyle \frac{q_{14}(q_{13}{q}_{21}-q_{11}{q}_{23})}{q_{23}}}& d_{4(36)}\\ q_{14}& -{\displaystyle \frac{q_{13}{q}_{25}}{q_{23}}}\\ -{\alpha }_1& 0\\ 0& -{\alpha }_2\end{array}),\end{array}$where $\begin{array}{rl}{d}_{4(15)}& =q_{14}\{({\alpha }_3-{\alpha }_1)[q_{12}{q}_{21}-q_{11}^2-\frac{q_{13}{q}_{21}(q_{11}+q_{22})}{q_{23}}]\\ & -{\tau }_3+{\tau }_2\frac{q_{11}{q}_{23}-q_{13}{q}_{21}}{q_{23}}\},\\ d_{4(16)}& =q_{25}\{({\alpha }_2-{\alpha }_3)[\frac{q_{13}}{q_{23}}(q_{22}^2-q_{12}{q}_{21})-q_{11}{q}_{12}-q_{12}{q}_{22}]\\ & +\frac{q_{13}}{q_{23}}{\tau }_3+\frac{{\tau }_2}{q_{23}}(q_{12}{q}_{23}-q_{13}{q}_{22})\},\\ d_{4(36)}& =\frac{q_{25}(q_{13}{q}_{22}-q_{12}{q}_{23})}{q_{23}}.\end{array}$Besides, an equivalent algebraic equation of (6.6) takes the following form $(Q_4{Q}_3{Q}_2{Q}_1)H_3^2(Q_4{Q}_3{Q}_2{Q}_1{)}^T+D_4(Q_4{Q}_3{Q}_2{Q}_1){\mathrm{\Sigma }}_3+{\mathrm{\Sigma }}_3(Q_4{Q}_3{Q}_2{Q}_1{)}^T{D}_4^T=0.$Let ${\mathrm{\Sigma }}_{03}=\frac{1}{{\rho }_3^2}(Q_4{Q}_3{Q}_2{Q}_1){\mathrm{\Sigma }}_3(Q_4{Q}_3{Q}_2{Q}_1{)}^T,$where ${\rho }_3={\sigma }_3{q}_{23}{q}_{36}{\omega }_3$, we derive (35) $G_0^2+D_4{\mathrm{\Sigma }}_{03}+{\mathrm{\Sigma }}_{03}{D}_4^T=0.$(35) By solving Equation (35(35) $G_0^2+D_4{\mathrm{\Sigma }}_{03}+{\mathrm{\Sigma }}_{03}{D}_4^T=0.$(35) ), we have $\begin{array}{rl}{\mathrm{\Sigma }}_{03}& =(\begin{array}{ccc}{\displaystyle \frac{{\alpha }_3^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_3{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0& -{\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}\\ 0& {\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0\\ -{\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0& {\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}\\ 0& -{\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ -{\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0& 0\\ 0& 0& 0\\ {\displaystyle \frac{{\alpha }_3{\tau }_1({\alpha }_3+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_3{\tau }_3{\mathrm{\Sigma }}_{03}^{\ast }}}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}),\end{array}$where ${\mathrm{\Sigma }}_{03}^{\ast }=2({\alpha }_3^3+{\alpha }_3^2{\tau }_1+{\alpha }_3{\tau }_2+{\tau }_3)({\tau }_1{\tau }_2-{\tau }_3)$.

It is easy to see that the matrix ${\mathrm{\Sigma }}_{03}$ is positive semi-definite and its submatrix $\begin{array}{rl}{\mathrm{\Sigma }}_{03}^{(4)}& =(\begin{array}{cc}{\displaystyle \frac{{\alpha }_3^2({\tau }_1{\tau }_2-{\tau }_3)+{\tau }_2({\alpha }_3{\tau }_2+{\tau }_3)}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0\\ 0& {\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}\\ -{\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0\\ 0& -{\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}\end{array}\\ & \begin{array}{cc}-{\displaystyle \frac{{\alpha }_3{\tau }_2+{\tau }_3}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0\\ 0& -{\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}\\ {\displaystyle \frac{{\alpha }_3+{\tau }_1}{{\mathrm{\Sigma }}_{03}^{\ast }}}& 0\\ 0& {\displaystyle \frac{{\alpha }_3{\tau }_1({\alpha }_3+{\tau }_1)+{\tau }_1{\tau }_2-{\tau }_3}{{\alpha }_3{\tau }_3{\mathrm{\Sigma }}_{03}^{\ast }}}\end{array})\end{array}$is positive definite. Therefore, the matrix ${\mathrm{\Sigma }}_3={\rho }_3^2(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{\mathrm{\Sigma }}_{03}[(Q_4{Q}_3{Q}_2{Q}_1{)}^{-1}{]}^T$ is also positive semi-definite and there is a positive constant ${\eta }_3$ such that ${\mathrm{\Sigma }}_3⪰{\eta }_3\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right).$Case 2. If ${\omega }_3=0$, consider the standardized transformation matrix $Q_5=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& -q_{23}& -q_{22}-{\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}}& -q_{21}& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{D}_5& =Q_5{D}_2{Q}_5^{-1}\\ & =(\begin{array}{ccc}-{\alpha }_3& 0& 0\\ q_{23}{q}_{36}& -(q_{22}+q_{33}+{\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}})& d_{5(23)}\\ 0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ & \begin{array}{ccc}0& 0& 0\\ d_{5(24)}& -q_{14}{q}_{21}& -q_{22}{q}_{25}\\ 0& 0& q_{25}\\ {\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}-q_{11}}& q_{14}& -{\displaystyle \frac{q_{13}{q}_{25}}{q_{23}}}\\ 0& -{\alpha }_1& 0\\ 0& 0& -{\alpha }_2\end{array}),\end{array}$where $\begin{array}{rl}{d}_{5(23)}& =-(q_{13}{q}_{31}+q_{22}{q}_{33}+q_{23}{q}_{32}+\frac{q_{13}{q}_{21}{q}_{33}}{q_{23}}),\\ d_{5(24)}& =q_{21}(q_{11}-q_{33}-\frac{q_{13}{q}_{21}}{q_{23}})-q_{23}{q}_{31}.\end{array}$Then continue to do a matrix similarity transformation to $D_5$, choose the matrix $Q_6=\left(\begin{array}{cccccc}{q}_{23}{q}_{36}& -(q_{22}+q_{33}+{\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}})& d_{5(23)}& d_{5(24)}& -q_{14}{q}_{21}& -q_{22}{q}_{25}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right),$such that $\begin{array}{rl}{D}_6& =Q_6{D}_5{Q}_6^{-1}\\ & =\left(\begin{array}{cccccc}-(l_1+{\alpha }_3)& -(l_1{\alpha }_3+l_2)& -l_2{\alpha }_3& l_3& d_{6(15)}& d_{6(16)}\\ 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& q_{25}\\ 0& 0& 0& {\displaystyle \frac{q_{13}{q}_{21}}{q_{23}}-q_{11}}& q_{14}& -{\displaystyle \frac{q_{13}{q}_{25}}{q_{23}}}\\ 0& 0& 0& 0& -{\alpha }_1& 0\\ 0& 0& 0& 0& 0& -{\alpha }_2\end{array}\right),\end{array}$where $\begin{array}{rl}{l}_1& =q_{22}+q_{33}+\frac{q_{13}{q}_{21}}{q_{23}},l_2=q_{13}{q}_{31}+q_{22}{q}_{33}+q_{23}{q}_{32}+\frac{q_{13}{q}_{21}{q}_{33}}{q_{23}},\\ l_3& =({\alpha }_3+\frac{q_{13}{q}_{21}}{q_{23}}-q_{11})[q_{21}(q_{11}-q_{33}-\frac{q_{13}{q}_{21}}{q_{23}})-q_{23}{q}_{31}],\\ d_{6(15)}& =q_{14}[q_{21}(q_{11}-q_{33}-\frac{q_{13}{q}_{21}}{q_{23}})-q_{23}{q}_{31}]+q_{14}{q}_{21}({\alpha }_1-{\alpha }_3),\\ d_{6(16)}& =q_{22}{q}_{25}({\alpha }_2-{\alpha }_3)-q_{25}(q_{22}{q}_{33}+q_{23}{q}_{32})-\frac{q_{13}{q}_{25}}{q_{23}}(q_{11}{q}_{21}-\frac{q_{13}{q}_{21}^2}{q_{23}}).\end{array}$Furthermore, an equivalent algebraic equation of (6.6) can be described by $\begin{array}{rl}& (Q_6{Q}_5{Q}_2{Q}_1)H_3^2(Q_6{Q}_5{Q}_2{Q}_1{)}^T+D_6(Q_6{Q}_5{Q}_2{Q}_1){\mathrm{\Sigma }}_3(Q_6{Q}_5{Q}_2{Q}_1{)}^T\\ & +(Q_6{Q}_5{Q}_2{Q}_1){\mathrm{\Sigma }}_3(Q_6{Q}_5{Q}_2{Q}_1{)}^T{D}_6^T=0.\end{array}$Denote by ${\tilde{\mathrm{\Sigma }}}_{03}=\frac{1}{{\tilde{\rho }}_3^2}(Q_6{Q}_5{Q}_2{Q}_1){\mathrm{\Sigma }}_3(Q_6{Q}_5{Q}_2{Q}_1{)}^T,{\tilde{\rho }}_3={\sigma }_3{q}_{23}{q}_{36},$then (36) $G_0^2+D_6{\tilde{\mathrm{\Sigma }}}_{03}+{\tilde{\mathrm{\Sigma }}}_{03}{D}_6^T=0.$(36) By solving Equation (36(36) $G_0^2+D_6{\tilde{\mathrm{\Sigma }}}_{03}+{\tilde{\mathrm{\Sigma }}}_{03}{D}_6^T=0.$(36) ), we obtain $\begin{array}{rl}{\tilde{\mathrm{\Sigma }}}_{03}& =(\begin{array}{cc}{\displaystyle \frac{{\alpha }_3{l}_1+l_2}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0\\ 0& {\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}\\ -{\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\\ & \begin{array}{cccc}-{\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0& 0& 0\\ 0& 0& 0& 0\\ {\displaystyle \frac{{\alpha }_3+l_1}{2{\alpha }_3{l}_1{l}_2({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}).\end{array}$It is clear that the matrix ${\tilde{\mathrm{\Sigma }}}_{03}$ is positive semi-definite and its submatrix ${\tilde{\mathrm{\Sigma }}}_{03}^{(4)}=\left(\begin{array}{cccc}{\displaystyle \frac{{\alpha }_3{l}_1+l_2}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0& -{\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0\\ 0& {\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0& 0\\ -{\displaystyle \frac{1}{2l_1({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0& {\displaystyle \frac{{\alpha }_3+l_1}{2{\alpha }_3{l}_1{l}_2({\alpha }_3{l}_1+l_2+{\alpha }_3^2)}}& 0\\ 0& 0& 0& 0\end{array}\right)$is also positive semi-definite. Hence the matrix ${\mathrm{\Sigma }}_3={\tilde{\rho }}_3^2(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{\tilde{\mathrm{\Sigma }}}_{03}[(Q_6{Q}_5{Q}_2{Q}_1{)}^{-1}{]}^T$ is positive semi-definite and there is a positive constant ${\tilde{\eta }}_3$ such that ${\mathrm{\Sigma }}_3⪰{\tilde{\eta }}_3\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right).$Now we prove the matrix $\mathrm{\Sigma }={\mathrm{\Sigma }}_1+{\mathrm{\Sigma }}_2+{\mathrm{\Sigma }}_3$ in Equation (28(28) $H^2+A\mathrm{\Sigma }+\mathrm{\Sigma }{A}^T=0.$(28) ) is positive definite. If ${\omega }_1=0$, ${\omega }_2=0$, ${\omega }_3=0$, then the covariance matrix $\begin{array}{rl}\mathrm{\Sigma }& ={\mathrm{\Sigma }}_1+{\mathrm{\Sigma }}_2+{\mathrm{\Sigma }}_3\\ & ⪰{\tilde{\eta }}_1\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)+{\tilde{\eta }}_2\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)\\ & +{\tilde{\eta }}_3\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)\\ & ⪰({\tilde{\eta }}_1\wedge {\tilde{\eta }}_2\wedge {\tilde{\eta }}_3)\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)\end{array}$is positive definite. Similarly, we can show that under the other seven cases, the covariance matrix Σ is also positive definite. This completes the proof.

7. Numerical simulations

In this section, we will present some numerical simulations to show the effectiveness of our theoretical findings. By employing the Milstein higher-order method introduced in [8], we obtain the corresponding discretization equations of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) which are given by: (37) $\{\begin{array}{l}{x}_1^{j+1}=x_1^j+x_1^j(a_1^j-b_{11}{x}_1^j-b_{12}{x}_2^j-b_{13}{x}_3^j)\mathrm{\Delta }t,x_2^{j+1}=x_2^j+x_2^j(a_2^j-b_{21}{x}_1^j-b_{22}{x}_2^j-b_{23}{x}_3^j)\mathrm{\Delta }t,x_3^{j+1}=x_3^j+x_3^j(-a_3^j+b_{31}{x}_1^j+b_{32}{x}_2^j-b_{33}{x}_3^j)\mathrm{\Delta }t,a_1^{j+1}=a_1^j+{\alpha }_1({\overline{a}}_1-a_1^j)\mathrm{\Delta }t+{\sigma }_1\sqrt{\mathrm{\Delta }t}{\epsilon }_{1,j},a_2^{j+1}=a_2^j+{\alpha }_2({\overline{a}}_2-a_2^j)\mathrm{\Delta }t+{\sigma }_2\sqrt{\mathrm{\Delta }t}{\epsilon }_{2,j},a_3^{j+1}=a_3^j+{\alpha }_3({\overline{a}}_3-a_3^j)\mathrm{\Delta }t+{\sigma }_3\sqrt{\mathrm{\Delta }t}{\epsilon }_{3,j},\end{array}$(37) where $(x_1^j,x_2^j,x_3^j,a_1^j,a_2^j,a_3^j{)}^T$ is the corresponding value of the jth iteration of the Equation (37(37) $\{\begin{array}{l}{x}_1^{j+1}=x_1^j+x_1^j(a_1^j-b_{11}{x}_1^j-b_{12}{x}_2^j-b_{13}{x}_3^j)\mathrm{\Delta }t,x_2^{j+1}=x_2^j+x_2^j(a_2^j-b_{21}{x}_1^j-b_{22}{x}_2^j-b_{23}{x}_3^j)\mathrm{\Delta }t,x_3^{j+1}=x_3^j+x_3^j(-a_3^j+b_{31}{x}_1^j+b_{32}{x}_2^j-b_{33}{x}_3^j)\mathrm{\Delta }t,a_1^{j+1}=a_1^j+{\alpha }_1({\overline{a}}_1-a_1^j)\mathrm{\Delta }t+{\sigma }_1\sqrt{\mathrm{\Delta }t}{\epsilon }_{1,j},a_2^{j+1}=a_2^j+{\alpha }_2({\overline{a}}_2-a_2^j)\mathrm{\Delta }t+{\sigma }_2\sqrt{\mathrm{\Delta }t}{\epsilon }_{2,j},a_3^{j+1}=a_3^j+{\alpha }_3({\overline{a}}_3-a_3^j)\mathrm{\Delta }t+{\sigma }_3\sqrt{\mathrm{\Delta }t}{\epsilon }_{3,j},\end{array}$(37) ). $\mathrm{\Delta }t$ is the time increment which is positive, ${\sigma }_i^2$ are the noise intensities, ${\epsilon }_{ij}$ are mutually independent Gaussian random variables which follow the distribution $\mathbb{N}(0,1)$ for $i=1,2,3;j=1,\dots ,n$. We choose real parameter values from the previous references, and all the values of biological parameters are given in Table 1.

Table 1. List of parameter values of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ).

Parameters	Values	Source
${\overline{a}}_1$	0.8	[17, 28]
${\overline{a}}_2$	0.6	[17, 28]
${\overline{a}}_3$	0.05	[17, 18, 28]
$b_{11}$	1.0	[17, 18, 28]
$b_{12}$	0.27	[17, 28]
$b_{13}$	0.4	[17, 28]
$b_{21}$	0.1	[17, 28]
$b_{22}$	2.0	[17, 18, 28]
$b_{23}$	0.2	[17, 28]
$b_{31}$	0.4	[17, 28]
$b_{32}$	0.1	[17, 28]
$b_{33}$	3.0	[17, 18, 28]
${\alpha }_1$	0.2	Estimated
${\alpha }_2$	0.1	Estimated
${\alpha }_3$	0.2	Estimated

Example 7.1

To show the existence of a stationary distribution numerically, we choose ${\sigma }_1={\sigma }_2={\sigma }_3=0.001$ and the other parameter values are given in Table 1. By simple computation, we obtain $\begin{array}{rl}& b_{12}{b}_{21}{b}_{33}+b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}-b_{11}{b}_{22}{b}_{33}-b_{11}{b}_{23}{b}_{32}-b_{13}{b}_{22}{b}_{31}=-6.2334<0,\\ & {\overline{a}}_2{b}_{12}{b}_{33}+{\overline{a}}_2{b}_{13}{b}_{32}+{\overline{a}}_3{b}_{12}{b}_{23}-{\overline{a}}_1{b}_{22}{b}_{33}-{\overline{a}}_1{b}_{23}{b}_{32}-{\overline{a}}_3{b}_{13}{b}_{22}=-4.3433<0,\\ & {\overline{a}}_1{b}_{21}{b}_{33}+{\overline{a}}_1{b}_{23}{b}_{31}+{\overline{a}}_3{b}_{13}{b}_{21}-{\overline{a}}_2{b}_{11}{b}_{33}-{\overline{a}}_2{b}_{13}{b}_{31}-{\overline{a}}_3{b}_{11}{b}_{23}=-1.6<0,\\ & {\overline{a}}_1{b}_{21}{b}_{32}+{\overline{a}}_2{b}_{12}{b}_{31}+{\overline{a}}_3{b}_{11}{b}_{22}-{\overline{a}}_1{b}_{22}{b}_{31}-{\overline{a}}_2{b}_{11}{b}_{32}-{\overline{a}}_3{b}_{12}{b}_{21}=-0.52855<0,\\ & b_{11}=1>0.0005=\frac{{\sigma }_1}{2},b_{22}=2>0.0005=\frac{{\sigma }_2}{2},b_{33}=3>0.0005=\frac{{\sigma }_3}{2},\\ & 4b_{13}{b}_{23}{b}_{31}{b}_{32}(b_{11}-\frac{{\sigma }_1}{2})(b_{22}-\frac{{\sigma }_2}{2})=0.0255808032>6.5536\times {10}^{-4}\\ & =(b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}{)}^2,\\ & \omega =\frac{b_{31}{\sigma }_1}{2b_{13}{\alpha }_1}+\frac{b_{32}{\sigma }_2}{2b_{23}{\alpha }_2}+\frac{{\sigma }_3}{2{\alpha }_3}=0.0075<0.02156952448\\ & =min\{\frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})-\frac{ϵ}{2}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}}){\overline{x}}_1^2,\frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})\\ & -\frac{1}{2ϵ}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}}){\overline{x}}_2^2,(b_{33}-\frac{{\sigma }_3}{2}){\overline{x}}_3^2\},\end{array}$and $\begin{array}{rl}\frac{b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}}{2b_{13}{b}_{32}(b_{22}-\frac{{\sigma }_2}{2})}& =0.16004001<0.2=ϵ=0.2<6.246875\\ & =\frac{2b_{23}{b}_{31}(b_{11}-\frac{{\sigma }_1}{2})}{b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}}.\end{array}$

That is to say, the conditions of Theorem 5.1 hold. According to Theorem 5.1, we get that system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) admits a stationary distribution $\pi (\cdot )$ which indicates that all the species are coexistent a.s. in the long time. See Figure 1.

Figure 1. The left column displays the trajectory of the species $x_1$, $x_2$ and $x_3$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) and their corresponding deterministic model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ) with noise intensities ${\sigma }_1={\sigma }_2={\sigma }_3=0.001$. The right column shows the frequency histograms and marginal density functions of $x_1$, $x_2$ and $x_3$ in system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ).

Example 7.2

To illustrate the existence of the probability density near the origin point $(0,0,0,0,0,0{)}^T$, we choose ${\sigma }_1={\sigma }_2={\sigma }_3=0.001$ and the other parameter values are presented in Table 1. Direct calculation gives that $(x_1^{\ast },x_2^{\ast },x_3^{\ast },{\overline{a}}_1,{\overline{a}}_2,{\overline{a}}_3{)}^T=(0.6968,0.2567,0.0848,0.8,0.6,0.05{)}^T$ and $\begin{array}{rl}& b_{12}{b}_{21}{b}_{33}+b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}-b_{11}{b}_{22}{b}_{33}-b_{11}{b}_{23}{b}_{32}-b_{13}{b}_{22}{b}_{31}=-6.2334<0,\\ & {\overline{a}}_2{b}_{12}{b}_{33}+{\overline{a}}_2{b}_{13}{b}_{32}+{\overline{a}}_3{b}_{12}{b}_{23}-{\overline{a}}_1{b}_{22}{b}_{33}-{\overline{a}}_1{b}_{23}{b}_{32}-{\overline{a}}_3{b}_{13}{b}_{22}=-4.3433<0,\\ & {\overline{a}}_1{b}_{21}{b}_{33}+{\overline{a}}_1{b}_{23}{b}_{31}+{\overline{a}}_3{b}_{13}{b}_{21}-{\overline{a}}_2{b}_{11}{b}_{33}-{\overline{a}}_2{b}_{13}{b}_{31}-{\overline{a}}_3{b}_{11}{b}_{23}=-1.6<0,\\ & {\overline{a}}_1{b}_{21}{b}_{32}+{\overline{a}}_2{b}_{12}{b}_{31}+{\overline{a}}_3{b}_{11}{b}_{22}-{\overline{a}}_1{b}_{22}{b}_{31}-{\overline{a}}_2{b}_{11}{b}_{32}-{\overline{a}}_3{b}_{12}{b}_{21}=-0.52855<0.\end{array}$In other words, the conditions of Theorem 6.1 hold. Accordingly, system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) has a Gaussian probability density near the origin point $(0,0,0,0,0,0{)}^T$. Hence, in view of Theorem 6.1, we calculate the specific form of the covariance matrix Σ, $\begin{array}{rl}\mathrm{\Sigma }& =(\begin{array}{ccc}0.00015103530& 0.03692245829& -0.00000640956\\ 0.03692245829& 9.52340771693& -0.00164333196\\ -0.00000640956& -0.00164333196& 0.00000214253\\ 0.00000784583& -0.00000003812& -0.00000017173\\ 0.00002207741& 0.00573886426& -0.00000102527\\ 0.00000023056& 0.00000007420& 0.00000046306\end{array}\\ & \begin{array}{ccc}0.00000784583& 0.00002207741& 0.00000023056\\ -0.00000003812& 0.00573886426& 0.00000007420\\ -0.00000017173& -0.00000102527& 0.00000046306\\ 0.00001377900& 0& 0\\ 0& 0.00001365688& 0\\ 0& 0& 0.00002734324\end{array}),\end{array}$and the corresponding probability density $\mathrm{\Phi }(u_1,u_2,u_3,u_4,u_5,u_6)$ takes the form $\begin{array}{rl}& \mathrm{\Phi }(u_1,u_2,u_3,u_4,u_5,u_6)\\ & =(2\pi {)}^{-3}|\mathrm{\Sigma }{|}^{-\frac{1}{2}}{\mathrm{e}}^{-\frac{1}{2}(u_1,u_2,u_3,u_4,u_5,u_6){\mathrm{\Sigma }}^{-1}(u_1,u_2,u_3,u_4,u_5,u_6{)}^T}\\ & =8390771386.35635060275{\mathrm{e}}^{-\frac{1}{2}(u_1,u_2,u_3,u_4,u_5,u_6){\mathrm{\Sigma }}^{-1}(u_1,u_2,u_3,u_4,u_5,u_6{)}^T}.\end{array}$Figure 2 shows this.

Figure 2. Numerical simulations for: (i) the frequency histogram fitting density curves of $x_1$, $x_2$, $x_3$, $a_1$, $a_2$ and $a_3$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) with 200,000 iteration points, respectively. (ii) The marginal probability densities of $x_1$, $x_2$, $x_3$, $a_1$, $a_2$ and $a_3$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ). All of the parameter values are the same as in Figure 1.

8. Conclusion

Predator–prey system is an important component of a real ecosystem. In fact, predator–prey interaction with multiple preys is common and more complicated in the nature. In this work, we concentrate on a stochastic predator–prey model with two competitive preys. In general, existing literature always simulated environmental fluctuations with white noise or coloured noise. While in this paper, we employ another way to introduce stochastic perturbations to simulate random interferences in the environment, i.e. the Ornstein–Uhlenbeck process. To ensure the mathematical and biological rationality of the stochastic system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ), we first show that the solution of the system is global. Then we study the dynamics of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) in detail. More precisely, analytical findings of this paper are presented as below:

• Given initial value $(x_1(0),x_2(0),x_3(0),a_1(0),a_2(0),a_3(0){)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$ and any integer $p\ge 3$, there exists a positive constant $M(p)$ such that the solution $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) has the property $\underset{t\to \mathrm{\infty }}{lim sup}\mathbb{E}[\frac{(\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t){)}^p}{p}+\frac{a_1^{2p}(t)}{2p}+\frac{a_2^{2p}(t)}{2p}+\frac{a_3^{2p}(t)}{2p}]\le M(p).$

• Given initial value $(x_1(0),x_2(0),x_3(0),a_1(0),a_2(0),a_3(0){)}^T\in {\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$, the solution $(x_1(t),x_2(t),x_3(t),a_1(t),a_2(t),a_3(t){)}^T$ of system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) has the property $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln (\frac{b_{31}}{b_{13}}{x}_1(t)+\frac{b_{32}}{b_{23}}{x}_2(t)+x_3(t))}{\ln t}\le \frac{1}{p}\mathrm{a}.\mathrm{s}.,\mathrm{foranyinteger}p\ge 3.$In addition $\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_1(t)}{t}\le 0,\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_2(t)}{t}\le 0\mathrm{and}\underset{t\to \mathrm{\infty }}{lim sup}\frac{\ln x_3(t)}{t}\le 0\mathrm{a}.\mathrm{s}.$

• Let Assumption (A) hold. If the following conditions hold $\begin{array}{rl}& b_{11}>\frac{{\sigma }_1}{2},b_{22}>\frac{{\sigma }_2}{2},b_{33}>\frac{{\sigma }_3}{2},\\ & 4b_{13}{b}_{23}{b}_{31}{b}_{32}(b_{11}-\frac{{\sigma }_1}{2})(b_{22}-\frac{{\sigma }_2}{2})>(b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}{)}^2,\\ & \omega <min\{\frac{b_{31}}{b_{13}}(b_{11}-\frac{{\sigma }_1}{2})-\frac{ϵ}{2}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}}){\overline{x}}_1^2,\frac{b_{32}}{b_{23}}(b_{22}-\frac{{\sigma }_2}{2})\\ & -\frac{1}{2ϵ}(\frac{b_{12}{b}_{31}}{b_{13}}+\frac{b_{21}{b}_{32}}{b_{23}}){\overline{x}}_2^2,(b_{33}-\frac{{\sigma }_3}{2}){\overline{x}}_3^2\},\end{array}$ and $ϵ>0$ is a sufficiently small constant satisfying $\frac{b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}}{2b_{13}{b}_{32}(b_{22}-\frac{{\sigma }_2}{2})}<ϵ<\frac{2b_{23}{b}_{31}(b_{11}-\frac{{\sigma }_1}{2})}{b_{12}{b}_{23}{b}_{31}+b_{13}{b}_{21}{b}_{32}},$then system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) has at least one stationary distribution $\pi (\cdot )$ on ${\mathbb{R}}_{+}^3\times {\mathbb{R}}^3$, where $\omega =b_{31}{\sigma }_1/(2b_{13}{\alpha }_1)+b_{32}{\sigma }_2/(2b_{23}{\alpha }_2)+{\sigma }_3/(2{\alpha }_3)$ and $({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3)$ is the positive equilibrium of system (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}{x}_1(t)}{\mathrm{d}t}=x_1(t)[a_1-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_2(t)}{\mathrm{d}t}=x_2(t)[a_2-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)],{\displaystyle \frac{\mathrm{d}{x}_3(t)}{\mathrm{d}t}=x_3(t)[-a_3+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)],}}}\end{array}$(1) ).

• Let $(u_1(t),u_2(t),u_3(t),u_4(t),u_5(t),u_6(t){)}^T$ be a solution of system (25(25) $\{\begin{array}{l}\mathrm{d}{u}_1=(-q_{11}{u}_1-q_{12}{u}_2-q_{13}{u}_3+q_{14}{u}_4)\mathrm{d}t,\\ \mathrm{d}{u}_2=(-q_{21}{u}_1-q_{22}{u}_2-q_{23}{u}_3+q_{25}{u}_5)\mathrm{d}t,\\ \mathrm{d}{u}_3=(q_{31}{u}_1+q_{32}{u}_2-q_{33}{u}_3-q_{36}{u}_6)\mathrm{d}t,\\ \mathrm{d}{u}_4=-{\alpha }_1{u}_4\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{u}_5=-{\alpha }_2{u}_5\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{u}_6=-{\alpha }_3{u}_6\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t),\end{array}$(25) ) with the initial value $$\begin{gathered} (u_1(0),u_2(0) \\ ,u_3(0),u_4(0),u_5(0),u_6(0){)}^T\in {\mathbb{R}}^6 \end{gathered}$$. If Assumption (A) holds, then there is a local density function of a multivariate normal distribution $\mathrm{\Phi }(u_1,u_2,u_3,u_4,u_5,u_6)$ near the origin point $(0,0,0,0,0,0{)}^T$, which takes the form $\mathrm{\Phi }(u_1,u_2,u_3,u_4,u_5,u_6)=(2\pi {)}^{-3}|\mathrm{\Sigma }{|}^{-\frac{1}{2}}{\mathrm{e}}^{-\frac{1}{2}(u_1,u_2,u_3,u_4,u_5,u_6){\mathrm{\Sigma }}^{-1}(u_1,u_2,u_3,u_4,u_5,u_6{)}^T},$where the exact form of covariance matrix Σ can be obtained in the proof process of Theorem 6.1.

However, some other topics of interest are worthy of further consideration. For example, in this paper, we only suppose that the parameters $a_i$ $(i=1,2,3)$ satisfy the Ornstein–Uhlenbeck process, as a matter of fact, if we suppose that other parameters involved in system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) satisfy the Ornstein–Uhlenbeck process, what kind of dynamic behaviour should system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) exhibit? Secondly, as far as we know, most of key parameters in population systems are usually affected by coloured noise [27]. Hence, the corresponding generalised stochastic system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) with discrete Markov switching should be studied. Finally, in this paper, we do not consider the effects of spatial diffusion factor, if we take into account the factor, i.e. we formulate a stochastic partial differential equation model, then the system (5(5) $\{\begin{array}{l}\mathrm{d}{x}_1(t)=x_1(t)[a_1(t)-b_{11}{x}_1(t)-b_{12}{x}_2(t)-b_{13}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_2(t)=x_2(t)[a_2(t)-b_{21}{x}_1(t)-b_{22}{x}_2(t)-b_{23}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{x}_3(t)=x_3(t)[-a_3(t)+b_{31}{x}_1(t)+b_{32}{x}_2(t)-b_{33}{x}_3(t)]\mathrm{d}t,\\ \mathrm{d}{a}_1(t)={\alpha }_1[{\overline{a}}_1-a_1(t)]\mathrm{d}t+{\sigma }_1\mathrm{d}{B}_1(t),\\ \mathrm{d}{a}_2(t)={\alpha }_2[{\overline{a}}_2-a_2(t)]\mathrm{d}t+{\sigma }_2\mathrm{d}{B}_2(t),\\ \mathrm{d}{a}_3(t)={\alpha }_3[{\overline{a}}_3-a_3(t)]\mathrm{d}t+{\sigma }_3\mathrm{d}{B}_3(t).\end{array}$(5) ) will show more abundant dynamic behaviours. These matters will be the subject of our future work.

 

Acknowledgments

The author read and approved the final manuscript.

Availability of data and materials

Data sharing is not applicable to this article as no datasets are generated or analysed during the current study.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The author thanks the National Natural Science Foundation of P.R. China [grant number 12001090] and the Jilin Provincial Science and Technology Development Plan Project [grant number YDZJ202201ZYTS633].