Loading [MathJax]/jax/output/SVG/jax.js
Research article

Survival analysis and probability density function of switching heroin model


  • Received: 14 January 2023 Revised: 15 March 2023 Accepted: 03 April 2023 Published: 08 June 2023
  • 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 coefficients 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.

    Citation: Hui Jiang, Ling Chen, Fengying Wei, Quanxin Zhu. Survival analysis and probability density function of switching heroin model[J]. Mathematical Biosciences and Engineering, 2023, 20(7): 13222-13249. doi: 10.3934/mbe.2023590

    Related Papers:

    [1] Ying He, Yuting Wei, Junlong Tao, Bo Bi . Stationary distribution and probability density function analysis of a stochastic Microcystins degradation model with distributed delay. Mathematical Biosciences and Engineering, 2024, 21(1): 602-626. doi: 10.3934/mbe.2024026
    [2] Aniello Buonocore, Luigia Caputo, Enrica Pirozzi, Amelia G. Nobile . A non-autonomous stochastic predator-prey model. Mathematical Biosciences and Engineering, 2014, 11(2): 167-188. doi: 10.3934/mbe.2014.11.167
    [3] Buyu Wen, Bing Liu, Qianqian Cui . Analysis of a stochastic SIB cholera model with saturation recovery rate and Ornstein-Uhlenbeck process. Mathematical Biosciences and Engineering, 2023, 20(7): 11644-11655. doi: 10.3934/mbe.2023517
    [4] Yan Xie, Zhijun Liu, Ke Qi, Dongchen Shangguan, Qinglong Wang . A stochastic mussel-algae model under regime switching. Mathematical Biosciences and Engineering, 2022, 19(5): 4794-4811. doi: 10.3934/mbe.2022224
    [5] Pierre Degond, Maxime Herda, Sepideh Mirrahimi . A Fokker-Planck approach to the study of robustness in gene expression. Mathematical Biosciences and Engineering, 2020, 17(6): 6459-6486. doi: 10.3934/mbe.2020338
    [6] Xi-Chao Duan, Xue-Zhi Li, Maia Martcheva . Dynamics of an age-structured heroin transmission model with vaccination and treatment. Mathematical Biosciences and Engineering, 2019, 16(1): 397-420. doi: 10.3934/mbe.2019019
    [7] Huili Wei, Wenhe Li . Dynamical behaviors of a Lotka-Volterra competition system with the Ornstein-Uhlenbeck process. Mathematical Biosciences and Engineering, 2023, 20(5): 7882-7904. doi: 10.3934/mbe.2023341
    [8] Zhongcai Zhu, Xiaomei Feng, Xue He, Hongpeng Guo . Mirrored dynamics of a wild mosquito population suppression model with Ricker-type survival probability and time delay. Mathematical Biosciences and Engineering, 2024, 21(2): 1884-1898. doi: 10.3934/mbe.2024083
    [9] Yijun Lou, Bei Sun . Stage duration distributions and intraspecific competition: a review of continuous stage-structured models. Mathematical Biosciences and Engineering, 2022, 19(8): 7543-7569. doi: 10.3934/mbe.2022355
    [10] H.Thomas Banks, Shuhua Hu . Nonlinear stochastic Markov processes and modeling uncertainty in populations. Mathematical Biosciences and Engineering, 2012, 9(1): 1-25. doi: 10.3934/mbe.2012.9.1
  • 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 coefficients 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.



    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:

    {dS(t)=[Λ(m(t))β1(m(t))S(t)U(t)μ(m(t))S(t)]dt+σ1(m(t))S(t)dB1(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)dB2(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)dB3(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; Bi(t)(i=1,2,3) are mutually independent standard Brownian motions defined on a complete probability space (Ω,F,P) with a filtration {Ft}t0, 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 t0 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 ij; otherwise, P{m(t+Δt)=j|m(t)=i}1+piiΔt+o(Δt) if i=j, where pij0 is the transition rate from state i to state j if ij while Nj=1pij=1.

    In this paper, we assume that pij>0 for i,j=1,,N with ij. 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 t0. That is, for each fixed kS, Λ(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 kS. Define Rn+={xRn:xi>0,1in}. For any vector g=(g(1),g(2),,g(N)), let ˆg=minkS{g(k)} and ˇg=maxkS{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.

    In this section, we give the generalized SDEs

    dX(t)=f(X(t),m(t))dt+g(X(t),m(t))dB(t),t0, (2.1)

    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 kS, let V(,k) be any twice continuously differentiable function, and the operator L can be defined by

    LV(X,k)=Ni=1fi(X,k)V(X,k)Xi+12Ni,j=1gij(X,k)2V(X,k)XiXj+lNpklV(X,l).

    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 t0, 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

    τr=inf{t[0,τe):min{S(t),U(t),T(t)}1rormax{S(t),U(t),T(t)}r}.

    Therefore, there exists an integer r1r0 such that P{τrl}ε for each integer rr1. Define a C2-function V:R3+R+ as follows:

    V(S,U,T)=SclnSc+U1lnU+T1lnT,

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

    LV=μ(k)S(μ(k)+δ1(k))U(μ(k)+δ2(k))TcΛ(k)Sβ1(k)Sβ2(k)Tp(k)UT+(cβ1(k)+β2(k))U+Λ(k)+(c+2)μ(k)+p(k)+δ1(k)+δ2(k)+12(cσ21(k)+σ22(k)+σ23(k))<ˇΛ+(c+2)ˇμ+ˇp+ˇδ1+ˇδ2+12(cˇσ21+ˇσ22+ˇσ23)+(cˇβ1+ˇβ2ˆμˆδ1)U.

    Choosing c such that cˇβ1+ˇβ2=ˆμ+ˆδ1,

    LVˇΛ+(c+2)ˇμ+ˇp+ˇδ1+ˇδ2+12(cˇσ21+ˇσ22+ˇσ23):=K.

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

    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

    limt1tt0ˇσ1S(s)dB1(s)=limt1tt0ˇσ2U(s)dB2(s)=0a.s.,
    limt1tt0ˇσ3T(s)dB3(s)=0a.s..

    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 t0, the solution (S(t),U(t),T(t),m(t)) of model (1.1) satisfies

    limt1t(S(t)+U(t)+T(t))=0a.s..

    Proof. Define W(N)=(1+N)ρ with ρ>1, which gives

    dW(N)=LW(N)dt+ρ(1+N)ρ1(σ1(k)SdB1(t)+σ2(k)UdB2(t)+σ3(k)TdB3(t)),

    where

    LW(N)=ρ(1+N)ρ1[Λ(k)μ(k)S(μ(k)+δ1(k))U(μ(k)+δ2(k))T]+ρ(ρ1)2(1+N)ρ2(σ21(k)S2+σ22(k)U2+σ23(k)T2)=ρ(1+N)ρ2[(1+N)(Λ(k)μ(k)S(μ(k)+δ1(k))U(μ(k)+δ2(k))T)+ρ12(σ21(k)S2+σ22(k)U2+σ23(k)T2)]ρ(1+N)ρ2[(1+N)(ˇΛˆμN)+ρ12(ˇσ21ˇσ22ˇσ23)N2]=ρ(1+N)ρ2[(ˆμρ12(ˇσ21ˇσ22ˇσ23))N2+(ˇΛˆμ)N+ˇΛ].

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

    Theorem 2. If the conditions

    ˆμ>12(ˇσ21ˇσ22ˇσ23)

    and

    max{ˇβ1,ˇβ2}ˆμˇΛ(ˆp+ˆμ+ˆδ1+12ˆσ22)<0

    hold, then

    limtlnU(t)t<0.

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

    Proof. Model (1.1) gives

    d(S+U+T)<[ˇΛˆμ(S+T)(ˆμ+ˆδ1)U]dt+ˇσ1SdB1(t)+ˇσ2UdB2(t)+ˇσ3TdB3(t). (2.2)

    Integrating (2.2) from 0 to t, we get

    A(t)t1tt0[ˇΛˆμ(S(s)+T(s))(ˆμ+ˆδ1)U(s)]ds+M(t)t. (2.3)

    Together with

    M(t)=t0ˇσ1S(s)dB1(s)+t0ˇσ2U(s)dB2(s)+t0ˇσ3T(s)dB3(s)

    and

    A(t)=S(t)+U(t)+T(t)S(0)U(0)T(0),

    the expression (2.3) further gives

    1tt0(S(s)+T(s))ds1ˆμ{1tt0[ˇΛ(ˆμ+ˆδ1)U(s)]ds+M(t)tA(t)t}. (2.4)

    From model (1.1), we can see

    1t(lnU(t)lnU(0))max{ˇβ1,ˇβ2}1tt0(S(s)+T(s))ds(ˆp+ˆμ+ˆδ1+12ˆσ22)+ˇσ2B2(t)t. (2.5)

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

    lnU(t)tmax{ˇβ1,ˇβ2}ˆμ{1tt0[ˇΛ(ˆμ+ˆδ1)U(s)]ds+M(t)tA(t)t}(ˆp+ˆμ+ˆδ1+12ˆσ22)+ˇσ2B2(t)t+lnU(0)t<max{ˇβ1,ˇβ2}ˆμ(M(t)tA(t)t)+ˇσ2B2(t)t+lnU(0)t+ˇΛmax{ˇβ1,ˇβ2}ˆμ(ˆp+ˆμ+ˆδ1+12ˆσ22).

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

    limt1tt0ˇσ2dB2(s)=0a.s.,

    and

    limt1t{max{ˇβ1,ˇβ2}ˆμ(M(t)A(t))+lnU(0)}=0a.s..

    If

    max{ˇβ1,ˇβ2}ˆμˇΛ(ˆp+ˆμ+ˆδ1+12ˆσ22)<0,

    this means that

    limtlnU(t)t<0.

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

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

    Rs0:=kSπkR0k, (2.6)

    where

    R0k=c1(k)Λ(k)p(k)μ(k)δ1(k)12σ22(k),

    and c1(k) is the solution of the linear system (2.7).

    Lemma 3. For each kS, the linear system

    c1(k)μ(k)β1(k)lSpklc1(l)=0 (2.7)

    has a unique solution c1=(c1(1),c1(2),,c1(N))T0; moreover,

    c2(k)(μ(k)+δ2(k))β2(k)lSpklc2(l)=0 (2.8)

    has a unique solution c2=(c2(1),c2(2),,c2(N))T0.

    Proof. The linear system (2.7) can be rewritten as the form of AV=β1, where VRN, β1=(β1(1),β1(2),,β1(N))T, and

    A=(μ(1)p11p12p1Np21μ(2)p22p2NpN1pN2μ(N)pNN).

    Obviously, AZN×N, and ZN×N={B=(bij)N×N:bij0,ij}. By Lemma 5.3 in [40], we obtain that determinant of (Ak) is positive for k=1,2,,N, where

    Ak=(μ(1)p11p12p1kp21μ(2)p22p2kpk1pk2μ(k)pkk).

    In other words, the leading principal minors of A are all positive, which means that A is a nonsingular M-matrix. For the vector β1RN, 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

    x(t)=lnS(t),y(t)=lnU(t),z(t)=lnT(t),

    and then model (1.1) is rewritten as follows:

    {dx(t)=[Λ(m(t))exβ1(m(t))ey(μ(m(t))+12σ21(m(t)))]dt+σ1(m(t))dB1(t),dy(t)=[β1(m(t))exp(m(t))+β2(m(t))ez(μ(m(t))+δ1(m(t))+12σ22(m(t)))]dt+σ2(m(t))dB2(t),dz(t)=[p(m(t))eyezβ2(m(t))ey(μ(m(t))+δ2(m(t))+12σ23(m(t)))]dt+σ3(m(t))dB3(t). (2.9)

    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 ij implies that condition (i) in Lemma 2.1 in [41] is satisfied.

    Step 2. The diffusion matrix

    D(x,k)=(σ21(k)000σ22(k)000σ23(k))

    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

    W(x,y,z,k)=1θ+1(ex+ey+ez)(θ+1)B[c1(k)(ex+ey)+c2(k)(ey+ez)+y+ωk]xz,

    such that θ(0,1) satisfying

    ˆμ0.5θˇσ21>0,ˆμ+ˆδ10.5θˇσ22>0,ˆμ+ˆδ20.5θˇσ23>0,

    and such that B>0 satisfying

    fu1+fu3BRs02,

    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:

    V(x,y,z,k)=W(x,y,z,k)W(x0,y0,z0,k). (2.10)

    Denote

    V1(x,y,z,k)=1θ+1(ex+ey+ez)(θ+1),V2(x,y,z,k)=c1(k)(ex+ey)c2(k)(ey+ez)yωk,V3(y,k)=x,V4(x,k)=z.

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

    (a+b+c)θ3θ(aθ+bθ+cθ),fora>0,b>0,c>0,

    we obtain

    LV1=(ex+ey+ez)θ[Λ(k)μ(k)ex(μ(k)+δ1(k))ey(μ(k)+δ2(k))ez]+θ2(ex+ey+ez)θ1(σ21(k)e2x+σ22(k)e2y+σ23(k)e2z)3θΛ(k)(eθx+eθy+eθz)+θ2(σ21(k)e(θ+1)x+σ22(k)e(θ+1)y+σ23(k)e(θ+1)z)(μ(k)+δ2(k))e(θ+1)zμ(k)e(θ+1)x(μ(k)+δ1(k))e(θ+1)y, (2.11)

    and picking the coefficients by terms gives that

    LV1(ˆμθ2ˇσ21)e(θ+1)x(ˆμ+ˆδ1θ2ˇσ22)e(θ+1)y(ˆμ+ˆδ2θ2ˇσ23)e(θ+1)z+3θˇΛ(eθx+eθy+eθz). (2.12)

    According to the similar discussion and Lemma 3, we obtain

    LV2=c1(k)[Λ(k)μ(k)exp(k)ey+β2(k)ey+z(μ(k)+δ1(k))ey]c2(k)[β1(k)ex+y(μ(k)+δ1(k))ey(μ(k)+δ2(k))ez]β1(k)exβ2(k)ez+p(k)+μ(k)+δ1(k)+12σ22(k)lSpklω(l)lSpklc1(l)(ex+ey)lSpklc2(l)(ey+ez)[c1(k)μ(k)β1(k)lSpklc1(l)]ex+[c2(k)(μ(k)+δ2(k))β2(k)lSpklc2(l)]ez+[c1(k)(p(k)+μ(k)+δ1(k))+c2(k)(μ(k)+δ1(k))lSpklc1(l)lSpklc2(l)]ey+p(k)+μ(k)+δ1(k)+12σ22(k)c1(k)Λ(k)lSpklω(l)=:R0klSpklω(l)+[c1(k)(p(k)+δ1(k))+c2(k)(δ1(k)δ2(k))+β1(k)+β2(k)]ey, (2.13)

    with

    R0k=c1(k)Λ(k)p(k)μ(k)δ1(k)12σ22(k).

    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

    Γω=(Nk=1πkR0k)1R0, (2.14)

    where 1 is a column vector in which all elements are one, which further implies

    R0k+lSpklω(l)=Nk=1πkR0k,

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

    LV2Rs0+[ˇc1(ˇp+ˇδ1)+ˇc2ˇδ1+ˇβ1+ˇβ2]ey. (2.15)

    By the same arguments, we derive

    LV3=Λ(k)ex+β1(k)ey+μ(k)+12σ21(k)ˆΛex+ˇβ1ey+ˇμ+12ˇσ21, (2.16)
    LV4ˆp(k)ez+ˇβ2(k)ey+ˇμ+ˇδ2+12ˇσ23. (2.17)

    Thus the following result is derived

    LV=LV1+BLV2+LV3+LV4<f(x,y,z)=f1(x)+f2(y)+f3(z), (2.18)

    where

    f1(x)=(ˆμθ2ˇσ21)e(θ+1)x+3θˇΛeθxˆΛex+2ˇμ+ˇδ2+12(ˇσ21+ˇσ23),f2(y)=(ˆμ+ˆδ1θ2ˇσ22)e(θ+1)y+3θˇΛeθy+(ˇβ1+ˇβ2)ey+B[Rs0+(ˇc1(ˇp+ˇδ1)+ˇc2ˇδ1+ˇβ1+ˇβ2)ey],f3(z)=(ˆμ+ˆδ2θ2ˇσ23)e(θ+1)z+3θˇΛeθzˆpez.

    Furthermore, we have

    f(x,y,z)f1(x)+fu2+fu3,ifx+orx,

    and the same arguments give that

    f(x,y,z)fu1+f2(y)+fu3,ify+,f(x,y,z)fu1+f2(y)+fu3fu1+fu3BRs02,ify,f(x,y,z)fu1+fu2+f3(z),ifz+orz.

    Therefore, we take ε>0 sufficiently large, and let

    U=(ε,ε)×(ε,ε)×(ε,ε)×(ε,ε).

    Then,

    LV(x,y,z,k)1,(x,y,z,k)Uc×S.

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

    Lemma 4. [42] Let Υ0 be a symmetric positive definite matrix, such that the three dimensional algebraic equation

    G20+A0Υ0+Υ0AT0=0 (2.19)

    holds, where G0=diag{1,0,0}, and

    A0=(c1c2c3100010)

    and also that c1>0,c3>0 and c1c2c3>0, then Υ0 follows

    Υ0=12(c1c2c3)(c20101010c1c3).

    Lemma 5. [42] Let Υ0 be a symmetric positive semi-definite matrix, such that the three-dimensional algebraic equation

    G20+˜A0Υ1+Υ1˜AT0=0

    holds, where G0=diag{1,0,0}, and

    ˜A0=(d1d2d310000d33.).

    Also, d1>0,d2>0, and thus Υ1 takes the form

    Υ1=diag{12d1,12d1d2,0}.

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

    {dS(t)=[Λβ1S(t)U(t)μS(t)]dt+σ1S(t)dB1(t),dU(t)=[β1S(t)U(t)pU(t)+β2U(t)T(t)(μ+δ1)U(t)]dt+σ2U(t)dB2(t),dT(t)=[pU(t)β2U(t)T(t)(μ+δ2)T(t)]dt+σ3T(t)dB3(t). (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

    g1=0,1+m1m2m1pβ1Λ>0, (2.21)

    or

    g1<0,Δ>0,β1Λm1pm1m2>0, (2.22)

    or

    g1<0,Δ=0,(β2m1β1m3)(m2+p)+β2(m1pβ1Λ)>0, (2.23)

    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

    {dz1=[Λez1β1ez2(μ+σ212)]dt+σ1dB1(t),dz2=[β1ez1p+β2ez3(μ+δ1+σ222)]dt+σ2dB2(t),dz3=[pez2z3β2ez2(μ+δ2+σ232)]dt+σ3dB3(t). (2.24)

    Then, we determine the unique local equilibrium point

    P=(S,U,T)=(ez1,ez2,ez3),

    by solving the following equations:

    {Λez1β1ez2(μ+σ212)=0,β1ez1p+β2ez3(μ+δ1+σ222)=0,pez2z3β2ez2(μ+δ2+σ232)=0, (2.25)

    in which

    S=pβ2T+m2β1>0,U=m3Tpβ2T>0,

    and T satisfies the following quadratic equation

    g1T2+g2T+g3=0, (2.26)

    with

    g1=β1β2m3β22m1,g2=2β2m1p+β2m1m2β1β2Λβ1m2m3β1m3p,g3=β1pΛm1m2pm1p2,

    and

    m1=μ+σ212,m2=μ+δ1+σ222,m3=μ+δ2+σ232.

    Next, we discuss the value of g1 by three cases.

    Case 1. If g1=0 (i.e., β1m3β2m1=0), and

    1+m1m2m1pβ1Λ>0, (2.27)

    together with

    g2=β2(m1pβ1Λ),g3=p(β1Λm1m2m1p),

    as (2.27) is valid, then we get

    g3g2=p(β1Λm1m2m1p)β2(m1pβ1Λ)=pβ2(1+m1m2m1pβ1Λ)<0,

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

    Case 2. If g1<0, then we get

    Δ=g224g1g3=[2β2m1p+β2m1m2β1β2Λβ1m2m3β1m3p]2+4β2p(β2m1β1m3)(β1Λm1pm1m2).

    Next, we turn to analyze the value of Δ. When

    Δ>0,β1Λm1pm1m2>0, (2.28)

    which further gives

    g3g1=β1pΛm1m2pm1p2β1β2m3β22m1=p(β1Λm1m2m1p)β2(β1m3β2m1)<0.

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

    Δ=0,(β2m1β1m3)(m2+p)+β2(m1pβ1Λ)>0, (2.29)

    which gives that

    g22g1=2β2m1p+β2m1m2β1β2Λβ1m2m3β1m3p2β2(β1m3β2m1)=(β2m1β1m3)(m2+p)+β2(m1pβ1Λ)2β2(β2m1β1m3)>0,

    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

    Λβ21m3β2p>0, (2.30)

    then, model (2.20) possesses a probability density function

    Φ(S,U,T)=(2π)32|Σ|12e12(lnSS,lnUU,lnTT)Σ1(lnSS,lnUU,lnTT)T,

    and the positive definite matrix Σ is presented later.

    Proof. Let xi=zizi for i=1,2,3, and the linearized equation of model (2.24) is written as

    {dx1=(a11x1a12x2+a13x3)dt+σ1dB1(t),dx2=(a21x1+a22x2+a23x3)dt+σ2dB2(t),dx3=(a31x1+a32x2a33x3)dt+σ3dB3(t), (2.31)

    where

    a11=Λez1,a12=β1ez2,a13=0,a21=β1ez1,a22=0,a23=β2ez3,a31=0,a32=pez2z3β2ez2,a33=pez2z3.

    Let X=(x1,x2,x3)T,B(t)=(B1(t),B2(t),B3(t))T,M=diag{σ1,σ2,σ3} and

    A=(a11a120a210a230a32a33).

    Therefore, Eq (2.31) can be equally rewritten as

    dX(t)=AX(t)dt+MdB(t).

    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):

    3r=1σ2i22Φx2i+x1[(a11x1a12x2+a13x3)Φ]+x2[(a21x1+a22x2+a23x3)Φ]+x3[a31x1+a32x2a33x3)Φ]=0. (2.32)

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

    Φ(X)=Φ(x1,x2,x3)=C0e12(x1,x2,x3)Q(x1,x2,x3)T,

    where C0 is a positive constant, which is determined by

    R3Φ(x1,x2,x3)dx1dx2dx3=1.

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

    QM2Q+QA+ATQ=0,

    such that Σ=Q1, and then we derive

    M2+AΣ+ΣAT=0. (2.33)

    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:

    M2k+AΣk+ΣkAT=0,k=1,2,3, (2.34)

    where

    M1=diag(σ1,0,0),M2=diag(0,σ2,0),M3=diag(0,0,σ3)

    with

    Σ=Σ1+Σ2+Σ3,M2=M21+M22+M23.

    Obviously, the characteristic polynomial of matrix A is

    λ3+p1λ2+p2λ+p3=0,

    with

    p1=a33+a11,p2=a11a33+a12a21a23a32,p3=a12a21a33a11a23a32. (2.35)

    We find that

    Δ1=p1=a33+a11>0,Δ2=p1p2p3=a211a33+a11a12a21+a11a233a23a32a33, (2.36)

    due to a211a33+a11a233>0, and direct substitution gives that

    a11a12a21a23a32a33=m3T(Λβ21β2m3p)pβ2T>0. (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

    M21+AΣ1+Σ1AT=0, (2.38)

    and our discussion will be separated into two cases according to the value of a32.

    Case 1.1. If a320, according to Li et al. [45], we select the standardized transformation matrix

    H1=(a32a21a32a33a233+a32a230a32a33001), (2.39)

    such that B1=H1AH11. By direct calculation, we obtain

    B1=(y1y2y3100010),

    where

    y1=a11+a33,y2=a11a33+a12a21a23a32,y3=a12a21a33a11a23a32.

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

    H1M21HT1+B1H1Σ1HT1+H1Σ1HT1BT1=0,

    letting

    Θ1=ϱ21H1Σ1HT1,ϱ1=a21a32σ1,

    and algebraic Eq (2.38) is converted as

    G20+B1Θ1+Θ1BT1=0. (2.40)

    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

    Θ1=12(y1y2y3)(y20101010y1y3).

    Therefore, Σ1=ϱ21H11Θ1(HT1)1.

    Case 1.2. If a32=0, we choose ˆH1 such that ˆB1=ˆH1AˆH11 with

    ˆH1=(a210a23010001),ˆB1=(b1b2b310000a33),

    where

    b1=a11,b2=a12a21,b3=a23a33a11a23.

    One can equivalently transform (2.38) into

    ˆH1M21ˆHT1+ˆB1ˆH1Σ1ˆHT1+ˆH1Σ1ˆHT1ˆBT1=0,

    letting

    ˆΘ1=ˆϱ21ˆH1Σ1ˆHT1,ˆϱ1=a21σ1.

    The algebraic Eq (2.38) becomes

    G20+ˆB1ˆΘ1+ˆΘ1ˆBT1=0, (2.41)

    with

    ˆΘ1=diag{12b1,12b1b2,0}. (2.42)

    Therefore, Σ1=ˆϱ21ˆH11ˆΘ1(ˆHT1)1.

    Step 2. Let us consider the following algebraic equation

    M22+AΣ2+Σ2AT=0, (2.43)

    we select the corresponding elimination matrix J2 and let A2=J2AJ12 with

    J2=(01000110a12a32),A2=(0a12a21a32+a23a21a32a3300k2a11),

    with

    k2=a11a12a32a12a33a32.

    Case 2.1. If k20, we then let B2=H2A2H12 with

    H2=(k2a32k2(a11+a33)a2110k2a11001),B2=(q1q2q3100010), (2.44)

    where

    q1=a11+a33,q2=a11a33+a12a21a23a32,q3=a12a21a33a11a23a32.

    Moreover, algebraic Eq (2.43) is equivalently transformed into

    (H2J2)M22(H2J2)T+B2[(H2J2)Σ2(H2J2)T]+[(H2J2)Σ2(H2J2)T]BT2=0,

    and letting

    Θ2=ϱ22(H2J2)Σ2(H2J2)T,ϱ2=k2a32σ2,

    algebraic Eq (2.43) is converted as

    G20+B2Θ2+Θ2BT2=0, (2.45)

    with

    Θ2=12(q1q2q3)(q20101010q1q3).

    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ˆH12 with

    ˆH2=(a32a330010001),ˆB2=(ω1ω2ω310000a11),

    where

    ω1=a33,ω2=a12a21a23a32,ω3=a21a32.

    One can equivalently transform (2.43) into

    (ˆH2J2)M22(ˆH2J2)T+ˆB2[(ˆH2J2)Σ2(ˆH2J2)T]+[(ˆH2J2)Σ2(ˆH2J2)T]ˆBT2=0,

    letting

    ˆΘ2=ˆϱ22(ˆH2J2)Σ2(ˆH2J2)T,ˆϱ2=a32σ2,

    which by Lemma 5 is simplified as

    G20+ˆB2ˆΘ2+ˆΘ2ˆBT2=0, (2.46)

    with

    ˆΘ2=diag{12ω1,12ω1ω2,0}.

    Therefore, Σ2=ˆϱ22(ˆH2J2)1ˆΘ2[(ˆH2J2)T]1.

    Step 3. Let us consider the algebraic equation

    M23+AΣ3+Σ3AT=0, (2.47)

    and we select J3 and let A3=J3AJ13 with

    J3=(001010100),A3=(a33a320a230a210a12a11).

    We find H3 such that B3=H3A3H13 with

    H3=(a12a23a11a12a211a12a210a12a11001),B3=(s1s2s3100010), (2.48)

    where

    s1=a11+a33,s2=a11a33+a12a21a23a32,s3=a12a21a33a11a23a32.

    So, (2.47) is equivalently transformed into

    (H3J3)M23(H3J3)T+B3[(H3J3)Σ3(H3J3)T]+[(H3J3)Σ3(H3J3)T]BT3=0,

    and letting

    Θ3=ϱ23(H3J3)Σ3(H3J3)T,ϱ3=a12a23σ3,

    algebraic Eq (2.47) is converted as

    G20+B3Θ3+Θ3BT3=0, (2.49)

    with

    Θ3=12(s1s2s3)(s20101010s1s3).

    In other words, Σ3=ϱ23(H3J3)1Θ3[(H3J3)T]1.

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

    Γ=(0.800.800.200.20).

    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:

    Si+1=Si+(Λ(k)β1(k)SiUiμ(k)Si)Δt+σ1(k)SiΔtvk,i+0.5σ21(k)S2i(v2k,i1)Δt,Ui+1=Ui+(β1(k)SiUi+β2(k)UiTi(μ(k)+p(k)+δ1(k))Ui)Δt+σ2(k)UiΔtvk,i+0.5σ22(k)U2i(v2k,i1)Δt,Ti+1=Ti+(p(k)β2(k)UiTi(μ(k)+δ2(k))Ti)Δt+σ3(k)TiΔtvk,i+0.5σ23(k)T2i(v2k,i1)Δt,i=0,1,2,. (3.1)

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

    Δt=104,h(Δt)=Δt13,u1(r)=r3,Xc=3×10233S2i+U2i+T2i,

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

    Si+1=Si+(Λ(k)f1Δt,iμ(k)Si)Δt+g1Δt,iΔtvk,i,Ui+1=Ui+(f2Δt,i(μ(k)+p(k)+δ1(k))Ui)Δt+g2Δt,iΔtvk,i,Ti+1=Ti+(p(k)f3Δt,i(μ(k)+δ2(k))Ti)Δt+g3Δt,iΔtvk,i, (3.2)

    where i=0,1,2, and

    f1Δt,i=(1Xc)β1(k)SiUi,f2Δt,i=(1Xc)(β1(k)SiUi+β2(k)UiTi),f3Δt,i=(1Xc)β2(k)UiTi,g1Δt,i=(1Xc)σ1(k)Si,g2Δt,i=(1Xc)σ2(k)Ui,g3Δt,i=(1Xc)σ3(k)Ti.

    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

    ˆμ=0.25>0.225=12(ˇσ21ˇσ22ˇσ23),max{ˇβ1,ˇβ2}ˆμˇΛ(ˆp+ˆμ+ˆδ1+12ˆσ22)=0.235<0.

    By (Ⅱ), we derive

    ˆμ=0.25>0.245=12(ˇσ21ˇσ22ˇσ23),max{ˇβ1,ˇβ2}ˆμˇΛ(ˆp+ˆμ+ˆδ1+12ˆσ22)=0.260<0.

    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.

    Table 1.  Values of parameters to model (1.1).
    Group k Λ(k) p(k) μ(k) δ1(k) δ2(k) β1(k) β2(k) σ21(k) σ22(k) σ23(k)
    (Ⅰ) 1 0.60 0.65 0.35 0.20 0.25 0.25 0.35 0.200 0.450 0.100
    2 0.40 0.55 0.25 0.10 0.20 0.15 0.25 0.100 0.350 0.050
    (Ⅱ) 1 0.60 0.65 0.35 0.20 0.25 0.25 0.35 0.200 0.490 0.100
    2 0.40 0.55 0.25 0.10 0.20 0.15 0.25 0.100 0.400 0.050
    (Ⅲ) 1 0.60 0.30 0.25 0.20 0.25 0.55 0.65 0.002 0.003 0.003
    2 0.40 0.20 0.15 0.10 0.20 0.45 0.55 0.001 0.001 0.002

     | Show Table
    DownLoad: CSV
    Figure 1.  The extinction of the untreated drug addicts to model (1.1) under (Ⅰ) with initial conditions (S(0),U(0),T(0))=(0.70,0.50,0.40) and σ21(1)=0.2,σ21(2)=0.1,σ22(1)=0.45,σ22(2)=0.35,σ23(1)=0.1,σ23(2)=0.05.
    Figure 2.  The extinction of the untreated drug addicts to model (1.1) under (Ⅱ) with initial conditions (S(0),U(0),T(0))=(0.70,0.50,0.40) and σ21(1)=0.2,σ21(2)=0.1,σ22(1)=0.49,σ22(2)=0.4,σ23(1)=0.1,σ23(2)=0.05.

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

    Rs0=kSπkR0k=0.707>0.

    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.

    Figure 3.  The stationary distributions with same Markov chain (left) under MHOM (middle) and PTEMM (right) respectively.
    Figure 4.  Histogram of S(t),U(t),T(t) to model (1.1) with 50000 sample paths.

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

    Table 2.  Values of parameters to model (2.20).
    Group Λ p μ δ1 δ2 β1 β2 σ21 σ22 σ23
    (Ⅳ) 0.60 0.35 0.15 0.125 0.115 0.400 0.532 0.125 0.045 0.035
    (Ⅴ) 0.60 0.35 0.05 0.150 0.055 0.400 0.380 0.170 0.080 0.044
    (Ⅵ) 0.60 0.35 0.17 0.050 0.070 0.261 0.450 0.188 0.045 0.035

     | Show Table
    DownLoad: CSV

    By (Ⅳ), the conditions

    1+m1m2m1pβ1Λ=0.618>0,Λβ21m3β2p=0.043>0

    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.

    Figure 5.  Persistence and density function of model (2.20) around (1.078, 0.861, 0.407).

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

    β1Λm1pm1m2=0.105>0,Λβ21m3β2p=0.106>0

    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

    (β2m1β1m3)(m2+p)+β2(m1pβ1Λ)=0.160>0,Λβ21m3β2p=0.079>0
    Figure 6.  Persistence and density function of model (2.20) around (0.746, 1.670, 0.768).

    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.

    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.

    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).

    The authors declare there is no conflict of interest.



    [1] S. A. Ochoa-Orozco, J. C. Gutiérrez-Segura, A. M. Coral-Leiton, E. A. Trejos-Orozco, I. Gutirrez-Sanjun, J. D. Carvajal-Guevara, Chasing the dragon: A fatal case report of toxic leucoencé phalopathie due to inhaled heroin, Rev. Colomb. Psiquiat., 49 (2020), 289–292. https://doi.org/10.1016/j.rcp.2019.06.003 doi: 10.1016/j.rcp.2019.06.003
    [2] National, Institute on drug abuse. Available from: https://www.drugabuse.gov/drug-topics/opioids
    [3] E. White, C. Comiskey, Heroin epidemics, treatment and ODE modelling, Math. Biosci., 208 (2007), 312–324. https://doi.org/10.1016/j.mbs.2006.10.008 doi: 10.1016/j.mbs.2006.10.008
    [4] S. Djilali, T. M. Touaoula, M. S. El-Hadi, A heroin epidemic model: Very general nonlinear incidence treat-age and global stability, Acta. Appl. Math., 152 (2017), 171–194. https://doi.org/10.1007/s10440-017-0117-2 doi: 10.1007/s10440-017-0117-2
    [5] J. Wang, J. Wang, T. Kuniya, Analysis of an age-structured multi-group heroin epidemic model, Appl. Math. Comput., 347 (2019), 78–100. https://doi.org/10.1016/j.amc.2018.11.012 doi: 10.1016/j.amc.2018.11.012
    [6] X. Duan, X. Li, M. Martcheva, Qualitative analysis on a diffusive age-structured heroin transmission model, Nonlinear Anal.-Real World Appl., 54 (2020), 103105. https://doi.org/10.1016/j.nonrwa.2020.103105 doi: 10.1016/j.nonrwa.2020.103105
    [7] J. Liu, T. Zhang, Global behaviour of a heroin epidemic model with distributed delays, Appl. Math. Lett., 24 (2011), 1685–1692. https://doi.org/10.1016/j.aml.2011.04.019 doi: 10.1016/j.aml.2011.04.019
    [8] G. Huang, A. Liu, A note on global stability for a heroin epidemic model with distributed delay, Appl. Math. Lett., 26 (2013), 687–691. https://doi.org/10.1016/j.aml.2013.01.010 doi: 10.1016/j.aml.2013.01.010
    [9] X. Abdurahman, Z. Teng, L. Zhang, Global dynamics in a heroin epidemic model with difffferent conscious stages and two distributed delays, Int. J. Biomath., 12 (2019), 1950038. https://doi.org/10.1142/S1793524519500384 doi: 10.1142/S1793524519500384
    [10] M. Ma, S. Liu, J. Li, Bifurcation of a heroin model with nonlinear incidence rate, Nonlinear Dyn., 88 (2017), 555–565. https://doi.org/10.1007/s11071-016-3260-9 doi: 10.1007/s11071-016-3260-9
    [11] L. Chen, F. Wei, Study on a susceptible-exposed-infected-recovered model with nonlinear incidence rate, Adv. Differ. Equations, 2020 (2020), 206. https://doi.org/10.1186/s13662-020-02662-5 doi: 10.1186/s13662-020-02662-5
    [12] S. Djilali, S. Bentout, T. M. Touaoula, A. Tridane, S. Kumar, Global behavior of Heroin epidemic model with time distributed delay and nonlinear incidence function, Results Phys., 31 (2021), 104953. https://doi.org/10.1016/j.rinp.2021.104953 doi: 10.1016/j.rinp.2021.104953
    [13] S. Bentout, Y. Chen, S. Djilali, Global dynamics of an SEIR model with two age structures and a nonlinear incidence, Acta. Appl. Math., 171 (2021), 1–27. https://doi.org/10.1007/s10440-020-00369-z doi: 10.1007/s10440-020-00369-z
    [14] S. Liu, L. Zhang, Y. Xing, Dynamics of a stochastic heroin epidemic model, J. Comput. Appl. Math., 351 (2019), 260–269. https://doi.org/10.1016/j.cam.2018.11.005 doi: 10.1016/j.cam.2018.11.005
    [15] S. Liu, Z. Liang, X. Zhang, A. Li, Dynamics of a stochastic heroin epidemic modelwith bilinear incidence and varying population size, Int. J. Biomath., 12 (2019), 1950005. https://doi.org/10.1142/S1793524519500050 doi: 10.1142/S1793524519500050
    [16] Y. Wei, Q. Yang, G. Li, Dynamics of the stochastically perturbed heroin epidemic model under non-degenerate noises, Physica A, 526 (2019), 120914. https://doi.org/10.1016/j.physa.2019.04.150 doi: 10.1016/j.physa.2019.04.150
    [17] F. Wei, H. Jiang, Q. Zhu, Dynamical behaviors of a heroin population model with standard incidence rates between distinct patches, J. Frankl. Inst.-Eng. Appl. Math., 358 (2021), 4994–5013. https://doi.org/10.1016/j.jfranklin.2021.04.024 doi: 10.1016/j.jfranklin.2021.04.024
    [18] J. Liu, S. Wang, Dynamics in a stochastic heroin model with seasonal variation, Phys. A, 532 (2019), 121873. https://doi.org/10.1016/j.physa.2019.121873 doi: 10.1016/j.physa.2019.121873
    [19] M. Jovanović, V. Jovanović, Stability of stochastic heroin model with two distributed delays, Discrete Cont. Dyn. Sys.-B, 25 (2020), 2407–2432. https://doi.org/10.3934/dcdsb.2020016 doi: 10.3934/dcdsb.2020016
    [20] Q. Luo, X. Mao, Stochastic population dynamics under regime switching, J. Math. Anal. Appl., 334 (2007), 69–84. https://doi.org/10.1016/j.jmaa.2006.12.032 doi: 10.1016/j.jmaa.2006.12.032
    [21] M. Slatkin, The dynamics of a population in a Markovian environment, Ecology, 59 (1978), 249–256. https://doi.org/10.2307/1936370 doi: 10.2307/1936370
    [22] X. Zou, K. Wang, The protection zone for biological population in random environment, Math. Method. Appl. Sci., 36 (2013), 707–721. https://doi.org/10.1002/mma.2621 doi: 10.1002/mma.2621
    [23] S. He, F. Liu, Optimal finite-time passive controller design for uncertain nonlinear Markovian jumping systems, J. Frankl. Inst.-Eng. Appl. Math., 351 (2014), 3782–3796. https://doi.org/10.1016/j.jfranklin.2013.03.006 doi: 10.1016/j.jfranklin.2013.03.006
    [24] X. Zhang, D. Jiang, A. Alsaedi, T. Hayat, Stationary distribution of stochastic SIS epidemic model with vaccination under regime switching, Appl. Math. Lett., 59 (2016), 87–93. https://doi.org/10.1016/j.aml.2016.03.010 doi: 10.1016/j.aml.2016.03.010
    [25] D. Greenhalgh, Y. Liang, X. Mao, Modelling the effect of telegraph noise in the SIRS epidemic model using Markovian switching, Phys. A, 462 (2016), 684–704. https://doi.org/10.1016/j.physa.2016.06.125 doi: 10.1016/j.physa.2016.06.125
    [26] Q. Lin, L. Chen, C. Wen, F. Wei, Asymptotic properties of a stochastic Lotka-Volterra model with infinite delay and regime switching, Adv. Differ. Equations, 2018 (2018), 155. https://doi.org/10.1186/s13662-018-1609-8 doi: 10.1186/s13662-018-1609-8
    [27] Q. Liu, D. Jiang, N. Shi, Threshold behavior in a stochastic SIQR epidemic model with standard incidence and regime switching, Appl. Math. Comput., 316 (2018), 310–325. https://doi.org/27.10.1016/j.amc.2017.08.042 doi: 10.1016/j.amc.2017.08.042
    [28] H. Wang, D. Jiang, T. Hayat, A. Alsaedi, A. Bashir, Stationary distribution of stochastic NP ecological model under regime switching, Phys. A, 549 (2020), 124064. https://doi.org/10.1016/j.physa.2019.124064 doi: 10.1016/j.physa.2019.124064
    [29] N. D. Phu, D. O'Regan, T. D. Tuong, Longtime characterization for the general stochastic epidemic SIS model under regime-switching, Nonlinear Anal.-Hybrid Syst., 38 (2020), 100951. https://doi.org/10.1016/j.nahs.2020.100951 doi: 10.1016/j.nahs.2020.100951
    [30] X. Zhang, H. Peng, Stationary distribution of a stochastic cholera epidemic model with vaccination under regime switching, Appl. Math. Lett., 102 (2020), 106095. https://doi.org/10.1016/j.aml.2019.106095 doi: 10.1016/j.aml.2019.106095
    [31] B. Zhou, B. Han, D. Jiang, T. Hayat, A. Alsaedi, Ergodic stationary distribution and extinction of a hybrid stochastic SEQIHR epidemic model with media coverage, quarantine strategies and pre-existing immunity under discrete Markov switching, Appl. Math. Comput., 410 (2021), 126388. https://doi.org/10.1016/j.amc.2021.126388 doi: 10.1016/j.amc.2021.126388
    [32] J. Xu, Y. Wang, Z. Cao, Dynamics of a stochastic SIRS epidemic model with standard incidence under regime switching, Int. J. Biomath., 15 (2022), 2150074. https://doi.org/10.1142/S1793524521500741 doi: 10.1142/S1793524521500741
    [33] G. Li, Q. Yang, Y. Wei, Dynamics of stochastic heroin epidemic model with Levy jumps, J. Appl. Anal. Comput., 8 (2018), 998–1010. https://doi.org/10.11948/2018.99 doi: 10.11948/2018.99
    [34] F. Wei, C. Wang, Survival analysis of a single-species population model with fluctuations and migrations between patches, Appl. Math. Model., 81 (2020), 113–127. https://doi.org/10.1016/j.apm.2019.12.023 doi: 10.1016/j.apm.2019.12.023
    [35] F. Wei, L. Chen, Psychological effect on single-species population models in a polluted environment, Math. Biosci., 290 (2017), 22–30. https://doi.org/10.1016/j.mbs.2017.05.011 doi: 10.1016/j.mbs.2017.05.011
    [36] L. Chen, F. Wei, Persistence and distribution of a stochastic susceptible-infected-removed epidemic model with varying population size, Phys. A, 483 (2017), 386–397. https://doi.org/10.1016/j.physa.2017.04.114 doi: 10.1016/j.physa.2017.04.114
    [37] F. Wei, R. Xue, Stability and extinction of SEIR epidemic models with generalized nonlinear incidence, Math. Comput. Simul., 170 (2020), 1–15. https://doi.org/10.1016/j.matcom.2018.09.029 doi: 10.1016/j.matcom.2018.09.029
    [38] R. Lu, F. Wei, Persistence and extinction for an age-structured stochastic SVIR epidemic model with generalized nonlinear incidence rate, Phys. A, 513 (2019), 572–587. https://doi.org/10.1016/j.physa.2018.09.016 doi: 10.1016/j.physa.2018.09.016
    [39] F. Wei, J. Liu, Long-time behavior of a stochastic epidemic model with varying population size, Phys. A, 470 (2017), 146–153. https://doi.org/10.1016/j.physa.2016.11.031 doi: 10.1016/j.physa.2016.11.031
    [40] X. Mao, C. Yuan, Stochastic Differential Equations with Markovian Switching, Imperial College Press, London, 2006. https://doi.org/10.1142/p473
    [41] H. Peng, X. Zhang, Dynamics of a stochastic rabies epidemic model with Markovian switching, Int. J. Biomath., 14 (2021), 2150032. https://doi.org/10.1142/S1793524521500327 doi: 10.1142/S1793524521500327
    [42] B. Zhou, X. Zhang, D. Jiang, Dynamics and density function analysis of a stochastic SVI epidemic model with half saturated incidence rate, Chaos Solitons Fract., 137 (2020), 109865. https://doi.org/10.1016/j.chaos.2020.109865 doi: 10.1016/j.chaos.2020.109865
    [43] C. W. Gardiner, Handbook of stochastic methods for physics, chemistry and the natural sciences, Springer, Berlin, 1986. https://doi.org/10.1007/978-3-662-02452-2
    [44] H. Roozen, An asymptotic solution to a two-dimensional exit problem arising in population dynamics, SIAM J. Appl. Math., 49 (1989), 1793–1810. https://doi.org/10.2307/2101938 doi: 10.2307/2101938
    [45] D. Li, F. Wei, X. Mao, Stationary distribution and density function of a stochastic SVIR epidemic model, J. Frankl. Inst.-Eng. Appl. Math., 359 (2022), 9422–9449. https://doi.org/10.1016/j.jfranklin.2022.09.026 doi: 10.1016/j.jfranklin.2022.09.026
    [46] D. J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev., 43 (2001), 525–546. https://doi.org/10.1137/S0036144500378302 doi: 10.1137/S0036144500378302
    [47] X. Mao, F. Wei, T. Wiriyakraikul, Positivity preserving truncated Euler-Maruyama method for stochastic Lotka-Volterra competition model, J. Comput. Appl. Math., 394 (2021), 113566. https://doi.org/10.1016/j.cam.2021.113566 doi: 10.1016/j.cam.2021.113566
    [48] Q. Guo, W. Liu, X. Mao, R. Yue, The partially truncated Euler-Maruyama method and its stability and boundedness, Appl. Numer. Math., 115 (2017), 235–251. https://doi.org/10.1016/j.apnum.2017.01.010 doi: 10.1016/j.apnum.2017.01.010
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1343) PDF downloads(65) Cited by(0)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog