Survival analysis and probability density function of switching heroin model

: We study a switching heroin epidemic model in this paper, in which the switching of supply of heroin occurs due to the flowering period and fruiting period of opium poppy plants. Precisely, we give three equations to represent the dynamics of the susceptible, the dynamics of the untreated drug addicts and the dynamics of the drug addicts under treatment, respectively, within a local population, and the coe ffi cients of each equation are functions of Markov chains taking values in a finite state space. The first concern is to prove the existence and uniqueness of a global positive solution to the switching model. Then, the survival dynamics including the extinction and persistence of the untreated drug addicts under some moderate conditions are derived. The corresponding numerical simulations reveal that the densities of sample paths depend on regime switching, and larger intensities of the white noises yield earlier times for extinction of the untreated drug addicts. Especially, when the switching model degenerates to the constant model, we show the existence of the positive equilibrium point under moderate conditions, and we give the expression of the probability density function around the positive equilibrium point.


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: dS (t) = Λ(m(t)) − β 1 (m(t))S (t)U(t) − µ(m(t))S (t) dt +σ 1 (m(t))S (t)dB 1 (t), dU(t) = β 1 (m(t))S (t)U(t) − p(m(t))U(t) + β 2 (m(t))U(t)T (t) −(µ(m(t)) + δ 1 (m(t)))U(t) dt + σ 2 (m(t))U(t)dB 2 (t), dT (t) = p(m(t))U(t) − β 2 (m(t))U(t)T (t) − (µ(m(t)) + δ 2 (m(t)))T (t) dt +σ 3 (m(t))T (t)dB 3 (t), (1.1) 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 (Ω, F , P) with a filtration {F t } t≥0 , which is increasing and right continuous while F 0 contains all P-null sets; and σ 2 i > 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 where p i j ≥ 0 is the transition rate from state i to state j if i j while N j=1 p i j = 1. In this paper, we assume that p i j > 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 N k=1 π k = 1 and π k > 0 for any k ∈ S. Define R n For any vector g = (g(1), g(2), · · · , g(N)), letĝ = min k∈S {g(k)} anď g = max k∈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.

Main results
In this section, we give the generalized SDEs with the initial values X(0) = X 0 , 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 R n × S to R n and R n×d with g(X, k)g T (X, k) = (g i j (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
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.

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.  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. where The remaining proof is referred as Lemma 2.3 of [28]. We omit the details.
Theorem 2. If the conditionsμ and 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 (2.5) 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 In other words, by Definition 3.2 in [34], the density of the untreated drug addicts declines to extinction exponentially. The proof is complete.

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 and c 1 (k) is the solution of the linear system (2.7).
Proof. The linear system (2.7) can be rewritten as the form of Obviously, A ∈ Z N×N , and [40], we obtain that determinant of (A k ) 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 ∈ R N , the linear system (2.7) has a solution c 1 = (c 1 (1), c 1 (2), · · · , c 1 (N)) T . Similarly, we show that the linear system (2.8) has a solution c 2 = (c 2 (1), c 2 (2), · · · , c 2 (N)) T .
Step 1. The assumption p i j > 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 C 2 -function such that θ ∈ (0, 1) satisfyinĝ here ω 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: (2.10) By using the generalized Itô's formula, together with the elementary equality we obtain 11) and picking the coefficients by terms gives that (2.12) According to the similar discussion and Lemma 3, we obtain We define a vector R 0 = (R 01 , R 02 , · · · , R 0N ) 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 2 2 e (θ+1)y + 3 θΛ e θy + (β 1 +β 2 )e y + B − R s 0 + č 1 (p +δ 1 ) +č 2δ1 +β 1 +β 2 e y , Furthermore, we have and the same arguments give that Therefore, we take ε > 0 sufficiently large, and let Then, LV(x, y, z, k) ≤ −1, (x, y, z, k) ∈ U c × S.

Lemma 5. [42]
Let Υ 0 be a symmetric positive semi-definite matrix, such that the three-dimensional algebraic equation Also, d 1 > 0, d 2 > 0, and thus Υ 1 takes the form Next, as a special case, we consider the degenerated model (2.20) as follows: (2.20) 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 hold, then, model (2.20) admits a positive equilibrium point P * , where g 1 , m 1 , m 2 , m 3 and ∆ could be found later.

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 ∆. When (2.28) which further gives Thus, Eq (2.26) has a unique positive root. When 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 β 1 will be very small, and the drug addicts U(t) and susceptible individuals S (t) are negative, so we omit this case.
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): (2.32) 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 Σ = Q −1 , and then we derive According to the finite independent superposition principle, we express Eq (2.33) as the sum of the solutions of the following algebraic sub-equations: Obviously, the characteristic polynomial of matrix A is λ 3 + p 1 λ 2 + p 2 λ + p 3 = 0, with p 1 = a 33 + a 11 , p 2 = a 11 a 33 + a 12 a 21 − a 23 a 32 , p 3 = a 12 a 21 a 33 − a 11 a 23 a 32 .
(2. 37) 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 a 32 . Case 1.1. If a 32 0, according to Li et al. [45], we select the standardized transformation matrix a 32 a 21 −a 32 a 33 a 2 33 + a 32 a 23 0 a 32 −a 33 0 0 1 such that B 1 = H 1 AH −1 1 . By direct calculation, we obtain where y 1 = a 11 + a 33 , y 2 = a 11 a 33 + a 12 a 21 − a 23 a 32 , y 3 = a 12 a 21 a 33 − a 11 a 23 a 32 .
Furthermore, algebraic Eq (2.38) can be converted to the equivalent 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, Θ 1 is positive definite and takes the form One can equivalently transform (2.38) intô The algebraic Eq (2.38) becomes Step 2. Let us consider the following algebraic equation we select the corresponding elimination matrix J 2 and let A 2 = J 2 AJ −1 2 with with k 2 = a 11 a 12 a 32 − a 12 a 33 a 32 .

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 s 0 > 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.