Survival analysis and probability density function of switching heroin model

1.   Model establishment

Heroin is a semi-synthetic opioid drug, which is mainly extracted from opium poppy. Heroin was originally developed as a drug to cure morphine addiction, but later it was found to be highly addictive, dependence causing and toxic ^[1]. Heroin became one of the most popular drugs in the world ^[2]. White and Comiskey ^[3] were the first to study the spreading of heroin by using an ordinary differential equation (ODE) compartmental model, and they separated the local population into three compartments based on the states of drug addicts: the susceptible individuals, the untreated drug addicts and drug addicts under treatment. Based on White's model, many scholars developed different mathematical models to discuss the transmission mechanisms of heroin, such as age structure models ^[4,5,6], distributed delay models ^[7,8,9] and nonlinear incidence models ^{[10,11,12,13]} as well. Within the above-mentioned works, the authors found that the consumption of heroin was transmitted from a drug addict to a non-drug addict, which was similar to the mechanism of the spreading of infectious diseases. They further discussed the basic reproduction number $R_{0}$ as the threshold, and they determined the stability of the drug-free equilibrium and the endemic equilibrium.

Meanwhile, environmental noises usually affected the dynamics of heroin models in ^{[14,15,16,17,18,19]}. More precisely, Liu et al. ^[14] proposed a stochastic heroin epidemic model, in which they obtained a threshold for the extinction of the drug addicts. Further, ^[15] studied a stochastic heroin epidemic model with the bilinear incidence within a varying population. Then, Wei et al. ^[16] analyzed the long-term dynamics of a perturbed heroin epidemic model under non-degenerate noise. Later, Wei et al. ^[17] established a heroin population model with the standard incidence rates between distinct patches, and by constructing suitable Lyapunov functions, they established the sufficient criteria for the existence of the addict elimination and the existence of an ergodic stationary distribution. The recent contributions in ^{[20,21,22,23,24,25,26,27,28,29,30,31,32,33]} governed the continuous-time Markov chains taking values in a finite-state space to describe the regime switchings, in which Markov-chains were memoryless, and the waiting time from one state to another state usually obeyed the exponential distribution. Therefore, in this paper, we consider the following stochastic heroin model with the bilinear incidence rate under regime switching:

where $S(t)$ is the number of the susceptible individuals; $U(t)$ is the number of the untreated drug addicts; and $T(t)$ is the number of the drug addicts under treatment at time $t$ respectively. Moreover, $N(t) = S(t) + U(t) + T(t)$ denotes the total population size at time $t$; $B_{i}(t)\; (i = 1, 2, 3)$ are mutually independent standard Brownian motions defined on a complete probability space $({\mathrm{\Omega}}, \mathcal{F}, \mathbb{P})$ with a filtration $\{\mathcal{F}_{t}\}_{t \geq 0}$, which is increasing and right continuous while $\mathcal{F}_{0}$ contains all $\mathbb{P}$-null sets; and $\sigma_{i}^{2} > 0\; (i = 1, 2, 3)$ denote the intensities of the white noises. $\Lambda$ is the population density entering the susceptible per unit of time, $\mu$ is the natural death rate of the total population, $p$ is the proportion of drug users who are under treatment, $\beta_{1}$ is the rate that an individual becomes a drug user, $\beta_{2}$ is the rate that drug users under treatment relapsed to the untreated, $\delta_{1}$ is the drug-related death rate, $\delta_{2}$ is the successful cure rate. We assume that all parameters of model (1.1) are non-negative.

Let $m(t)$ be a right-continuous Markov chain on the complete probability space $({\mathrm{\Omega}}, \mathcal{F}, \mathbb{P})$ taking values in a finite state space $\mathbb{S} = \{1, 2, \cdots, N\}$ for $t \geq 0$ and $\Delta t > 0$, which is generated by the transition matrix $\Gamma = (p_{ij})_{N \times N}$, i.e., $\mathbb{P}\{m(t+\Delta t) = j|m(t) = i\}\leq p_{ij}\Delta t + o(\Delta t)$ if $i \neq j$; otherwise, $\mathbb{P}\{m(t+\Delta t) = j|m(t) = i\}\leq 1+p_{ii}\Delta t + o(\Delta t)$ if $i = j$, where $p_{ij} \geq 0$ is the transition rate from state $i$ to state $j$ if $i \neq j$ while $\sum_{j = 1}^N p_{ij} = 1.$

In this paper, we assume that $p_{ij} > 0$ for $i, j = 1, \cdots, N$ with $i \neq j$. In model (1.1), the parameters $\Lambda$, $p$, $\mu$, $\beta_{1}$, $\beta_{2}$, $\delta_{1}$, $\delta_{2}$, $\sigma_{i}\; (i = 1, 2, 3)$ are not constants; instead they are generated by a homogeneous continuous-time Markov chain $m(t)$ for $t \geq 0$. That is, for each fixed $k \in \mathbb{S}$, $\Lambda(k)$, $p(k)$, $\mu(k)$, $\beta_{1}(k)$, $\beta_{2}(k)$, $\delta_{1}(k)$, $\delta_{2}(k)$ and $\sigma_{i}(k)\; (i = 1, 2, 3)$ are all positive constants. We assume that the Markov chain $m(t)$ is irreducible, which means that the system can switch from one regime to another regime. It implies that the Markov chain $m(t)$ has a unique stationary distribution $\pi = (\pi_{1}, \pi_{2}, \cdots, \pi_{N})$ which can be determined by the equation $\pi \Gamma = 0$ subject to $\sum_{k = 1}^{N} \pi_{k} = 1$ and $\pi_{k} > 0$ for any $k \in \mathbb{S}.$ Define $\mathbb{R}_{+}^{n} = \{x\in \mathbb{R}^{n}: x_{i} > 0, 1\leq i \leq n\}$. For any vector $g = (g(1), g(2), \cdots, g(N))$, let $\hat{g} = \min_{k \in \mathbb{S}}\{g(k)\}$ and $\check{g } = \max_{k \in \mathbb{S}}\{g(k)\}$. Next, we will show the existence and uniqueness of a global positive solution. Then, we will discuss the survival dynamics including the extinction and persistence of the untreated drug addicts for the switching model (1.1). Further, we will investigate the probability density function of the degenerated model (2.20) under some sufficient conditions.

2.   Main results

In this section, we give the generalized SDEs

with the initial values $X(0) = X_{0}, m(0) = m$, where $B(\cdot)$ and $m(\cdot)$ are the $d$-dimensional Brownian motions and the right-continuous Markov chains, respectively. $f (\cdot, \cdot)$ and $g(\cdot, \cdot)$ respectively map $\mathbb{R}^{n} \times \mathbb{S}$ to $\mathbb{R}^{n}$ and $\mathbb{R}^{n\times d}$ with $g(X, k)g^{T}(X, k) = (g_{ij}(X, k))_{n\times n}$. For each $k \in \mathbb{S}$, let $V(\cdot, k)$ be any twice continuously differentiable function, and the operator $\mathcal{L}$ can be defined by

2.1.   Existence-and-uniqueness of the solution

We first of all consider the existence and uniqueness of a global positive solution before investigating other long-term properties of model (1.1) in this section.

Theorem 1. For any initial value $(S(0), U(0), T(0), m(0)) \in \mathbb{R}^{3}_{+}\times \mathbb{S}$, there exists a unique solution $(S(t), U(t), T(t), m(t))$ of model $(1.1)$ on $t\geq 0$, and the solution will remain in $\mathbb{R}^{3}_{+}\times \mathbb{S}$ with probability one.

Proof. We write down similar lines as we did in ^[34,35] and define the stopping time

Therefore, there exists an integer $r_1\geq r_0$ such that $\mathbb{P}\{\tau_r\leq l\}\geq \varepsilon$ for each integer $r\geq r_1$. Define a $C^2$-function $V: \mathbb R_+^{3}\rightarrow \mathbb R_+$ as follows:

where $c$ is a positive constant to be determined later. Let $l > 1$ be arbitrary, for any $0\leq t \leq \tau_{r}\wedge l = \min\{\tau_{r}, l\}$, and applying Itô's formula to $V$, we get

Choosing $c$ such that $c\check{\beta}_{1}+\check{\beta}_{2} = \hat{\mu}+\hat{\delta}_{1}$,

The rest of the proof is similar to Theorem 1 in ^[17], so we omit it. The proof is complete.

2.2.   Extinction of the untreated drug addicts within local population

For a long time, extinction always refers to the disappearance of infectious diseases in epidemiology. So, the most important concern of the dynamical behaviors for the stochastic heroin model is to control the spreading of heroin and the number of the untreated drug addicts. By the approaches given in ^{[34,35,36,37,38,39]}, together with constructing several Lyapunov functions, combining generalized Itô's formula and the strong law of large numbers, we derive the moderate conditions for the extinction of the untreated drug addicts to model (1.1). With these conditions, we find that the spreading of heroin ultimately vanishes in the local population, in other words, the number of the untreated drug addicts declines to zero.

Lemma 1. Assume that $\hat{\mu} > \frac{1}{2}(\check{\sigma}_{1}^2 \vee \check{\sigma}_{2}^2 \vee \check{\sigma}_{3}^2)$, and the solution $(S(t), U(t), T(t), m(t))$ of model $(1.1)$ satisfies

Proof. We write down similar lines by the same approach as in Lemma 2.2 of ^[28], so the proof is easy to check, and we omit the details.

Lemma 2. For $t \geq 0$, the solution $(S(t), U(t), T(t), m(t))$ of model $(1.1)$ satisfies

Proof. Define $W(N) = (1+N)^\rho$ with $\rho > 1$, which gives

where

The remaining proof is referred as Lemma 2.3 of ^[28]. We omit the details.

Theorem 2. If the conditions

and

hold, then

That is, the density of the untreated drug addicts will decline to zero with an exponential rate.

Proof. Model (1.1) gives

Integrating (2.2) from 0 to $t$, we get

Together with

and

the expression (2.3) further gives

From model (1.1), we can see

Combining expressions (2.4) and (2.5), we have

By the strong law of large numbers for local martingales, together with Lemma 1 and Lemma 2, we obtain

and

If

this means that

In other words, by Definition 3.2 in ^[34], the density of the untreated drug addicts declines to extinction exponentially. The proof is complete.

2.3.   Persistence of the untreated drug addicts within local population

Next, we investigate the sufficient conditions of the existence of an ergodic stationary distribution for model (1.1). Define

where

and $c_{1}(k)$ is the solution of the linear system (2.7).

Lemma 3. For each $k \in \mathbb{S}$, the linear system

has a unique solution $c_{1} = (c_{1}(1), c_{1}(2), \cdots, c_{1}(N))^{{\rm{ T}}} \gg 0$; moreover,

has a unique solution $c_{2} = (c_{2}(1), c_{2}(2), \cdots, c_{2}(N))^{{\rm{ T}}} \gg 0$.

Proof. The linear system (2.7) can be rewritten as the form of $AV = \beta_{1}$, where $V \in \mathbb{R}^ {N}$, $\beta_{1} = (\beta_{1}(1), \beta_{1}(2), \cdots, \beta_{1}(N))^{{\rm{ T}}}$, and

Obviously, $A \in Z^{N \times N}$, and $Z^{N \times N} = \big\{ B = (b_{ij})_{N \times N}: b_{ij} \leqslant 0, i \neq j \big\}$. By Lemma 5.3 in ^[40], we obtain that determinant of $(A_{k})$ is positive for $k = 1, 2, \cdots, N$, where

In other words, the leading principal minors of $A$ are all positive, which means that $A$ is a nonsingular $M$-matrix. For the vector $\beta_{1} \in \mathbb{R}^{N}$, the linear system (2.7) has a solution $c_{1} = (c_{1}(1), c_{1}(2), \cdots, c_{1}(N))^{{\rm{ T}}}$. Similarly, we show that the linear system (2.8) has a solution $c_{2} = (c_{2}(1), c_{2}(2), \cdots, c_{2}(N))^{{\rm{ T}}}$.

Theorem 3. If $R_{0}^{s} > 0,$ then model $(1.1)$ admits a unique ergodic stationary distribution.

Proof. Let

and then model (1.1) is rewritten as follows:

Equivalently, we study the stationary distribution of model (2.9) by using Lemma 2.1 in ^[41] (also referred as Lemma 5.1 in ^[36]).

Step $1$. The assumption $p_{ij} > 0$ for $i \neq j$ implies that condition (i) in Lemma 2.1 in ^[41] is satisfied.

Step $2$. The diffusion matrix

of model (2.9) is positive definite, which implies that condition (ii) in Lemma 2.1 in ^[41] holds.

Step $3$. We define a $C^2$-function

such that $\theta \in (0, 1)$ satisfying

and such that $B > 0$ satisfying

