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

Optimal control strategies for a heroin epidemic model with age-dependent susceptibility and recovery-age

  • Correction on: AIMS Mathematics 6: 7318-7319
  • In the present manuscript, an age-structured heroin epidemic model is formulated with the assumption that susceptibility and recovery are age-dependent. Keeping in view some control measures for heroin addiction, a control problem for further analysis is presented. The main results are the existence of control variables, sensitivities, adjoint system and the setting of an optimal control problem. We used the techniques of weak derivatives and a general principle of Pontryagin's type for obtaining the optimal control problem. To compare our results, we demonstrated sample simulations which show the effect of control on the entire population.

    Citation: Asaf Khan, Gul Zaman, Roman Ullah, Nawazish Naveed. Optimal control strategies for a heroin epidemic model with age-dependent susceptibility and recovery-age[J]. AIMS Mathematics, 2021, 6(2): 1377-1394. doi: 10.3934/math.2021086

    Related Papers:

    [1] Asaf Khan, Gul Zaman, Roman Ullah, Nawazish Naveed . Correction: Optimal control strategies for a heroin epidemic model with age-dependent susceptibility and recovery-age. AIMS Mathematics, 2021, 6(7): 7318-7319. doi: 10.3934/math.2021429
    [2] Xin Yi, Rong Liu . An age-dependent hybrid system for optimal contraception control of vermin. AIMS Mathematics, 2025, 10(2): 2619-2633. doi: 10.3934/math.2025122
    [3] Houssine Zine, Abderrahim El Adraoui, Delfim F. M. Torres . Mathematical analysis, forecasting and optimal control of HIV/AIDS spatiotemporal transmission with a reaction diffusion SICA model. AIMS Mathematics, 2022, 7(9): 16519-16535. doi: 10.3934/math.2022904
    [4] Feliz Minhós, Ali Raza, Umar Shafique . An efficient computational analysis for stochastic fractional heroin model with artificial decay term. AIMS Mathematics, 2025, 10(3): 6102-6127. doi: 10.3934/math.2025278
    [5] Aqeel Ahmad, Cicik Alfiniyah, Ali Akgül, Aeshah A. Raezah . Analysis of COVID-19 outbreak in Democratic Republic of the Congo using fractional operators. AIMS Mathematics, 2023, 8(11): 25654-25687. doi: 10.3934/math.20231309
    [6] Huihui Liu, Yaping Wang, Linfei Nie . Dynamical analysis and optimal control of an multi-age-structured vector-borne disease model with multiple transmission pathways. AIMS Mathematics, 2024, 9(12): 36405-36443. doi: 10.3934/math.20241727
    [7] Shiwen Jing, Hairong Lian, Yiming Tang, Zhaohai Ma . Traveling waves for a nonlocal dispersal SIRS epidemic model with age structure. AIMS Mathematics, 2024, 9(4): 8001-8019. doi: 10.3934/math.2024389
    [8] Aziz Belmiloudi . Time-varying delays in electrophysiological wave propagation along cardiac tissue and minimax control problems associated with uncertain bidomain type models. AIMS Mathematics, 2019, 4(3): 928-983. doi: 10.3934/math.2019.3.928
    [9] Muhammad Altaf Khan, E. Bonyah, Yi-Xia Li, Taseer Muhammad, K. O. Okosun . Mathematical modeling and optimal control strategies of Buruli ulcer in possum mammals. AIMS Mathematics, 2021, 6(9): 9859-9881. doi: 10.3934/math.2021572
    [10] Folashade Agusto, Daniel Bond, Adira Cohen, Wandi Ding, Rachel Leander, Allis Royer . Optimal impulse control of West Nile virus. AIMS Mathematics, 2022, 7(10): 19597-19628. doi: 10.3934/math.20221075
  • In the present manuscript, an age-structured heroin epidemic model is formulated with the assumption that susceptibility and recovery are age-dependent. Keeping in view some control measures for heroin addiction, a control problem for further analysis is presented. The main results are the existence of control variables, sensitivities, adjoint system and the setting of an optimal control problem. We used the techniques of weak derivatives and a general principle of Pontryagin's type for obtaining the optimal control problem. To compare our results, we demonstrated sample simulations which show the effect of control on the entire population.



    Heroin is highly addictive and an illegal opioid drug and commonly known as hell dust. Usually, it is used in combination with other drugs or with alcohol which results from the risk of overdose. The overdose of heroin can cause a lot of troubles such as coma, shallow and slow breathing and death, etc [1]. It is observed that out of 10 heroin addicts, about 9 users used at least one other drug. Typically, people inject heroin, however, sometimes it may be snorted or smoked. Consequences of injecting heroin into a body include the risk of viral infection like HIV, HBV, HCV and bacterial infection of bloodstream and skin, etc. People whose age is between 18 and 25 or addicted to opioid painkiller, are more vulnerable to heroin addiction [2]. It is reported that the use of heroin is more than doubled among young adults and it tends to increase. Further, more than 45% of heroin users are addicted to prescription of opioid painkiller.

    Heroin addiction develops and propagates in a community like an epidemic disease and almost follows the same pattern. More precisely, people are susceptible to heroin addiction initially, and then they become active users of heroin and finally get recovery through some control measures or by self ceasing the use of heroin. The mathematical study of epidemic diseases has a rich history, and the readers are suggested to see the works of [3,4,5,6,7]. From the last two decades, the same modeling approaches have been extended to drug addiction problems for understanding its dynamics [8,9,10,11,12,13]. Heroin epidemic models comprises of ordinary differential equations and delay differential equations have been investigated and analyzed in [14,15,16,17,18,19,20,21] and references cited therein.

    In the context of partial differential equations (PDEs), mathematicians have developed some novel models called age-structured models to explore several insights into the heroin epidemic. Feng et. al. assumed the age-dependent susceptibility and qualitatively analyzed an age-structured heroin model [22]. They studied the existence of steady states and their local and global stabilities. Another paper by Feng et. al. considered treat-age while formulating the model and the authors presented the dynamical aspects of the model [12]. Subsequently, the theory was extended by Yang et al. [23] and Djilali et al. [13] by incorporating different types of incident rates. A recent paper of Wang et al. [24] dealt with the multi-group heroin epidemic that includes both treat and susceptibility age. Again the work was restricted only to the dynamical behavior of the concerned model.

    To reduce heroin addiction, government officials and health providers should take some agammaessive steps. These control measures include providing educational training to health care providers in prescribing opioid painkillers, expansion of medication-assisted treatments, making informed decisions, behavioral therapy, development and distribution of life-saving drugs, educating vulnerable people about heroin addiction and its consequences, etc [2]. It is also suggested to identify those societies who are at most risk to heroin addiction and provide them effective prevention strategies. For active heroin users, it is suggested to treat them with medications or behavioral therapy [1]. Also, heroin addiction takes place during the age of 18-25 among young people, thus, it is necessary to make aware vulnerable especially young people about the harm of heroin use. The key aim of modeling epidemics is to identify how to control such outbreaks and to eliminate the same from the community. Since disease epidemics and drugs have the same dynamics; therefore, it is necessary to apply the techniques of control theory to heroin models for controlling heroin addiction.

    The age-structured models cited earlier are indeed a great contribution to understanding the heroin epidemic and predicting its future behavior. However, it is worthy to notice that they have given less attention to controlling heroin addiction. Though the authors take into account treatment classes, it doesn't reflect the fact that what type of treatment strategy should be adapted which is more economical. In the current study, the authors aimed to include class-age in susceptible and recovered classes, whereas, heroin users are assumed to depend only on time. Further, to control heroin addiction in a community, we will use two control measures namely; treatment of active heroin users through medications and educating susceptible people about the risk of heroin use.

    The rest of the work is organized as follows. Section 1, is about the formulation and existence of the solution to the proposed without control model. The development of the control problem, the existence of optimal control variables, and the characterization of controls are carried out in Section 2. In Section 3, we presented a detailed comparison of both analytical and numerical results. Besides, we simulated both the control and without control problems and showed that the application of control measures could markedly reduce heroin addiction. Finally, concluding remarks are given at the end of the paper.

    Like infectious diseases, modeling the heroin epidemic follows the same transmission pattern [13,24]. While modeling the heroin epidemic, one must keep in mind some terms such as carriers of heroin use, relapse, time since heroin initiation to habitual use, age of susceptibility and recovery, treatment, and many more factors. It is also suggested to interview heroin addicts and note feedback from professionals in drug abuse-related fields. We assumed that the entire population is composed of individuals susceptible to heroin addiction, heroin users and recovered from drug-life which are respectively denoted by S(t),U(t) and R(t). The function s(t,a) denotes the density of susceptible population at any time t and susceptibility-age a. Thus the total susceptible population at any time t in terms of density function becomes

    S(t)=0s(t,a)da.

    In the same way, r(t,a) represent the density of recovered people at time t having recovery-age a and thus R(t) in terms of r(t,a) can be written as

    R(t)=0r(t,a)da.

    Further, we assumed that the addicts who have treated or they quit the use of heroin by self, again may start the use of heroin. These facts and many others motivate us to formulate the model of the form:

    s(t,a)t+s(t,a)a=β(a)U(t)s(t,a)μs(t,a), (2.1a)
    dU(t)dt=U(t)0β(a)s(t,a)da(μ+δ+p)U(t)+0γ(a)r(t,a)da, (2.1b)
    r(t,a)t+r(t,a)a=(μ+γ(a))r(t,a), (2.1c)
    s(t,0)=Λ, (2.1d)
    r(t,0)=pU(t), (2.1e)
    s(0,a)=s0(a),U(0)=U00,r(0,a)=r0(a). (2.1f)

    The parameters μ,δ,Λ and β(a) represent the natural death rate, heroin-induced removal rate, recruitment rate, and the age-specific transmission rate. The term p stand for self-cure (users who cease the use of heroin by self-decision) rate, whereas, γ(a) is the age-specific relapse rate. The terms s0(a) and r0(a) are non-negative functions of a and it represent the initial age distributions of susceptible and recovered population, respectively. Besides, the parameter used in the model are subject to the following hypothesis:

    H1. μ,p,δ,Λ>0.

    H2. β(a),γ(a)L+ with an essential upper bound ˉβ and ˉγ, respectively.

    H3. The function β(a) as well as γ(a) are Lipschitz continuous on R0, with Lipschitz coefficients Mβ and Mγ, respectively.

    For the sake of simplicity, let us consider some notations

    Q(t)=0β(a)s(t,a)da,L(t)=0γ(a)r(t,a)da,Γ(a)=ea0(μ+γ(τ))dτ,Ω(a)=ea0(β(η)U(ta+η)+μ)dη.

    Equations (2.1a) and (2.1c) can be solved by using the method of characteristics, thus

    s(t,a)={ΛΩ(a),t>a0;s0(at)Ω(a)Ω(at),at0, (2.2)

    and

    r(t,a)={pU(ta)Γ(a),t>a0;r0(at)Γ(a)Γ(at),at0. (2.3)

    Similarly, by applying direct integration to Eq (2.1b) and considering (2.2) and (2.3), we get solution of system (2.1) of the form

    S(t)=0s(t,a)da=t0ΛΩ(a)da+ts0(at)Ω(a)Ω(at)da, (2.4a)
    U(t)=t0[U(ta)Q(ta)+L(ta)]e(μ+δ+p)ada+U0e(μ+δ+p)t, (2.4b)
    R(t)=0r(t,a)da=t0pU(ta)Γ(a)da+tr0(at)Γ(a)Γ(at)da. (2.4c)

    Applying standard fixed point theorem of [25] to system (2.4), we can prove that there exist a nonnegative solution of system (2.1) over the interval [0,).

    The desired functional space for system (2.1) is of the form

    Z=L1+×R+×L1+,

    where L1+ denotes the space of positive and Lebesgue integrable function over [0,). The functions from the space L1+ are equipped with the norm

    (ϕ,x,ψ)Z=0|ϕ|da+|x|+0|ψ|da.

    Biologically, this norm is equivalent to total population. The initial conditions for system (2.1) can be written as

    Z0:=(s(0,),U0,r(0,))=(s0(),U0,r0())Z.

    Existence and uniqueness theorems along with hypothesis of the model suggest that system (2.1) has a unique nonnegative solution using conditions Z0Z. Further, a continuous semi-flow χ:R+×ZZ is defined by system (2.1) and solutions have a non-empty omega limit sets because of compact closure [26].

    Naturally, heroin is highly addictive and can cause a lot of damage if not addressed properly. To overcome its spread, agammaessive steps must be taken by health providers and government officials. It is observed that behavioral therapy and the administration of pharmacological medications are the most effective treatments for heroin addiction. Although, the use of medicines for treatment purpose have some side effects (such as nausea, bone pain, and vomiting, etc) that can be treated with additional supplements. Further, it is noted that the treatment of heroin addicts through medicines has a positive impact and as a result, the rate of crimes and heroin consumption tends to decrease. The present study aims to include two types of measures as control variables; education campaign and pharmacological medications. The former control variable somehow similar to behavioral therapy and will be implemented upon the susceptible population. The educational campaign is widely used in giving-up smoking models as a control variable and it has shown a great effect. The later control will be imposed on heroin addicts. The strategy will be to minimize the susceptible and heroin user population and to maximize the recovered people with the least cost.

    From the dynamical perspective, optimal control theory is the most powerful and useful tool than the calculus of variation, particularly in the epidemic problems. Since, heroin addiction is closely related to mathematical epidemiology, therefore, it is proposed to implement the techniques of control theory to minimize the number of susceptible and heroin users. As the problem under discussion lying in infinite-dimensional spaces, thus, we will use the maximum principle type results introduced by Li and Yong [27] and suggested by [28].

    The proposed control problem consist of three state variables S(t),U(t),R(t) and two control variables v1(t),v2(t) from the admissible control set

    Vad={v1,v2L2(0,T)|0v1(t)l110v2(t)l21,t[0,T]}. (3.1)

    Physical interpretation of the control v1(t) is educating susceptible population regarding heroin addiction and its consequences, whereas, v2(t) means to treat active heroin users. Since, v1(t)s(t,a) denotes the density of susceptible who have been educated. Therefore, the term v1(t)s(t,a) needs to subtracted from the equation which represent the dynamics of the susceptible population and the same will be added to the equation of recovered people. Likewise, the term v2(t)U(t) represent the number of heroin addicts who were treated through medicines, thus, it should be subtracted from the equation which describe the dynamics of heroin users. Since, the equation for recovered population is a PDE, therefore, the term v2(t)U(t) will be included to the boundary condition r(t,0) which will make it a type of boundary control problem. Further, nature of the problem demands that vi(t)0 for i=1,2. The case vi(t)=0 means implementation of no control measures, whereas, vi(t)=1 shows that the entire susceptible are educated and the heroin users are treated. Due to wide used of quadratic term for the control variable, our objective functional will also depends quadratically upon the controls. Thus, the objective functional will takes the following form

    J(v1,v2)=T0[S(t)+U(t)+B12v21(t)+B22v22(t)]dt, (3.2)

    subject to the model

    s(t,a)t+s(t,a)a=β(a)U(t)s(t,a)μs(t,a)v1(t)s(t,a), (3.3a)
    dU(t)dt=U(t)0β(a)s(t,a)da(μ+δ+p+v2(t))U(t)+0γ(a)r(t,a)da, (3.3b)
    r(t,a)t+r(t,a)a=(μ+γ(a))r(t,a)+v1(t)s(t,a), (3.3c)
    s(t,0)=Λ, (3.3d)
    r(t,0)=(p+v2(t))U(t), (3.3e)
    s(0,a)=s0(a),U(0)=U00,r(0,a)=r0(a). (3.3f)

    The terms B1,B2 are positive constants and it is called balancing factors. The constant B1 and B2 represent the level of acceptance of education campaign and medication by a susceptible and an addict, respectively. The theorem stated and proved below is about the existence of optimal control variables.

    Theorem 3.1. There exists a pair of optimal control variables v1,v2Vad such that, J(v1,v2)=minv1,v2VadJ(v1,v2) subject to the system (3.3).

    proof. Clearly, the states and the control variables are nonnegative and bounded. The set Vad of all the controls is convex and closed. Further, the boundedness of the control variables suggests the compactness of the control system that is necessary for the existence of an optimal control. Also in the objective functional, the integrand is a convex function on the control set Vad. In addition, there exists two positive real numbers η1 and η2 and a constant ξ>1 such that

    J(v1,v2)η1+η2(|v1|2+|v2|2)ξ2.

    We noticed that our control set and the objective functional fulfill all the hypotheses of Theorem (13.8.1) in [29] that guarantee the existence of the optimal system. Thus, there exists an optimal control problem and hence the proof.

    In order to obtain the necessary criteria of optimality, we will use the Gateaux derivative rule (also called the weak derivatives). For this purpose, we will introduce other controls vϵi(t)=vi(t)+ϵhi(t),i=1,2, where hi(t)0 are variation functions and 0<ϵ<1. As the state variables are continuously dependent on the control variables, thus, we can write S(t)=S(v1(t),v2(t)),U(t)=U(v1(t),v2(t)) and R(t)=R(v1(t),v2(t)). By perturbing the controls variables will disturb the state as well and hence state variable will takes the form Sϵ(t)=S(vϵ1,vϵ2),Uϵ(t)=U(vϵ1,vϵ2) and Rϵ(t)=R(vϵ1,vϵ2), where

    Sϵ(t)=0sϵ(t,a)da,andRϵ(t)=0rϵ(t,a)da.

    Accordingly due to vϵ1(t),vϵ2(t), system (3.3) becomes

    sϵ(t,a)t+sϵ(t,a)a=β(a)Uϵ(t)sϵ(t,a)μsϵ(t,a)vϵ1(t)sϵ(t,a), (3.4a)
    dUϵ(t)dt=(0β(a)sϵ(t,a)da(μ+δ+p+vϵ2(t)))Uϵ(t)+0γ(a)rϵ(t,a)da, (3.4b)
    rϵ(t,a)t+rϵ(t,a)a=(μ+γ(a))rϵ(t,a)+vϵ1(t)sϵ(t,a), (3.4c)
    sϵ(t,0)=Λ, (3.4d)
    rϵ(t,0)=(p+vϵ2(t))Uϵ(t), (3.4e)
    sϵ(0,a)=s0(a),Uϵ(0)=U0,rϵ(0,a)=r0(a). (3.4f)

    To form the weak derivatives, first we will subtract system (3.3) from (3.4) and will make the difference quotients sϵ(t,a)s(t,a)ϵ, Uϵ(t)U(t)ϵ and rϵ(t,a)r(t,a)ϵ. Next, we will find the sensitivities ˉs(t,a),ˉU(t) and ˉr(t,a) which can be obtain by applying limϵ0 to weak derivatives. The following linearized version of system (3.3) must be satisfied by the sensitivities

    ˉs(t,a)t+ˉs(t,a)a=β(a)[U(t)ˉs(t,a)+ˉU(t)s(t,a)]μˉs(t,a)v1(t)ˉs(t,a)h1(t)s(t,a), (3.5a)
    dˉU(t)dt=ˉU(t)0β(a)s(t,a)da+U(t)0β(a)ˉs(t,a)da(μ+δ+p)ˉU(t)v2(t)ˉU(t)h2(t)U(t)+0γ(a)ˉr(t,a)da, (3.5b)
    ˉr(t,a)t+ˉr(t,a)a=(μ+γ(a))ˉr(t,a)+v1(t)ˉs(t,a)+h1(t)s(t,a), (3.5c)
    ˉs(t,0)=0, (3.5d)
    ˉr(t,0)=pˉU(t)+v2(t)ˉU(t)+h2(t)U(t), (3.5e)
    ˉs(0,a)=0,ˉU(0)=0,ˉr(0,a)=0. (3.5f)

    Next, we will use the sensitivities to obtain the adjoint system. In the case of ODEs, the adjoint are chosen such that the sensitivity terms drop out. However in PDE models, we need the sensitivity equations for finding the adjoint system. Each equation in the sensitivity PDEs can be written in a simplified form as

    LiˉΨi=Fi(s,U,r,v1,v2),

    with necessary initial and boundary conditions and ˉΨi's are the corresponding variables ˉs,ˉU and ˉr for i=1,2,3 respectively. Here, the linear operators Li comes from the linearization of the corresponding state PDE operators. To obtain a characterization of the control variables, it is necessary to derive the adjoint system. The operators in the adjoint system are the adjoint operators of the Li acting on ˉs,ˉU and ˉr in the sensitivity PDEs. The operators Li and the adjoint operator Li are related by

    λi,LiˉΨi=Liλi,ˉΨi,

    where , is the L2-inner product. By using the integration by parts in multidimensions, we can throw the derivatives on each state from the operators Li onto the derivatives of the adjoint variables λi in the operators Li. Once we find the adjoint operators, then one can easily obtain the adjoint equations by using

    Liλi=(integrand ofJ)x,

    where x stand for the state variables. To do so, first we will consider (3.5a) as

    0=ˉs(t,a)t+ˉs(t,a)a+[β(a)U(t)+μ+v1(t)]ˉs(t,a)+[β(a)ˉU(t)+h1(t)]s(t,a),λ1(t,a),=ˉs(t,a),tλ1(t,a)aλ1(t,a)+β(a)U(t)λ1(t,a)+μλ1(t,a)+v1(t)λ1(t,a)+ˉU(t),β(a)s(t,a)λ1(t,a)+T00h1(t)s(t,a)λ1(t,a)da, (3.6)

    with subsidiary

    λ1(t,)=0,λ1(T,a)=0. (3.7)

    Here the notation , is defined as f1,f2=T00f1f2dadt. On the other hand, Equation (3.5d) will take the form

    0=dˉU(t)dtˉU(t)0β(a)s(t,a)daU(t)0β(a)ˉs(t,a)da+(μ+δ+p)ˉU(t)+v2(t)ˉU(t)+h2(t)U(t)0γ(a)ˉr(t,a)da,λ2(t),=ˉU(t),dλ2(t)dtλ2(t)0β(a)s(t,a)da+(μ+δ+p)λ2(t)+v2(t)λ2(t)ˉs(t,a),U(t)β(a)λ2(t)ˉr(t,a),γ(a)λ2(t)+T0λ2(t)h2(t)U(t)dt, (3.8)

    with

    λ2(T)=0, (3.9)

    where f1,f2=T0f1f2dt. In the same way, Equation (3.5c) can be written in the form of

    0=ˉr(t,a)t+ˉr(t,a)a+(μ+γ(a))ˉr(t,a)v1(t)ˉs(t,a)h1(t)s(t,a),λ3(t,a),0=ˉr(t,a),λ3(t,a)tλ3(t,a)a+(μ+γ(a))λ3(t,a)ˉs(t,a),v1(t)λ3(t,a)ˉU(t),[p+v2(t)]λ3(t,0)T00h1(t)s(t,a)λ3(t,a)dadtT0h2(t)U(t)λ3(t,0)dt, (3.10)

    along the condition

    λ3(t,)=0,λ3(T,a)=0. (3.11)

    For future use, we will compute the weak derivative of the objective functional J(v1,v2) as

    0J(v1,v2)=T0(ˉS(t)+ˉU(t)+B1h1(t)v1(t)+B2h2(t)v2(t))dt. (3.12)

    By rewriting and suitable arrangements, Equations (3.6), (3.8) and (3.10) yields the following adjoint system

    λ1(t,a)t+λ1(t,a)a=β(a)U(t)λ1(t,a)+[μ+v1(t)]λ1(t)β(a)U(t)λ2(t)v1(t)λ3(t,a)1, (3.13a)
    dλ2(t)dt=λ2(t)0β(a)s(t,a)da+(μ+δ+p)λ2(t)+v2(t)λ2(t)[p+v2(t)]λ3(t,0)+0β(a)s(t,a)λ1(t,a)da1, (3.13b)
    λ3(t,a)t+λ3(t,a)a=[μ+γ(a)]λ3(t,a)γ(a)λ2(t). (3.13c)

    with transversality conditions

    λ1(T,a)=0,λ2(T)=0,λ3(T,a)=0, (3.14)

    and boundary conditions

    λ1(t,)=0,λ3(t,)=0. (3.15)

    The following theorem will conclude our main discussion regarding the existence of an optimal control problem.

    Theorem 3.2. If v1,v2 are in Vad is a pair of optimal control variables which minimizing J(v1,v2) and (S(t),U(t),R(t)) and (λ1(t,a),λ2(t),λ3(t,a)) are the corresponding state and adjoint variables, respectively, then

    v1(t)=min{l1,max{0,0(λ1(t,a)λ3(t,a))s(t,a)daB1}}, (3.16)

    and

    v2(t)=min{l2,max{0,(λ2(t)λ3(t,0))U(t)B2}}. (3.17)

    proof. With the help of system (3.13), we can write the weak derivative of the objective functional defined in (3.12) as

    0J(v)=T0[ˉS(t)+ˉU(t)+B1h1(t)v1(t)+B2h2(t)v2(t)]dt,=T0ˉS(t)[λ1(t,a)tλ1(t,a)a+(β(a)U(t)+μ+v1(t))λ1(t,a)β(a)U(t)λ2(t)v1(t)λ3(t,a)]dt+T0ˉU(t)[dλ2(t)dt+(μ+δ+p0β(a)s(t,a)da)λ2(t)+v2(t)λ2(t)(p+v2(t))λ3(t,0)+0β(a)s(t,a)λ1(t,a)da]dt+T0ˉR(t)[λ3(t,a)tλ3(t,a)a+(μ+γ(a))λ3(t,a)γ(a)λ2(t)]dt+T0(B1h1(t)v1(t)+B2h2(t)v2(t))dt, (3.18)

    where

    ˉS(t)=0ˉs(t,a)da,andˉR(t)=0ˉr(t,a)da.

    Keeping in mind Eqs (3.6), (3.8) and (3.10), we can write the inequality (3.18) as

    0T0[0(λ3(t,a)λ1(t,a))s(t,a)da+B1v1(t)]h1(t)dt+T0[(λ3(t,0)λ2(t))U(t)+B2v2(t)]h2(t)dt, (3.19)

    for all v1(t),v2(t)Vad. As the variation functions h1(t),h2(t)0 on Vad, thus the rest of the integrand must be equal to zero. Hence

    v1(t)=0(λ1(t,a)λ3(t,a))s(t,a)daB1,

    and

    v2(t)=(λ2(t)λ3(t,0))U(t)B2.

    By considering the lower and upper bounds of the control variables as well as the optimal values of the state variables, we have

    v1(t)=min{l1,max{0,0(λ1(t,a)λ3(t,a))s(t,a)daB1}}, (3.20)

    and

    v2(t)=min{l2,max{0,(λ2(t)λ3(t,0))U(t)B2}}. (3.21)

    Formulas developed in relations (3.20) and (3.21) are the required characterization of the optimal controls.

    By using these optimal value of the control variables, the desired optimal system can be written as

    s(t,a)t+s(t,a)a=β(a)U(t)s(t,a)μs(t,a)v1(t)s(t,a), (3.22a)
    dU(t)dt=U(t)0β(a)s(t,a)da(μ+δ+p+v2(t))U(t)+0γ(a)r(t,a)da, (3.22b)
    r(t,a)t+r(t,a)a=(μ+γ(a))r(t,a)+v1(t)s(t,a), (3.22c)
    s(t,0)=Λ, (3.22d)
    r(t,0)=(p+v2(t))U(t), (3.22e)
    s(0,a)=s0(a),U(0)=U00,r(0,a)=r0(a). (3.22f)

    In the previous section, we have presented analytically an optimal control problem keeping in view some predetermined goals. In the present section, we will try to simulate both the optimal control and without control systems and to compare the obtained results. The numerical method to be used will be the finite difference method as the problem containing both ordinary and partial differential equations. This method approximates the derivative by difference quotients and it approximates the solution at each point of the domain. Let the maximum class-age is A and the maximum time is T both of which are finite. The age and time directions are divided into equal space with the same step size Δa=Δt due to together progression. Thus, all of the points in time and age directions can be computed as ak=kΔt and tn=nΔt, respectively. Assume the M and N are the number of points to be incorporated in age time directions so that A=MΔt and T=NΔt. The numerical method computes approximations to the solutions at the mesh points:

    snks(ak,tn),UnU(tn),rnkr(ak,tn).

    Adopting the similar procedure presented in [30], system (3.22) will gives the following simulation scheme

    sn+1k+1=snk1+(βkUn+1+μ+vm+11)Δt,Un+1=Un+ΔtMk=1γkrn+1kΔt1+(μ+δ+p+vn+12Mk=1βksn+1kΔt)Δt,rn+1k+1=rnk+vn+11sn+1k1+(μ+γk)Δt,un+10=Sn+1Ki=1βiun+1iΔt,sn+10=Λ,rn+10=(p+vn+12)Un+1,s0k=ϕ1k,U0=U0,r0k=ψ1k.

    Using the same procedure, the schemes for adjoint system (3.13) and system without control (2.1) are developed and then coded in MATLAB along with characterization of control variables defined in (3.20) and (3.21). For simulation purpose, we assumed that Δt=0.05, maximum class-age a=40 years and maximum time T for the control strategy is 40 days. Furthermore, the initial data related to all classes were chosen according to assumption of the model. The values and justification of parameter used in simulation is presented in Table 1. Besides, at initial stage of the control strategy we assumed v1(t)=v2(t)=0.

    Table 1.  Parameters baseline values used in numerical simulation.
    Parameters Values References
    Λ 0.75 Assumed using [19]
    μ 0.085 [19]
    δ 0.00014 Calculated using [19]
    p 0.0081 Assumed
    l1 1 [31]
    l2 1 [31]
    B1 100 Assumed with the help of [31]
    B2 150 Assumed
    β(a) 0.05e0.1a Assumed with the help of [23]
    γ(a) e0.1a Assumed

     | Show Table
    DownLoad: CSV

    Figures 1 represent the dynamics of all classes with respect to class-age and time. Sub-figure 1a shows the density of susceptible population without control, whereas, Figure 1b shows the same population density with control. Comparison reveal that there are more susceptible to heroin addiction when the control was not implemented. As the control strategy was launched, the control shows his effect and susceptible people tend to decrease. As a result, fewer people will tend to initiate the use of heroin. Figure 1c is self-explanatory and depict that both the controls are effective for reducing the number of heroin addicts in a community. Figures 1d and 1e show the recovered population and again it is clear from the graphs that the recovered population tends to increase as we impose the control strategy.

    Figure 1.  The plot shows the effectiveness of control strategy in all classes.

    To understand the effectiveness of education control, in Figure 2 we have chosen some sample curves (solution profiles) from the density function of the susceptible population at fixed time and age. Figure 2a represents the solution profile of s(t,a) at fixed age a=15. The figure shows that although education campaign has reverse effect on young people at initial stages but with the passage of time it will work. After 20 days of strategy, the population of young susceptible tends to decline and finally there will be no susceptible whose age of susceptibility is between 15-22. Likewise, Figure 2b represent the solution profile of those susceptible whose age of susceptibility is 30. Again this figure suggest that the education campaign is much more effective in the long run. Similarly, figures 2c and 2d represent solution profiles of s(t,a) at fixed time t=15 and 20, respectively. These figures also indicate the effectiveness of awareness among susceptible. Again, these figures suggest that such type of campaign should be launched for the long run to achieve the desired outcomes.

    Figure 2.  The plot shows solution profiles of susceptible population at fixed time and class-age both with and without control.

    In the same way, we have plotted sample solution curves from r(t,a) at fixed age and time in Figures 3a-3d. These figures shows a much clear effect of both awareness and treatment. Further, it shows that these control measures should be implemented in pairs; the former will control the susceptible people and the later will recover the active heroin users. Finally, Figures 4a and 4b represent the dynamics of the control variables v1(t) and v2(t), respectively.

    Figure 3.  The plot shows solution profiles of recovered population at fixed time and class-age both with and without control.
    Figure 4.  The figure represent the dynamics of control variables.

    In this manuscript, we have developed an age-structured heroin epidemic model with continuous age-structure in susceptible and recovered classes using previously formulated models on the heroin epidemic. Using two control variables and a suitable objective functional, we investigated a control problem. Using the generalized maximum principle of Pontryagin's type and the Gateaux derivative rule, the authors presented an optimal control problem and characterized the optimal values of the control variables. To compare the theoretical results, we presented some simulation. Simulation of the model suggests that to control the use of heroin in a community, one must impose these two control measures in combination.

    The authors are grateful to the anonymous reviewers for their valuable comments.

    All authors carried out the proof. All authors conceived of the study, and participated in its design and coordination.

    All authors declare that they have no competing interests.



    [1] NIDA, Heroin, National Institute on Drug Abuse, 27 June 2019. Available from: https: //www.drugabuse.gov/publications/drugfacts/heroin.
    [2] Today's Heroin Epidemic, Center for Disease Control and Prevention (CDC), July 2015.
    [3] M. Iannelli, Mathematical Theory of Age-Structured Population Dynamics, Giardini Editori e Stampatori, Pisa, 1994.
    [4] H. R. Thieme, Mathematics in Populations Biology, Princeton University Press, Princeton, 2003.
    [5] H. R. Thieme, C. Castillo-Chavez, How many infection-age-dependent infectivity affect the dynamics of HIV/AIDS?, SIAM J. Appl. Math., 53 (2014), 14-47.
    [6] S. P. Rajasekar, M. Pitchaimani, Quanxin Zhu, Dynamic threshold probe of stochastic SIR model with saturated incidence rate and saturated treatment function, Physica A: Statistical Mechanics and its Applications, 535 (2019), 122300.
    [7] G. Zamana, Y. Saitob, M. Khanc, Optimal vaccination of an endemic model with variable infectivity and infinite delay, Zeitschrift für Naturforschung A, 68 (2013), 677-685.
    [8] M. Ma, S. Liu, J. Li, Bifurcation of a heroin model with nonlinear incidence rate, Nonlinear Dyn., 88 (2017), 555-565. doi: 10.1007/s11071-016-3260-9
    [9] G. P. Samanta, Dynamic behaviour for a nonautonomous heroin epidemic model with time delay, J. Appl. Math. Comput., 35 (2011), 161-178. doi: 10.1007/s12190-009-0349-z
    [10] I. M. Wangari, L. Stone, Analysis of a heroin epidemic model with saturated treatment function, J. Appl. Math., 2017 (2017), 1-22.
    [11] J. Liu, T. Zhang, Global behaviour of a heroin epidemic model with distributed delays, Appl. Math. Lett., 24 (2011), 1685-1692. doi: 10.1016/j.aml.2011.04.019
    [12] B. Fang, X. Li, M. Martcheva, L. Cai, Global asymptotic properties of a heroin epidemic model with treat-age, Appl. Math. Comput., 263 (2015), 315-331.
    [13] S. Djilali, T. M. Touaoula, S. E. H. Miri, A heroin epidemic model: very general non linear incidence, treat-age, and global stability, Acta Appl. Math., 152 (2017), 171-194.
    [14] E. White, C. Comiskey, Heroin epidemics, treatment and ODE modelling, Math. Biosci., 208 (2007), 312-324.
    [15] G. Mulone, B. Straughan, A note on heroin epidemics, Math. Biosci., 218 (2009), 138-141.
    [16] J. Mushanyu, F. Nyabadza, G. Muchatibaya, A. G. R. Stewart, Modelling multiple relapses in drug epidemics, Ricerche Mat., 65 (2016), 37-63. doi: 10.1007/s11587-015-0241-0
    [17] B. Fang, X. Li, M. Martcheva, L. Cai, Global stability for a heroin model with two distributed delays, Disc. Cont. Dyna. Sys., 19 (2014), 715-733.
    [18] Y. Muroya, H. Li, T. Kuniya, Complete global analysis of an SIRS epidemic model with graded cure and incomplete recovery rates, J. Math. Anal. Appl., 410 (2014), 719-732. doi: 10.1016/j.jmaa.2013.08.024
    [19] F. Nyabadza, S. D. Hove-Musekwa, From heroin epidemics to methamphetamine epidemics: Modelling substance abuse in a South African province, Math. Biosci., 225 (2010), 132-140.
    [20] G. Huang, A. Liu, A note on global stability for a heroin epidemic model with distributed delay, Appl. Math. Lett., 26 (2013), 687-691. doi: 10.1016/j.aml.2013.01.010
    [21] X. Liu, J. Wang, Epidemic dynamics on a delayed multi-group heroin epidemic model with nonlinear incidence rate, J. Nonlinear Sci. Appl., 9 (2016), 2149-2160. doi: 10.22436/jnsa.009.05.20
    [22] B. Fang, X. Li, M. Martcheva, L. Cai, Global stability for a heroin model with age-dependent susceptibility, J. Syst. Sci. Complex., 28 (2015), 1243-1257. doi: 10.1007/s11424-015-3243-9
    [23] J. Yang, X. Li, F. Zhang, Global dynamics of a heroin epidemic model with age structure and nonlinear incidence, Int. J. Biomath., 9 (2016), 1650033. doi: 10.1142/S1793524516500339
    [24] J. Wang, J. Wang, T. Kuniya, Analysis of an age-structured multi-group heroin epidemic model, Appl. Math. Comput., 347 (2019), 78-100.
    [25] G. Grippenberg, S. O. Londen, O. Staffans, Volterra Integral and Functional Equations, Cambridge University Press, Cambridge, 1990.
    [26] G. F. Webb, Theory of Nonlinear Age-dependent Population Dynamics, Marcel Dekker: New York, 1985.
    [27] X. Li, J. Yong, Optimal Control Theory for Infinite Dimensional Systems, Birkhäuser: Boston, 1995.
    [28] A. Khan, G. Zaman, Optimal control strategy of SEIR endemic model with continuous agestructure in the exposed and infectious classes, Optim. Contr. Appl. Meth., 39 (2018), 1716-1727. doi: 10.1002/oca.2437
    [29] D. L. Lukes, Differential Equations: Classical to Controlled, Mathematics in Science and Engineering, Academic Press: New York, 1982.
    [30] M. Martcheva, Age-Structured Epidemic Models. In: An Introduction to Mathematical Epidemiology, Texts in Applied Mathematics, Vol 61, Springer, Boston, MA, 2015.
    [31] M. Thater, K. Chudej, H. J. Pesch, Optimal vaccination strategies for an SEIR model of infectious diseases with logistic growth, Math. Biosc. Eng., 15 (2017), 485-505. doi: 10.3934/mbe.2018022
  • This article has been cited by:

    1. Soufiane Bentout, Sunil Kumar, Salih Djilali, Hopf bifurcation analysis in an age-structured heroin model, 2021, 136, 2190-5444, 10.1140/epjp/s13360-021-01167-8
    2. Asaf Khan, Gul Zaman, Roman Ullah, Nawazish Naveed, Correction: Optimal control strategies for a heroin epidemic model with age-dependent susceptibility and recovery-age, 2021, 6, 2473-6988, 7318, 10.3934/math.2021429
    3. Ali Raza, Yu-Ming Chu, Mohd Yazid Bajuri, Ali Ahmadian, Nauman Ahmed, Muhammad Rafiq, Soheil Salahshour, Dynamical and nonstandard computational analysis of heroin epidemic model, 2022, 34, 22113797, 105245, 10.1016/j.rinp.2022.105245
    4. Eminugroho Ratna Sari, Fajar Adi-Kusumo, Lina Aryati, Mathematical analysis of a SIPC age-structured model of cervical cancer, 2022, 19, 1551-0018, 6013, 10.3934/mbe.2022281
    5. Joan C. Micó, A Population Pyramid Dynamics Model and Its Analytical Solution. Application Case for Spain, 2022, 10, 2227-7390, 3443, 10.3390/math10193443
    6. P. T. Sowndarrajan, L. Shangerganesh, A. Debbouche, D. F. M. Torres, Optimal control of a heroin epidemic mathematical model, 2022, 71, 0233-1934, 3107, 10.1080/02331934.2021.2009823
    7. A. Abidemi, J. O. Akanni, Dynamics of illicit drug use and banditry population with optimal control strategies and cost-effectiveness analysis, 2022, 41, 2238-3603, 10.1007/s40314-022-01760-2
    8. SARA SALEM ALZAID, BADR SAAD T. ALKAHTANI, GLOBAL ANALYSIS OF DIFFERENT COMPARTMENTS IN A GIVING-UP SMOKING MODEL, 2022, 30, 0218-348X, 10.1142/S0218348X22401284
    9. Chelsea Spence, Mary E. Kurz, Thomas C. Sharkey, Bryan Lee Miller, Scoping Literature Review of Disease Modeling of the Opioid Crisis, 2024, 0279-1072, 1, 10.1080/02791072.2024.2367617
    10. Xinxin Wang, Xiaoyun Wang, Fengqin Zhang, DYNAMIC ANALYSIS OF A DRUG TRANSMISSION MODEL WITH ANTI-DRUG EDUCATION AND MEDIA COVERAGE, 2023, 13, 2156-907X, 2184, 10.11948/20220430
    11. Lucas Böttcher, Tom Chou, Maria R D’Orsogna, Sandro Galea, Forecasting drug-overdose mortality by age in the United States at the national and county levels, 2024, 3, 2752-6542, 10.1093/pnasnexus/pgae050
    12. Salih Djilali, Yuming Chen, Soufiane Bentout, Dynamics of a delayed nonlocal reaction–diffusion heroin epidemic model in a heterogenous environment, 2024, 0170-4214, 10.1002/mma.10327
    13. Eric Numfor, Necibe Tuncer, Maia Martcheva, Optimal control of a multi-scale HIV-opioid model, 2024, 18, 1751-3758, 10.1080/17513758.2024.2317245
    14. Miguel Vivas-Cortez, Abu Bakar, M.S. Alqarni, Nauman Raza, Talat Nazir, Muhammad Farman, Global stability and modeling with a non-singular kernel for fractional order heroin epidemic model: Insights from different population studies, 2024, 36, 10183647, 103329, 10.1016/j.jksus.2024.103329
    15. Lucas Böttcher, Tom Chou, Maria R. D’Orsogna, Modeling and forecasting age-specific drug overdose mortality in the United States, 2023, 232, 1951-6355, 1743, 10.1140/epjs/s11734-023-00801-z
    16. Feliz Minhós, Ali Raza, Umar Shafique, An efficient computational analysis for stochastic fractional heroin model with artificial decay term, 2025, 10, 2473-6988, 6102, 10.3934/math.2025278
  • Reader Comments
  • © 2021 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(4234) PDF downloads(220) Cited by(16)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog