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 R0 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; Bi(t)(i=1,2,3) are mutually independent standard Brownian motions defined on a complete probability space (Ω,F,P) with a filtration {Ft}t≥0, which is increasing and right continuous while F0 contains all P-null sets; and σ2i>0(i=1,2,3) denote the intensities of the white noises. Λ is the population density entering the susceptible per unit of time, μ is the natural death rate of the total population, p is the proportion of drug users who are under treatment, β1 is the rate that an individual becomes a drug user, β2 is the rate that drug users under treatment relapsed to the untreated, δ1 is the drug-related death rate, δ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 (Ω,F,P) taking values in a finite state space S={1,2,⋯,N} for t≥0 and Δt>0, which is generated by the transition matrix Γ=(pij)N×N, i.e., P{m(t+Δt)=j|m(t)=i}≤pijΔt+o(Δt) if i≠j; otherwise, P{m(t+Δt)=j|m(t)=i}≤1+piiΔt+o(Δt) if i=j, where pij≥0 is the transition rate from state i to state j if i≠j while ∑Nj=1pij=1.
In this paper, we assume that pij>0 for i,j=1,⋯,N with i≠j. In model (1.1), the parameters Λ, p, μ, β1, β2, δ1, δ2, σi(i=1,2,3) are not constants; instead they are generated by a homogeneous continuous-time Markov chain m(t) for t≥0. That is, for each fixed k∈S, Λ(k), p(k), μ(k), β1(k), β2(k), δ1(k), δ2(k) and σ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 π=(π1,π2,⋯,πN) which can be determined by the equation πΓ=0 subject to ∑Nk=1πk=1 and πk>0 for any k∈S. Define Rn+={x∈Rn:xi>0,1≤i≤n}. For any vector g=(g(1),g(2),⋯,g(N)), let ˆg=mink∈S{g(k)} and ˇg=maxk∈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)=X0,m(0)=m, where B(⋅) and m(⋅) are the d-dimensional Brownian motions and the right-continuous Markov chains, respectively. f(⋅,⋅) and g(⋅,⋅) respectively map Rn×S to Rn and Rn×d with g(X,k)gT(X,k)=(gij(X,k))n×n. For each k∈S, let V(⋅,k) be any twice continuously differentiable function, and the operator 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))∈R3+×S, there exists a unique solution (S(t),U(t),T(t),m(t)) of model (1.1) on t≥0, and the solution will remain in R3+×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 r1≥r0 such that P{τr≤l}≥ε for each integer r≥r1. Define a C2-function V:R3+→R+ as follows:
where c is a positive constant to be determined later. Let l>1 be arbitrary, for any 0≤t≤τr∧l=min{τr,l}, and applying Itô's formula to V, we get
Choosing c such that cˇβ1+ˇβ2=ˆμ+ˆδ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 ˆμ>12(ˇσ21∨ˇσ22∨ˇσ23), 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≥0, the solution (S(t),U(t),T(t),m(t)) of model (1.1) satisfies
Proof. Define W(N)=(1+N)ρ with ρ>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 c1(k) is the solution of the linear system (2.7).
Lemma 3. For each k∈S, the linear system
has a unique solution c1=(c1(1),c1(2),⋯,c1(N))T≫0; moreover,
has a unique solution c2=(c2(1),c2(2),⋯,c2(N))T≫0.
Proof. The linear system (2.7) can be rewritten as the form of AV=β1, where V∈RN, β1=(β1(1),β1(2),⋯,β1(N))T, and
Obviously, A∈ZN×N, and ZN×N={B=(bij)N×N:bij⩽0,i≠j}. By Lemma 5.3 in [40], we obtain that determinant of (Ak) is positive for k=1,2,⋯,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 β1∈RN, the linear system (2.7) has a solution c1=(c1(1),c1(2),⋯,c1(N))T. Similarly, we show that the linear system (2.8) has a solution c2=(c2(1),c2(2),⋯,c2(N))T.
Theorem 3. If Rs0>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 pij>0 for i≠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 C2-function
such that θ∈(0,1) satisfying
and such that B>0 satisfying
here ωk will be determined later. Obviously, there exists a point (x0,y0,z0,k) at which the minimum value W(x0,y0,z0,k) is taken. We define a non-negative C2-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 R0=(R01,R02,⋯,R0N)T, since the generator matrix Γ is irreducible, there exists a solution of the Poisson system ω=(ω1,⋯,ωN)T such that
where →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 ε>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 Υ0 be a symmetric positive definite matrix, such that the three dimensional algebraic equation
holds, where G0=diag{1,0,0}, and
and also that c1>0,c3>0 and c1c2−c3>0, then Υ0 follows
Lemma 5. [42] Let Υ0 be a symmetric positive semi-definite matrix, such that the three-dimensional algebraic equation
holds, where G0=diag{1,0,0}, and
Also, d1>0,d2>0, and thus Υ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 g1,m1,m2,m3 and Δ could be found later.
Proof. Let (z1,z2,z3)T=(lnS,lnU,lnT)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∗ satisfies the following quadratic equation
with
and
Next, we discuss the value of g1 by three cases.
Case 1. If g1=0 (i.e., β1m3−β2m1=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 g1<0, then we get
Next, we turn to analyze the value of Δ. 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 g1>0, if the drug addicts under treatment T(t) has a unique positive root, the value of β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 Σ is presented later.
Proof. Let xi=zi−z∗i for i=1,2,3, and the linearized equation of model (2.24) is written as
where
Let X=(x1,x2,x3)T,B(t)=(B1(t),B2(t),B3(t))T,M=diag{σ1,σ2,σ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 Φ(X) around the positive equilibrium point P∗ 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 C0 is a positive constant, which is determined by
Also, the real symmetric inverse matrix Q meets the subsequent algebraic equation
such that Σ=Q−1, and then we derive
Furthermore, we have C0=(2π)−32|Σ|−12.
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 a211a33+a11a233>0, and direct substitution gives that
We derive that A is a Hurwitz matrix. Next, we will prove that Σ is positive definite.
Step 1. We consider the algebraic equation
and our discussion will be separated into two cases according to the value of a32.
Case 1.1. If a32≠0, according to Li et al. [45], we select the standardized transformation matrix
such that B1=H1AH−11. 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 B1 is a Hurwitz matrix. By Lemma 4, Θ1 is positive definite and takes the form
Therefore, Σ1=ϱ21H−11Θ1(HT1)−1.
Case 1.2. If a32=0, we choose ˆH1 such that ˆB1=ˆH1AˆH−11 with
where
One can equivalently transform (2.38) into
letting
The algebraic Eq (2.38) becomes
with
Therefore, Σ1=ˆϱ21ˆH−11ˆΘ1(ˆHT1)−1.
Step 2. Let us consider the following algebraic equation
we select the corresponding elimination matrix J2 and let A2=J2AJ−12 with
with
Case 2.1. If k2≠0, we then let B2=H2A2H−12 with
where
Moreover, algebraic Eq (2.43) is equivalently transformed into
and letting
algebraic Eq (2.43) is converted as
with
In other words, Σ2=ϱ22(H2J2)−1Θ2[(H2J2)T]−1.
Case 2.2. If k2=0, then we select ˆH2 and let ˆB2=ˆH2A2ˆH−12 with
where
One can equivalently transform (2.43) into
letting
which by Lemma 5 is simplified as
with
Therefore, Σ2=ˆϱ22(ˆH2J2)−1ˆΘ2[(ˆH2J2)T]−1.
Step 3. Let us consider the algebraic equation
and we select J3 and let A3=J3AJ−13 with
We find H3 such that B3=H3A3H−13 with
where
So, (2.47) is equivalently transformed into
and letting
algebraic Eq (2.47) is converted as
with
In other words, Σ3=ϱ23(H3J3)−1Θ3[(H3J3)T]−1.
3.
Examples and numerical experiments
We assume that the Markov chain m(t) takes values in the state space 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 π=(π1,π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,⋯ and
vk,i are the Gaussian random variables, which follow the standard normal distribution 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 Rs0>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 Σ 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.