here $\omega_{k}$ will be determined later. Obviously, there exists a point $(x_{0}, y_{0}, z_{0}, k)$ at which the minimum value $W(x_{0}, y_{0}, z_{0}, k)$ is taken. We define a non-negative $C^2$-Lyapunov function as follows:

Denote

By using the generalized Itô's formula, together with the elementary equality

we obtain

and picking the coefficients by terms gives that

According to the similar discussion and Lemma 3, we obtain

with

We define a vector $R_{0} = (R_{01}, R_{02}, \cdots, R_{0N})^{{\rm{ T}}}$, since the generator matrix $\Gamma$ is irreducible, there exists a solution of the Poisson system $\omega = (\omega_{1}, \cdots, \omega_{N})^{{\rm{ T}}}$ such that

where $\vec{1}$ is a column vector in which all elements are one, which further implies

and together with (2.6), the expression (2.13) turns into

By the same arguments, we derive

Thus the following result is derived

where

Furthermore, we have

and the same arguments give that

Therefore, we take $\varepsilon > 0$ sufficiently large, and let

Then,

Hence condition (iii) of Lemma 2.1 in ^[41] is verified.

2.4.   Probability density function of model $(2.20)$ within local population

Lemma 4. ^[42] Let $\Upsilon_{0}$ be a symmetric positive definite matrix, such that the three dimensional algebraic equation

holds, where $G_{0} = diag\{1, 0, 0\}$, and

and also that $c_{1} > 0, c_{3} > 0$ and $c_{1}c_{2}-c_{3} > 0$, then $\Upsilon_{0}$ follows

Lemma 5. ^[42] Let $\Upsilon_{0}$ be a symmetric positive semi-definite matrix, such that the three-dimensional algebraic equation

holds, where $G_{0} = diag\{1, 0, 0\}$, and

Also, $d_{1} > 0, d_{2} > 0$, and thus $\Upsilon_{1}$ takes the form

Next, as a special case, we consider the degenerated model (2.20) as follows:

Then, we investigate the existence of the probability density function of model (2.20). First of all, we consider the existence of the positive equilibrium point to model (2.20).

Theorem 4. If the conditions

or

or

hold, then, model $(2.20)$ admits a positive equilibrium point $P^*$, where $g_1, m_1, m_2, m_3$ and $\Delta$ could be found later.

Proof. Let $(z_{1}, z_{2}, z_{3})^{{\rm{ T}}} = (\ln S, \ln U, \ln T)^{{\rm{ T}}}$, by using Itô's formula, the following is derived from model (2.20) that

Then, we determine the unique local equilibrium point

by solving the following equations:

in which

and $T^{\ast}$ satisfies the following quadratic equation

with

and

Next, we discuss the value of $g_{1}$ by three cases.

Case 1. If $g_{1} = 0$ (i.e., $\beta_{1}m_{3}-\beta_{2}m_{1} = 0$), and

together with

as (2.27) is valid, then we get

and thus Eq (2.26) has a unique positive root.

Case 2. If $g_{1} < 0$, then we get

Next, we turn to analyze the value of $\Delta$. When

which further gives

Thus, Eq (2.26) has a unique positive root. When

which gives that

and thus Eq (2.26) has a unique positive root.

Case 3. When $g_{1} > 0$, if the drug addicts under treatment $T(t)$ has a unique positive root, the value of $\beta_{1}$ will be very small, and the drug addicts $U(t)$ and susceptible individuals $S(t)$ are negative, so we omit this case.

Theorem 5. If the conditions of Theorem $4$ are satisfied, and

then, model $(2.20)$ possesses a probability density function

and the positive definite matrix $\Sigma$ is presented later.

Proof. Let $x_{i} = z_{i}-z_{i}^{\ast}$ for $i = 1, 2, 3$, and the linearized equation of model (2.24) is written as

where

Let $X = (x_{1}, x_{2}, x_{3})^{{\rm{ T}}}, B(t) = (B_{1}(t), B_{2}(t), B_{3}(t))^{{\rm{ T}}}, M = {\rm{diag}}\{ \sigma_{1}, \sigma_{2}, \sigma_{3}\}$ and

Therefore, Eq (2.31) can be equally rewritten as

According to the relative theory in Gardiner ^[43], there is a unique density function $\Phi(X)$ around the positive equilibrium point $P^{\ast}$ which satisfies the following equation (i.e., Fokker-Planck equation):

On the basis of Roozen ^[44], we can approximate it with a Gaussian distribution

where $C_{0}$ is a positive constant, which is determined by

Also, the real symmetric inverse matrix $Q$ meets the subsequent algebraic equation

such that $\Sigma = Q^{-1}$, and then we derive

Furthermore, we have $C_{0} = (2\pi)^{-\frac{3}{2}}|\Sigma|^{-\frac{1}{2}}.$

According to the finite independent superposition principle, we express Eq (2.33) as the sum of the solutions of the following algebraic sub-equations:

where

with

Obviously, the characteristic polynomial of matrix $A$ is

with

We find that

due to $a_{11}^{2}a_{33}+a_{11}a_{33}^{2} > 0$, and direct substitution gives that

We derive that $A$ is a Hurwitz matrix. Next, we will prove that $\Sigma$ is positive definite.

Step $1$. We consider the algebraic equation

and our discussion will be separated into two cases according to the value of $a_{32}$.

Case 1.1. If $a_{32}\neq0$, according to Li et al. ^[45], we select the standardized transformation matrix

such that $B_{1} = H_{1}AH_{1}^{-1}$. By direct calculation, we obtain

where

Furthermore, algebraic Eq (2.38) can be converted to the equivalent

letting

and algebraic Eq (2.38) is converted as

We notice that the real parts of the eigenvalues of $A$ are all negative, so $B_{1}$ is a Hurwitz matrix. By Lemma 4, $\Theta_{1}$ is positive definite and takes the form

Therefore, $\Sigma_{1} = \varrho_{1}^{2}H_{1}^{-1}\Theta_{1}(H_{1}^{T})^{-1}$.

Case 1.2. If $a_{32} = 0$, we choose $\hat{H}_{1}$ such that $\hat{B}_{1} = \hat{H}_{1}A\hat{H}_{1}^{-1}$ with

where

One can equivalently transform (2.38) into

letting

The algebraic Eq (2.38) becomes

with

Therefore, $\Sigma_{1} = \hat{\varrho}_{1}^{2}\hat{H}_{1}^{-1}\hat{\Theta}_{1}(\hat{H}_{1}^{T})^{-1}$.

Step $2$. Let us consider the following algebraic equation

we select the corresponding elimination matrix $J_{2}$ and let $A_{2} = J_{2}AJ_{2}^{-1}$ with

with

Case 2.1. If $k_{2}\neq0$, we then let $B_{2} = H_{2}A_{2}H_{2}^{-1}$ with

where

Moreover, algebraic Eq (2.43) is equivalently transformed into

and letting

algebraic Eq (2.43) is converted as

with

In other words, $\Sigma_{2} = \varrho_{2}^{2}(H_{2}J_{2})^{-1}\Theta_{2}[(H_{2}J_{2})^{T}]^{-1}$.

Case 2.2. If $k_{2} = 0$, then we select $\hat{H}_{2}$ and let $\hat{B}_{2} = \hat{H}_{2}A_{2}\hat{H}_{2}^{-1}$ with

where

One can equivalently transform (2.43) into

letting

which by Lemma 5 is simplified as

with

Therefore, $\Sigma_{2} = \hat{\varrho}_{2}^{2}(\hat{H}_{2}J_{2})^{-1}\hat{\Theta}_{2}[(\hat{H}_{2}J_{2})^{T}]^{-1}$.

Step $3$. Let us consider the algebraic equation

and we select $J_{3}$ and let $A_{3} = J_{3}AJ_{3}^{-1}$ with

We find $H_{3}$ such that $B_{3} = H_{3}A_{3}H_{3}^{-1}$ with

where

So, (2.47) is equivalently transformed into

and letting

algebraic Eq (2.47) is converted as

with

In other words, $\Sigma_{3} = \varrho_{3}^{2}(H_{3}J_{3})^{-1}\Theta_{3}[(H_{3}J_{3})^{T}]^{-1}$.

3.   Examples and numerical experiments

We assume that the Markov chain $m(t)$ takes values in the state space $\mathbb{S} = \{1, 2\}$ with the generator

The initial value is $(S(0), U(0), T(0)) = (0.70, 0.50, 0.40)$, and the unique stationary distribution of $m(t)$ is $\pi = (\pi_{1}, \pi_{2}) = (0.20, 0.80)$, respectively. We next apply two methods to simulate the sample paths of model (1).

Milstein's higher order method (MHOM). The discretization equations of model (1.1) by MHOM in ^[46] are written as follows:

Partially truncated Euler-Maruyama method (PTEMM). The PTEMM in ^[47] is written as follows:

and the discretization equations of model (1.1) are written in (3.2), so the verifications of assumptions in ^[48] are straightforward

where $i = 0, 1, 2, \cdots$ and

$v_{k, i}$ are the Gaussian random variables, which follow the standard normal distribution $\mathcal{N}(0, 1)$. Next, we use PTEMM to simulate the figures in Examples 1–3.

Example 1 We choose (Ⅰ) and (Ⅱ) of Table 1 to simulate the extinction in Theorem 2. By (Ⅰ), we obtain

By (Ⅱ), we derive

Compare the trajectories of solutions under conditions (Ⅰ) and (Ⅱ), and the time spent in Figure 2 under (Ⅱ) is shorter than that in Figure 1 under (Ⅰ) when the intensities of the white noises increase.

Example 2 We choose (Ⅲ) of Table 1 to present the results in Theorem 3. In fact, the following condition is valid:

As shown in Figure 3, the densities of the susceptible, the untreated drug addicts, and the drug addicts under treatment are stationary over time. The related simulations are demonstrated by MHOM in the middle and by PTEMM on the right. Moreover, for 50000 sample paths in total, the distributions of frequency for the solution of model (1.1) are carried out in Figure 4.

Example 3 We choose the data in Table 2 to verify the conditions of Theorem 4 and Theorem 5.

By (Ⅳ), the conditions

hold, we derive the equilibrium point $P^{*} = (1.078, 0.861, 0.407)$ by Theorem 4. Meanwhile, the stochastic persistence of density function of model (2.20) is demonstrated in Figure 5.

Or, we take parameter (Ⅴ) to compute the following conditions

then, the equilibrium point $P^{*} = (0.746, 1.670, 0.768)$ is followed. Further, the stochastic persistence of density function of model (2.20) is shown in Figure 6. Or, by selecting parameter (Ⅵ), the following conditions

hold, we obtain the positive equilibrium point $P^* = (0.967, 0.845, 0.679)$. So, the same dynamical properties appear, and we omit this case hereby.

4.   Conclusions and discussion

Heroin is an addictive drug made from the various opium poppy plants around the world. The price and spreading of heroin depend on the flowering period (usually May–July for a year) and the fruiting period (usually June–August for a year). So, we give an SUT epidemic model with regime switching to describe the flowering period and fruiting period of opium poppy plants in this paper. We are motivated by the switching between flowering period and fruiting period of opium poppy plants in years, and the recent contributions ^[17,30,31] on epidemic models. We focus on the survival analysis of switching model (1.1) and its probability density function of constant model (2.20) for investigating their long-time dynamical properties.

For the switching SUT epidemic model (1.1), the existence and uniqueness is first derived with probability one in Theorem 1 by contradiction and stochastic analysis. Further, Theorem 2, Figures 1 and 2 verify the extinction of the switching SUT model under moderate conditions in theoretical and numerical aspects. The simulations therein also reveal that the larger intensities of the white noises make the time of extinction earlier. As a consequence of theoretical investigation, we derive the important index $R_{0}^s > 0$ of the existence and uniqueness of the ergodic stationary distribution in Theorem 3. The corresponding sample paths and histogram frequencies are demonstrated in Figures 3 and 4, respectively, in which Milstein's higher order method and partially truncated Euler-Maruyama method both verify well under the same Markovian chain.

For the constant SUT epidemic model (2.20), we aim at the existence of the positive equilibrium point in Theorem 4 and the existence of probability density function in Theorem 5, respectively. One of three types of sufficient conditions is required for determining a positive equilibrium point, and details could be found in Example 3. The sample paths of model (2.20) under distinct positive equilibrium points are demonstrated in Figure 5. Further, the expression of probability density function around the positive equilibrium point is obtained in Theorem 5 after we prove that coefficient matrix $A$ is a Hurwitz matrix and diffusion matrix $\Sigma$ is positive definite by using the Fokker-Planck equation.

Acknowledgments

The research is supported by the Natural Science Foundation of Fujian Province of China (2021J01621), Special Projects of the Central Government Guiding Local Science and Technology Development (2021L3018) and Education and Research Project for Middle and Young Teachers in Fujian Province (JAT220307).

Conflict of interest

The authors declare there is no conflict of interest.