# Stability Analysis and Optimal Control Strategies of Giving Up Relapse Smoking Model with Bilinear and Harmonic Mean Type of Incidence Rates

Authors: Badr Saad T. Alkahtani
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

In this manuscript, a giving-up smoking model is develop using bilinear incidence rate, harmonic mean type of incidence rate and keeping in view the relapse factor associated to smoking. It is shown that the model' solution is bounded and positive for the appropriate initial data. The equilibria of the model are obtained, and it is proved that the smoking-free equilibrium is both locally and globally asymptotically stable forR0 less than unity. It is shown that the model has a positive light smoker present equilibrium whenever R0 is greater than one, and it is locally asymptotically stable if we have 1<R0<1+2β1μ+d+δ/Λβ2. Conditions for the global stability of light smoker present equilibrium are rigorously investigated. Also, it is proved that the model has a smoking present equilibrium when R0 satisfies some condition which is investigated both for local and global behavior. By considering a few control measures, optimal control strategies are achieved with the help of Pontryagin’s maximum principle. The analytical results are verified numerically, and effectiveness of the control program is presented.


---

## Body

## 1. Introduction

Among other diseases, infectious diseases are the most alarming threat to humanity from the last few centuries. The plague of Athens is considered to be the first ever epidemic which affects human life to a great extent [1, 2]. The life in Egypt and Roman was completely demolished by the smallpox in 165-180 C.E. in which millions of people died [3]. The epidemic Black Death in Europe is observed to be the first well-documented epidemic which killed more than 50 millions of people that region and Mediterranean. Various other epidemics affected human life throughout the history, and the current corona virus disease is the last epidemic whose disasters are in front of us. Thus, it is very crucial to understand the dynamics of such infectious disease (particularly the emerging epidemics) and to control its spread in early transmission phase. In the context of these extensive illnesses, the tobacco pandemic is one of the world’s most serious health threats. According to statistics, up to 50% of smokers die as a result of their habit, a tobacco-related mortality occurs every 8 seconds, and 10 percent of the adult population dies as a result of tobacco-related illnesses [4].Smoking habit grows and spreads across a society in the same way as an epidemic illness does, and generally, it nearly follows the same process. To be specific, people are prone to smoking at first, then become active users, and eventually recover through certain control methods or self-abandonment of its usage. In the early twentieth century, William Hamer’s work on mathematical modeling of infectious illnesses achieved substantial progress. Perhaps, Hamer was the first who introduces mass action law in epidemic modeling. However, Sir Ronald Ross is considered the father of today’s mathematical biology due to his work on malaria. In 1911, he published a book regarding malaria in which he used mathematical models for describing the dynamics of the infection and calculated the threshold parameters. The formal mathematical biology started with the work of Kermack and McKendrick [5], and then, significant contributions were made in the subject, for instance, one can see [6–11].To study the transmission of smoking habits and its control in the society, numerous researchers used the tools of mathematical modeling and optimal control theory. Castillo-Garrsow et al. [12] used a model of SIR type and well-presented the qualitative aspects of the model under discussion. Sharomi and Gumel [13] extended the work of Castillo-Garrsow et al. and introduced temporary quit class into the model. Zaman [14] incorporated the light smoker compartment and rigorously investigated the model from mathematical perspectives. Subsequently, considering various features and characteristics of smoking, sophisticated models have been developed and investigated, for example, [15–18].Incidence rate plays a crucial role while investigating the dynamics of any epidemic disease. The bilinear and saturated incidence rates were widely used in case of epidemic diseases and smoking models, for instance, one can see [12, 19, 20]. The authors in [17, 21, 22] used the square root incidence rate and discussed the models from various mathematical aspects. In the sequel, Rahman et al. [4] used harmonic mean type of incidence rate in his giving-up smoking model, investigated the dynamics of the model, and set the controlling strategies for reducing the smoking habit. Similarly, Alzaid and Alkahtani used the same incidence rate and studied the effect of relapse [23]. This work as well as the work of Rahman et al. assumed that potential smokers will start smoking only if they make a contact with light smoker. However, this is usually not the case. A person (potential smoker) may start smoking when he/she comes into contact with a smoker. In this work, the authors intend to overcome this gap. We shall assume that smoking habit may spread in the population via two roots; (i) contacts of potential smoker and light smokers and (ii) contact between smoker and light smokers. For the former type of spread, we will use the usual harmonic mean type of incidence rate, and for the later, the standard bilinear incidence rate will be utilized. By doing so, the present model will cover models [4, 23] as a subcase.The rest of the work is organized as follows. The model formulation and essential biological features to the model are discussed in Section2. The smoking generation number and equilibria of the model are presented in Section 3. We investigated the local stability of each equilibrium point in Section 4, whereas the criteria for global dynamics of smoking-free, light smoker present, and smoking-present equilibria are derived in Section 5. The local sensitivity will be performed in Section 6. Based on sensitivity analysis, we take into account some control variables and formulated a control problem for further analysis in Section 7. The desired goals were obtained, and the results are verified through simulations in Section 8. Finally, we presented the conclusion of the work.


## 2. Model Formulation and Well-Posedness

To formulate the model, we will divide the entire population into four compartments, namely, the potential smokersP=Pt, the light (or occasional) smokers L=Lt, the smokers S=St, and the quit smokers Q=Qt. That is, if Tt denotes the total population, and then, Tt=P+L+S+Q. We assumed that smoking habit spread in the population through (i) contact between potential smokers and light smokers and (ii) contact between the light smokers and chain smokers. The former contact is mathematically described by harmonic mean type of incidence rate. The reason for taking such type of contact is that the potential smokers may start smoking at a slower rate whenever they come into contact with the light smokers. Due to tendency of light smokers towards smoking, the light smokers will come into contact with smokers frequently and the rate will obey the mass action law. Hence, we assumed the bilinear incidence rate between light smokers and smokers. Furthermore, the harmonic mean type of incidence rate takes into account the behavioral changes of potential smokers and the crowding effect of the light smokers which prevent the unboundedness of the contact rate by choosing suitable parameters. Also, it is very hard to become a light smoker (because many social/cultural variables can restrict the contacts of light smokers with potential smokers), and mathematically, this phenomenon can be modeled by using the harmonic mean type of incidence rate instead of using bilinear incidence rate (because PL/P+L≤PL whenever P and L are positive). Besides this type of spread, we take into account the spread arising from contact between potential smokers and smokers which is described by bilinear incidence rate. It is also assumed that the quitting of smoking is not permanent, and an individual may start smoking again. These assumptions lead to the following system of ODEs representing the dynamics of smoking habit in a population.
(1)dPtdt=Λ−2β1PtLtPt+Lt−d+μPt+γQt,dLtdt=2β1PtLtPt+Lt−β2LtSt−d+μLt,dStdt=β2LtSt−d+μ+δSt,Qtdt=δSt−d+μ+γQt,with initial conditions
(2)P0>0,L0≥0,S0≥0,Q0≥0,where P0, L0, S0, and Q0 stand for the initial size of the population of the respective compartments. The description of parameters used in model (1) is shown in Table 1.Table 1
Interpretation of parameters of model (1).
ParameterInterpretationΛConstant recruitment into the potential smokersβ1The contact rate of occasional and potential smokersβ2The contact rate between smokers and light smokersμThe natural mortality rated0Death rate due to smokingδTemporarily quitting rate of smokersBy adding the governing equations of model (1), we obtain the conservation law
(3)dTtdt=Λ−d0+μTt,with T0=T0=P0+L0+S0+Q0.Theorem 1.
For the nonnegative initial data, if a solution to model (1) exists, it must be positive.Proof.
Consider the second equation in model (1), and its solution is given by
(4)Lt=L0e2β1Pt/Pt+Lt−β2St−d+μ≥0.
By using this relation in the third equation of system (1), we have St≥0 which assures the positivity of Qt. Finally, using these facts in the first equation, we have Pt≥0.Theorem 2.
For nonnegative initial conditions that are not identically zero on any interval, all possible solutions of system (1) are bounded in the region
(5)D=Pt,Lt,St,Qt∈R+4:Pt+Lt+St+Qt≤Λμ.Since each compartment denotes the size of population and thus must be nonnegative, the desired solution space isD=Pt,Lt,St,Qt∈R+4. Thus, the each component in the vector Pt,Lt,St,Qt is bounded below by zero, and thus, we have to find the upper bound for the solution. From the conservation equation (3), we have
(6)dTtdt=Λ−d0+μTt≤Λ−μTt.By solving this differential inequality, we get(7)Tt≤T0e−μt+Λμ.Lettingt⟶∞, we get Tt≤Λ/μ. Thus, all solutions of the model are bounded by 0 and Λ/μ, and hence, the desired feasible set is (4).Theorem 3.
For system (1), the nonnegative space R+4 is positive invariant.Proof.
Consider the vectorΨ=P,L,S,Q; then, in matrix form, we can write model (1) in the form
(8)dΨtdt=AΨ+B,where
(9)A=−2β1P+LL+d0+μ00γ2β1P+LL−d0+μ−β2S000β2S−d0+μ+δ000δ−d0+μ+γ,B=Λ000.
Clearly, the nondiagonal entries of matrixA are nonnegative which guarantees that A is Metzler matrix [24]. Further, B≥0, and thus, system (1) is positively invariant in the desired space. In the next section, we intended to find the smoking generation number (an important parameter describing the behavior of a system) as well as the possible equilibria of the proposed system.


## 3. Equilibria of the Model and Smoker Generation Number

In order to find the equilibria of the proposed model (1), we set
(10)Λ−2β1PtLtPt+Lt−d+μPt+γQt=0,(11)2β1PtLtPt+Lt−β2LtSt−d+μLt=0,(12)β2LtSt−d+μ+δSt=0,(13)δSt−d+μ+γQt=0.System of equations (10)–(13) always has a solution of the form
(14)E0=Λμ+d,0,0,0.This fixed point is known as smoking-free equilibrium of the model.To calculate the smoker generation number, we will follow the procedure of [25]. For this purpose, consider the smoker classes from system (1) and assume y=Lt,StT, that is,
(15)dydt=F−W,where
(16)F=2β1P+LPL−β2LSβ2LS,W=d+μLtd+μ+δSt.Considering the Jacobians of these two matrices at DFE, we get(17)F~=2β1000,W~=d+μ00d+μ+δ,and thus,
(18)W~−1=1d+μ001d+μ+δ.Let us considerG=F~W~−1; thus,
(19)G=2β1d+μ000.SinceR0 (the smoking generation) is the dominant eigenvalue of the next generation matrix G, hence
(20)R0=2β1d+μ.By using the numberR0, we can obtain other possible equilibria of system (1) using equation (3).Theorem 4.
(1)
IfR0>1, then model (1) has a positive light smoker present equilibrium given by(21)El=Pl,Ll,Sl,Ql=Λ2β1,Λ2β1R0−1,0,0.(2)
IfR0<1, then the model has a smoking-free equilibrium given by (14)(3)
WheneverR0=1, then again, the model has the smoking-free equilibrium given by (14)Proof.
By solving systems (10)–(13) and keeping in view that there are no smokers (i.e., St=0) in the community, we obtained the light smoker present equilibrium (21). However, for R0<1, the light smoker becomes negative and hence, the model will have just the smoking-free equilibrium (14). In the similar way, if we assume R0=1, then the fixed point (21) becomes the smoking-free equilibrium as we have Λ/2β1=Λ/μ+dR0.Theorem 5.
There exists a positive smoking present equilibrium of model (1) blue if R0≥1 andβ2S∗<μ+d.Proof.
System (1) has a positive solution of the form
(22)E∗=P∗,L∗,S∗,Q∗,where
(23)P∗=Λμ+d+γδS∗μ+dμ+d+γ−μ+d+δβ2μ+dβ2S∗−μ+d,L∗=μ+d+δβ2,Q∗=δS∗μ+d+γ,and S∗ is a positive root of
(24)β2d+μ+γ+δd+μ+γS∗2−d+μR0−1d+μ+γ+δd+μ+γ+Λβ2μ+dS∗+ΛR0−1+R0μ+dμ+d+δβ2=0,
which surely exist by Descartes’ rule of signs wheneverR0>1 and β2S∗<μ+d. Further, if we set R0=1, then both solutions of equation (24) are positive and real and hence the theorem.


## 4. Local Stability Analysis of Equilibria

To find out the local stability of system (1) at each equilibrium point, first, we will calculate the Jacobian matrix
(25)J=−2β1L2L+P2−d+μ−2β1P2L+P20γ2β1L2L+P22β1P2L+P2−β2S−d+μ−β2L00β2Sβ2L−d+μ+δ000δ−d+μ+γ.Theorem 6.
The SFE (14) of system (1) is locally asymptotically stable for R0<1.
From relation (25), the Jacobian of system (1) at E0 is given by
(26)JE0=−d+μ2β10γ02β1−d+μ0000−d+μ+δ000δ−d+μ+γ.Matrix (26) has eigenvalues t1=−μ+d, t2=−μ+d+δ, t3=2β1−μ+d, and t4=−μ+d+γ. As the parameters of the model assume positive values, thus the eigenvalues t1, t2, and t4 are negative and t3=μ+d2β1/μ+d−1=μ+dR0−1 is negative only if R0<1. Hence, for R0<1, all of the eigenvalues of the Jacobian matrix at E0 are negative and so E0 is locally asymptotically stable. Further, if R0=1, then t3=0. However, if we assume R0>1, then E0 becomes unstable as one of the eigenvalues is positive in such case.Theorem 7.
The local smoking present equilibrium (LSPE) (21) of system (1) is locally asymptotically stable for 1<R0<1+2β1μ+d+δ/Λβ2.Proof.
From relation (25), the Jacobian of system (1) at El is given by
(27)JEl=−2β1Ll2Ll+Pl2−d+μ−2β1Pl2Ll+Pl20γ2β1Ll2Ll+Pl22β1Pl2Ll+Pl2−d+μ−β2Ll000β2Ll−d+μ+δ000δ−d+μ+γ.
The eigenvalues of matrix (30) are t1=−μ+d,t2=−μ+dR0−1, t3=−μ+d+γ, and t4=R0−1Λβ2/2β1−μ+d+δ. Now, for positive parameters of the model and for R0>1, the eigenvalues t1, t2, and t3 are negative. However, t4 is negative only for R0<1+2β1μ+d+δ/Λβ2. Hence, if 1<R0<1+2β1μ+d+δ/Λβ2, all of the eigenvalues of the Jacobian matrix at El are negative and so El is locally asymptotically stable. Moreover, for R0=1, we have t2=0. In the case of R0<1, t2 becomes positive and thus, the equilibrium El will be unstable.Theorem 8.
The smoking present equilibrium (SPE) is locally asymptotically stable of if1<R0<1+β2S∗/μ+d.Proof.
The Jacobian at SPE is given by(28)J=−a11−a120γa21−a22−β2L∗00a32β2L∗−d+μ+δ000δ−a44,where
(29)a11=2β1L∗2L∗+P∗2+d+μ,(30)a12=2β1P∗2L∗+P∗2,(31)a21=2β1L∗2L∗+P∗2,(32)a22=β2S∗+d+μ−2β1P∗2L∗+P∗2,(33)a32=β2S∗,(34)a44=d+μ+γ.
The characteristic equation of matrix (28) is given by
(35)τ4+a44+a22+a11τ3+a32d+δ+μ+a11a22+a44+a12a21+a22a44τ2+da32+δa32+μa32×a11+a44+a44a11a22+a12a21τ+d+δa11a32a44−δa32a21γ+μa11a32a44=0.
Clearly,a11, a12, a21, a32, and a44 are positive and
(36)a22=β2S∗+d+μ−2β1P∗2L∗+P∗2,≥β2S∗+d+μ−2β1,=β2S∗−d+μR0−1>0,
only ifβ2S∗>d+μR0−1. Thus, all of the coefficients of the characteristic equation (35) are positive, and hence, by Descartes’ rule of signs, this equation has no positive real solution. By using −τ in place of τ in equation (35) and utilizing Descartes’ rule of signs, we get that all of the eigenvalues of this equation are negative or complex conjugates with dominant negative real part. Therefore, the SPE is locally asymptotically stable whenever R0<1+β2S∗/μ+d. Further, for the existence of SPE, we must assume that R0>1. Moreover, the case R0=1 does not affect the stability of the equilibrium point; however, this condition will challenge the existence of smoking-present equilibrium.


## 5. Global Stability Analysis of the Equilibria

Theorem 9.
LetR0≤1; then, the smoking-free equilibrium (SFE) E0 of model (1) is globally asymptotically stable (GAS).Proof.
For proving the required result, the Lyapunov function of the following form is assumed:(37)V1P,L,S,Q=Lt+St+Qt.
By considering the derivative of (37) with respect to time t and using model (1), we have
(38)dV1dt=dLtdt+dStdt+dQtdt,=2β1PtLtPt+Lt−β2LtSt−d+μLt+β2LtSt−d+μ+δSt+δSt−d+μ+γQt,≤2β1Lt−d+μLt−d+μ+δSt+δSt−d+μ+γQt,because Pt/Pt+Lt≤1. Thus,
(39)dV1dt≤2β1−d+μLt−d+μSt−d+μ+γQt,=μ+dR0−1Lt−d+μSt−d+μ+γQt.
Hence,(40)dV1dt<0,only if R0≤1, and thus, by LaSalle’s invariant principle [26], E0 is globally asymptotically stable in the feasible region.Theorem 10.
If1<R0<β2Λ+2β1d+μ/β2Λ, then the local smoking present equilibrium (LSPE) El of model (1) is globally asymptotically stable (GAS).Proof.
To prove the main result, we will take into consideration the Lyapunov function of the form(41)V2P,L,S,Q=L−Ll−LllnLtLl+St+Qt.
By considering the derivative of (37) with respect to time t and using model (1), we have
(42)dV2dt=dLtdt−LlLtdLtdt+dStdt+dQtdt,=2β1PtLtPt+Lt−β2LtSt−d+μLt−LlLt2β1PtLtPt+Lt−β2LtSt−d+μLt+β2LtSt−d+μ+δSt+δSt−d+μ+γQt,≤μ+d1−R0Lt−μ+dR0−1Ll+β2Ll−d+μSt−d+μ+γQt,<0,for R0>1 and so by LaSalle’s invariant principle, the light smoker present equilibrium El is globally asymptotically stable. By considering relation (42) and the case of R0=1, we will reach to the conclusion of Theorem 9 (as the case of R0=1 does not guarantee the existence of LSPE).Theorem 11.
IfR0>1, β2S∗<μ+d and R0P∗2/P∗+L∗2<1, then (1) is GAS at SPE E∗.Proof.
By considering the Jacobean (25) of the system at the SPE, we have
(43)J∗=−d+μ+2β1L∗2P∗+L∗2−2β1P∗2P∗+L∗20γ2β1L∗2P∗+L∗2R0P∗2P∗+L∗2−1d+μ−β2S∗−β2L∗00β2S∗β2L∗−d+μ+δ000δ−d+μ+γ.
The third additive compound matrix ofJ∗ is given by
(44)J∗∣3∣=A1100γδA22−β2L∗00β2S∗A33−2β1P∗2P∗+L∗2002β1L∗2P∗+L∗2A44,
where(45)A11=−d+μ+2β1L∗2P∗+L∗2+R0P∗2P∗+L∗2−1d+μ−β2S∗+β2L∗−d+μ+δ,A22=−d+μ+2β1L∗2P∗+L∗2+R0P∗2P∗+L∗2−1d+μ−β2S∗−d+μ+γ,A33=−d+μ+2β1L∗2P∗+L∗2+β2L∗−d+μ+γ−d+μ+δ,A44=R0P∗2P∗+L∗2−1d+μ−β2S∗+β2L∗−d+μ+γ−d+μ+δ.
Considerρχ=diagQ,S,L,P such that ρχ−1=diag1/Q,1/S,1/L,1/P and time derivative is ρfχ=diagQ̇,Ṡ,L̇,Ṗ. Therefore,
(46)ρfρ−1=diagṖP,L̇L,ṠS,Q̇Q,and hence,
(47)ρJ∣3∣ρ−1=A1100γQPδSQA22−β2S00β2LA33−2β1LPP+L2002β1LPP+L2QSA44.
By definingB=ρfρ−1+∣2∣ρ−1 so that
(48)B=b1100b14b21b22b2300b32b33b3400b43b44,where
(49)b11=Q̇Q−d+μ+2β1L∗2P∗+L∗2+R0P∗2P∗+L∗2−1d+μ−β2S∗+β2L∗−d+μ+δ,b14=γQP=ṖP−ΛP+2β1LP+L+d+μ,b21=δSQ=Q̇Q+d+μ+γ,b22=ṠS−d+μ+2β1L∗2P∗+L∗2+R0P∗2P∗+L∗2−1d+μ−β2S∗−d+μ+γ,b23=−β2S,b32=β2L,b33=L̇L−d+μ+2β1L∗2P∗+L∗2+β2L∗−d+μ+δ−d+μ+γ,b34=−2β1LPP+L2,b43=2β1LPP+L2,b44=ṖP+R0P∗2P∗+L∗2−1d+μ−β2S∗+β2L∗−d+μ+δ−d+μ+γ.
Consequently,(50)h1t=b11+∑j=24b1j,h1t=Q̇Q−d+μ+2β1L∗2P∗+L∗2+R0P∗2P∗+L∗2−1d+μ−β2S∗+β2L∗−d+μ+δ+γQP,≤ṖP+Q̇Q+R0P∗2P∗+L∗2−1d+μ,h2t=b22+∑j=1andj≠24b2j,h2t=ṠS−d+μ+2β1L∗2P∗+L∗2+R0P∗2P∗+L∗2−1d+μ−β2S∗−d+μ+γ+δSQ+β2S+Q̇Q+d+μ+γ,≤ṠS+Q̇Q+R0P∗2P∗+L∗2−1d+μ,h3t=b33+∑j=1andj≠34b3j,h3t=L̇L−d+μ+2β1L∗2P∗+L∗2+β2L∗−d+μ+δ−d+μ+γ+β2L+2β1LPP+L2,≤L̇L−d+μ+γ,h4t=b44+∑j=13b4j,h4t=ṖP+R0P∗2P∗+L∗2−1d+μ−β2S∗+β2L∗−d+μ+δ−d+μ+γ+2β1LPP+L2,≤ṖP+R0P∗2P∗+L∗2−1d+μ.
Now we have(51)limt⟶∞supsup1t∫lim0th1tdt<limt⟶∞1tlogPtP0+limt⟶∞1tlogQtQ0+R0P∗2P∗+L∗2−1d+μ,<R0P∗2P∗+L∗2−1d+μ,limt⟶∞supsup1t∫lim0th2tdt<limt⟶∞1tlogStS0+limt⟶∞1tlogQtQ0+R0P∗2P∗+L∗2−1d+μ,<R0P∗2P∗+L∗2−1d+μ,limt⟶∞supsup1t∫lim0th3tdt<limt⟶∞1tlogLtL0−d+μ+γ,<−d+μ+γ,limt⟶∞supsup1t∫lim0th4tdt<limt⟶∞1tlogPtP0+R0P∗2P∗+L∗2−1d+μ,<R0P∗2P∗+L∗2−1d+μ.
By combining the preceding four inequalities, we have the following inequality.(52)q¯=limt⟶∞supsup1t∫0tμBdt<0,and we denote the Lozinskii measure by μB=hi,i=1,2,3,4, and hence, it shows that the SPE is GAS.


## 6. Local Sensitivity Analysis

The termR0 is generally influenced by the inconsistencies in data gathering and estimated values. Sensitivity analysis is used to assess the relative impact of epidemic factors for disease propagation and control.Definition 12.
For the basic reproduction number, the normalized sensitivity index of a parameterψ (which depends on the partial derivative of R0w.r.tψ) is of the following form:
(53)Sψ=ψR0∂R0∂ψ.By calculating the sensitivity index, we may determineR0’s responsiveness. To find the sensitivity of each parameter for R0, we will employ equation (53) which was suggested by Tilahun et al. [27]. On the basis of parameter values β1=0.006, d=0.00004, and μ=0.08, we have the sensitivity indices of the form
(54)Sβ1=β1R0∂R0∂β1=β1R02μ+d=1,Sμ=μR0∂R0∂μ=−μμ+d=−0.999500<0,Sd=dR0∂R0∂d=−dμ+d=−0.00049975<0.These indices enable us to determine the relevance of numerous variables involved in disease transmission, as well as the relative change in the number of reproductions as a function of parameter changes. Using such indices, we may identify the characteristics that have a significant impact onR0 and are critical for disease prevention.Table2 suggests that there is a positive relation between R0 and β1 while negative influence of μ and d on R0. The same argument is supported by Figure 1. This means that if you raise or reduce the value of parameter β1 by 10%, it will rise or decrease R0 by 10%. Similarly, if we reduce or increase μ and d by 1 percent, its relative inverse impact on R0 will be 99.9% and 0.04%, respectively. To prevent smoking from the population, we must consider those parameters which have high sensitivity indices. Based on the sensitivity indices, we concluded that β1 will directly affect R0 (100%) whereas μ affects R0 inversely about 99%. Thus, to reduce smoking habit, we must focus on the contacts between potential smokers and light smokers. It is simple to develop a control mechanism program for eradication of smoking habit based on this study.Table 2
Sensitivity indices ofR0 w.r.t and the parameters β1, μ, and d.
ParametersSensitivity indexValueβ1Sβ11μSμ−0.999500dSd−0.00049975Figure 1
Sensitive parametersβ, μ, and d and its impact on R0.
(a)
R0 verses μ and d at fixed β1=0.006(b)
R0 verses β and d at fixed μ=0.08


## 7. Formulation and Analysis of Optimal Control Problem

In this part, we simulate the eradication of smoking habit from the population by using the tools of control theory [9, 28]. We focus on reducing transmission, which has a sensitivity index of 1, to develop a control plan based on local sensitivity analysis. Further, to make the control program more effective, we take into consideration the following four control measures.The first control measure is the education campaign which is denoted byu1t, and its aim is to reduce the size of potential smokers. The second control variable u2t shows physically antismoking gum, and it will reduce the number of light smokers. The third and fourth control variables are antinicotive drug u3t and ban on smoking particularly in public places by the government u4t. These measures will be imposed on the smoking compartment. It is observed that whenever the law enforcement personnel came into contact with light smokers (which smoke in public places), the light smokers will tend to quit smoking in public places. That is, by increasing the law enforcement personnel in public places, the size of light smokers will tend to decrease and hence, small number of individuals will tend to quit smoking. Therefore, we will use the term 1−u4t which will reduce the contacts between the law enforcement personnel with the light smokers. The control variables are subject to some conditions; uit∈0,uimax for i=1,⋯,4 and u1≤u4≤u2≤u3. These control variables are the main measures which could help in reducing the smoking habit [29]. Based on the control measures, we have the following governing equations for the control problem:
(55)dPtdt=Λ−2β1PtLtPt+Lt−d+μ+u1tPt+γQt,dLtdt=2β1PtLtPt+Lt−β2LtSt−d+μ+u2t+1−u4tLt,dStdt=β2LtSt−d+μ+δ+u3t+qu4tSt,Qtdt=δSt−d+μ+γQt+u3t+qu4tSt+u2t+1−u4tLt+u1tPt,with the conditions(56)P0>0,L0≥0,S0≥0,Q0≥0,Our optimum control strategy is to reduce/minimize the number of potential smokers, light smokers, and smokers while increasing the number of people who quit smoking. Keeping in view the control problem (55), our cost functional is given by the following:
(57)Ju1,u2,u3,u4=∫0TD1L+D2S+D3P−D4Q+12F1v12t+F2v22t+F3v32t+F4v42tdt.In the cost functional (57), Di’s are positive constants describing the balancing factors, whereas Fi’s are the cost associated to the control measures. Clearly, the cost functional has a goal to reduce smoker population and to enhance the size of quit smokers. Our main focus is to find a set of functions u1⋆,u2⋆,u3⋆,u4⋆ in such a way that
(58)Ju1⋆,u2⋆,u3⋆,u4⋆=minui∈VJui,where the set V is the admissible control set and is define by
(59)V=uit∈L20,T:0≤uit≤uimax≤1fori=1,⋯,4,where uit are the control variables described above, and these are Lebesgue measurable functions. To identify such control measures, we must first establish that actually they exist.


### 7.1. Existence of Solution to the Control Problem

For proving the desired result, we will take into consideration problems (55) and (56) and will show that actually this system has a solution. It is worthy to notice that nonnegative bounded solutions to the state system exist if we have bounded Lebesgue measurable control measures and nonnegative initial conditions [28]. Let
(60)dψdt=Lψ+Xψ,and here,
(61)ψ=PtLtStQt,L=−d+μ+u1t00γ0−d+μ+u2t+1−u4t0000−d+μ+δ+u3t+pu4t0u1t1−u4tδ+u3t+pu4t−d+μ+γ,Xψ=Λ−2β1PtLtPt+Lt2β1PtLtPt+Lt−β2LtStβ2LtSt0.Equation (60) denotes a nonlinear system (the proposed control problem) with bounded coefficients. Set
(62)K=Lψ+Xψ.By considering the following,(63)Xψ1−Xψ2≤2β1P1−P2+2β1L1−L2+β2L1S1−L2S2,≤N1P1−P2+L1−L2+L1S1−L2S2,whereN1=max2β1,β2 is free from the state’s variable. Also, one can write
(64)Kψ1−Kψ2≤≤M1ψ1−ψ2,where ∞>M1=maxN1,L; thus, X is uniformly continuous function in the sense of Lipschitz. Further, for the definition of the control measures and conditions on the states (Pt>0,Lt≥0,St≥0 and Rt>0), we can observe that a solution to the control problem (55) does really exist. The following conclusion holds for the existence of control variables in the optimum control problem.Theorem 13.
There exists a control vectoru⋆=u1⋆,u2⋆,u3⋆,u4⋆∈U which minimizes the objective functional.Proof.
In order to show that actually such control variables exist, we need to follow [28], as
(a)
Both the state and control functions are nonnegative(b)
Convexity and closedness properties hold by the set of admissible controls(c)
The boundedness of the control model assures the compactness(d)
The function inside the integral in the objective functional (57) is convex in the control measures
Therefore, there exist the control variablesuit for i=1,2,3,4 which minimize the objective functional.


### 7.2. Optimality Conditions

For deriving characterization of the control from the control problem (55) subject to the cost functional (57), first of all, we will define the Lagrangian and the Hamiltonian for this problem. Let z=P,L,S,Q denote the vector whose components are the state functions and u=u1,u2,u3,u4 is the control vector. The Lagrangian L is given by
(65)Lz,u=D1L+D2S+D3P−D4Q+12F1v12t+F2v22t+F3v32t+F4v42t,and the Hamiltonian H is given by
(66)Hz,u,λ=Lz,u+λfz,u,whereλ=λ1,λ2,λ3,λ4 and fz,u=f1z,u+f2z,u+f3z,u+f4z,u with
(67)f1z,u=Λ−2β1PtLtPt+Lt−d+μ+u1tPt+γQt,f2z,u=2β1PtLtPt+Lt−β2LtSt−d+μ+u2t+1−u4tLt,f3z,u=β2LtSt−d+μ+δ+u3t+qu4tSt,f4z,u=δSt−d+μ+γQt+u3t+qu4tSt+1−u4tLt+u1tPt.Thus,(68)Hz,u,λ=D1L+D2S+D3P−D4Q+12F1v12t+F2v22t+F3v32t+F4v42t+λ1Λ−2β1PtLtPt+Lt−d+μ+u1tPt+γQt+λ22β1PtLtPt+Lt−β2LtSt−d+μ+u2t+1−u4tLt+λ3β2LtSt−d+μ+δ+u3t+qu4tSt+λ4δSt−d+μ+γQt+u3t+qu4tSt+u2t+1−u4tLt+u1tPt.The Pontryagin maximum principle [28] is used to determine the best solution for our given optimum control problem. If z⋆,u⋆ is the best solution for the control problem (55), then there exists a nonzero vector function λ such that the Hamiltonian system satisfies
(69a)dz⋆tdt=∂H∂λz⋆t,u⋆t,λt,(69b)∂H∂uz⋆t,u⋆t,λt=0,(69c)dλtdt=−∂H∂zz⋆t,u⋆t,λt,while holding the maximality and transversality conditions
(70)Hz⋆t,u⋆t,λt=maxuHz⋆t,u,λt,(71)λT=0.Theorem 14.
LetP⋆, L⋆, S⋆, and Q⋆ be optimal states associated with the optimal controls u1⋆,u2⋆,u3⋆,u4⋆ for problem (55). Then, there exist adjoint variables λ1t, λ2t, λ3t, and λ4t satisfying
(72)dλ1tdt=2β1L2tPt+Lt2+d+μ+u1tλ1t−2β1L2tPt+Lt2λ2t−u1tλ4t−D3,(73)dλ2tdt=2β1P2tPt+Lt2λ1t−2β1P2tPt+Lt2λ2t+β2Stλ2t+d+μ+u2t+1−u4tλ2t−β2Stλ3t−u2t+1−u4tλ4t−D1,(74)dλ3tdt=β2Ltλ2t−β2Ltλ3t+d+μ+δ+u3t+qu4tλ3t+δλ4t−u3t+qu4tλ4t−D2,(75)dλ4tdt=−γλ1t+d+μ+γλ4t+D4,with terminal conditions
(76)λ1T=0,λ2T=0,λ3T=0,λ4T=0.Further, the control measures are characterized by(77)u1⋆t=maxminc1,λ1t−λ4tF1Pt,0,(78)u2⋆t=maxminc2,λ2−λ4F2Lt,0,(79)u3⋆t=maxminc3,λ3−λ4F3St,0,(80)u4⋆t=maxminc4,λ4−λ2Lt+qλ3−λ4StF4,0.Proof.
The adjoint system (72) is obtained by using the adjoint equation (69c) in the maximum principles. The terminal conditions of (76) were derived as a direct consequence of equation (71). We took the partial derivatives of H with respect to the control measure ui for i=1,⋯,4 in turns and used (69b) for the derivation of (77)–(80).


## 7.1. Existence of Solution to the Control Problem

For proving the desired result, we will take into consideration problems (55) and (56) and will show that actually this system has a solution. It is worthy to notice that nonnegative bounded solutions to the state system exist if we have bounded Lebesgue measurable control measures and nonnegative initial conditions [28]. Let
(60)dψdt=Lψ+Xψ,and here,
(61)ψ=PtLtStQt,L=−d+μ+u1t00γ0−d+μ+u2t+1−u4t0000−d+μ+δ+u3t+pu4t0u1t1−u4tδ+u3t+pu4t−d+μ+γ,Xψ=Λ−2β1PtLtPt+Lt2β1PtLtPt+Lt−β2LtStβ2LtSt0.Equation (60) denotes a nonlinear system (the proposed control problem) with bounded coefficients. Set
(62)K=Lψ+Xψ.By considering the following,(63)Xψ1−Xψ2≤2β1P1−P2+2β1L1−L2+β2L1S1−L2S2,≤N1P1−P2+L1−L2+L1S1−L2S2,whereN1=max2β1,β2 is free from the state’s variable. Also, one can write
(64)Kψ1−Kψ2≤≤M1ψ1−ψ2,where ∞>M1=maxN1,L; thus, X is uniformly continuous function in the sense of Lipschitz. Further, for the definition of the control measures and conditions on the states (Pt>0,Lt≥0,St≥0 and Rt>0), we can observe that a solution to the control problem (55) does really exist. The following conclusion holds for the existence of control variables in the optimum control problem.Theorem 13.
There exists a control vectoru⋆=u1⋆,u2⋆,u3⋆,u4⋆∈U which minimizes the objective functional.Proof.
In order to show that actually such control variables exist, we need to follow [28], as
(a)
Both the state and control functions are nonnegative(b)
Convexity and closedness properties hold by the set of admissible controls(c)
The boundedness of the control model assures the compactness(d)
The function inside the integral in the objective functional (57) is convex in the control measures
Therefore, there exist the control variablesuit for i=1,2,3,4 which minimize the objective functional.


## 7.2. Optimality Conditions

For deriving characterization of the control from the control problem (55) subject to the cost functional (57), first of all, we will define the Lagrangian and the Hamiltonian for this problem. Let z=P,L,S,Q denote the vector whose components are the state functions and u=u1,u2,u3,u4 is the control vector. The Lagrangian L is given by
(65)Lz,u=D1L+D2S+D3P−D4Q+12F1v12t+F2v22t+F3v32t+F4v42t,and the Hamiltonian H is given by
(66)Hz,u,λ=Lz,u+λfz,u,whereλ=λ1,λ2,λ3,λ4 and fz,u=f1z,u+f2z,u+f3z,u+f4z,u with
(67)f1z,u=Λ−2β1PtLtPt+Lt−d+μ+u1tPt+γQt,f2z,u=2β1PtLtPt+Lt−β2LtSt−d+μ+u2t+1−u4tLt,f3z,u=β2LtSt−d+μ+δ+u3t+qu4tSt,f4z,u=δSt−d+μ+γQt+u3t+qu4tSt+1−u4tLt+u1tPt.Thus,(68)Hz,u,λ=D1L+D2S+D3P−D4Q+12F1v12t+F2v22t+F3v32t+F4v42t+λ1Λ−2β1PtLtPt+Lt−d+μ+u1tPt+γQt+λ22β1PtLtPt+Lt−β2LtSt−d+μ+u2t+1−u4tLt+λ3β2LtSt−d+μ+δ+u3t+qu4tSt+λ4δSt−d+μ+γQt+u3t+qu4tSt+u2t+1−u4tLt+u1tPt.The Pontryagin maximum principle [28] is used to determine the best solution for our given optimum control problem. If z⋆,u⋆ is the best solution for the control problem (55), then there exists a nonzero vector function λ such that the Hamiltonian system satisfies
(69a)dz⋆tdt=∂H∂λz⋆t,u⋆t,λt,(69b)∂H∂uz⋆t,u⋆t,λt=0,(69c)dλtdt=−∂H∂zz⋆t,u⋆t,λt,while holding the maximality and transversality conditions
(70)Hz⋆t,u⋆t,λt=maxuHz⋆t,u,λt,(71)λT=0.Theorem 14.
LetP⋆, L⋆, S⋆, and Q⋆ be optimal states associated with the optimal controls u1⋆,u2⋆,u3⋆,u4⋆ for problem (55). Then, there exist adjoint variables λ1t, λ2t, λ3t, and λ4t satisfying
(72)dλ1tdt=2β1L2tPt+Lt2+d+μ+u1tλ1t−2β1L2tPt+Lt2λ2t−u1tλ4t−D3,(73)dλ2tdt=2β1P2tPt+Lt2λ1t−2β1P2tPt+Lt2λ2t+β2Stλ2t+d+μ+u2t+1−u4tλ2t−β2Stλ3t−u2t+1−u4tλ4t−D1,(74)dλ3tdt=β2Ltλ2t−β2Ltλ3t+d+μ+δ+u3t+qu4tλ3t+δλ4t−u3t+qu4tλ4t−D2,(75)dλ4tdt=−γλ1t+d+μ+γλ4t+D4,with terminal conditions
(76)λ1T=0,λ2T=0,λ3T=0,λ4T=0.Further, the control measures are characterized by(77)u1⋆t=maxminc1,λ1t−λ4tF1Pt,0,(78)u2⋆t=maxminc2,λ2−λ4F2Lt,0,(79)u3⋆t=maxminc3,λ3−λ4F3St,0,(80)u4⋆t=maxminc4,λ4−λ2Lt+qλ3−λ4StF4,0.Proof.
The adjoint system (72) is obtained by using the adjoint equation (69c) in the maximum principles. The terminal conditions of (76) were derived as a direct consequence of equation (71). We took the partial derivatives of H with respect to the control measure ui for i=1,⋯,4 in turns and used (69b) for the derivation of (77)–(80).


## 8. Simulation Results

For solving the problems (both with and without controls), we used the standard RK method of fourth order. To verify numerically the key theorems on dynamical analysis, we utilized parameter values from Tables3 and 4. Specifically, values of parameters from Table 3 were used to show the stability of smoking-free equilibrium, and sequentially, parameter values from Table 4 guarantee the stability of light smoker present equilibrium. For presenting the effect of control measures, we assumed values from Table 5. For simulating the control problem, firstly, we solved model (55) forward in time and then used backward RK4 method for solving the adjoint system (72) with the help of terminal conditions and characterization of the control variables (77)–(80). The simulation clearly illustrates the results on dynamical analysis (Figures 2 and 3) as well as on control theory (Figures 4 and 5).Table 3
Parameter values for verifying the global stability of smoking-free equilibrium.
ParametersValuesParametersValuesΛ10.25μ0.08γ0.006β10.008β20.00038d0.0019δ0.000274Table 4
Parameter values for verifying the global stability of light smoker present equilibrium.
ParametersValuesParametersValuesΛ10.25μ0.0111γ0.006β10.038β20.00038d0.0019δ0.041Table 5
Parameter values are used for showing the effectiveness of the control program.
ParametersValuesParametersValuesΛ10.25μ0.0111γ0.006β10.9β20.00038d0.0019δ0.000274q0.2Figure 2
Different initial size of compartments and its long-term behavior forR0<1.
(a)
Dynamics of potential smokers wheneverR0<1(b)
The case ofR0<1 and evolution of light smokers(c)
Behavior of smoker population forR0<1Figure 3
The plot describes stability of smoking present equilibriumEå when R0>1.
(a)
The plot shows thatP⟶P∗ when R0>1(b)
Behavior of the light smokers forR0>1Figure 4
The plot presents a comparison of solution curves both with and without control.
(a)
The plot shows the overall effect of control variables on the class of potential smokers(b)
Comparison of curvesLt both with and without control(c)
Dynamics of smokers in the presence and absence of control variablesFigure 5
The plot describes stability of light smoker present equilibrium whenR0>1.
(a)(b)(c)By using values of parameter from Table3, calculate the smoking generation number R0=0.1954 and SFE E0=125.1526,0,0,0. It is clear from Figure 2 that each solution curve in the subplots tends to its SFE irrespective of the initial size of the population whenever R0<1.To show numerically the global stability of light smoker present equilibriumEl (21), we assumed the same values of Table 3 except β1=0.08 which gives R0=1.9536 and as a result, El=Pl,Ll,Sl,Ql=64.0625,61.0901,0,0. Thus, Figures 6(a)–6(d) shows a clear interpretation of Theorems 7 and 10, that is, if R0>1, then P,L,S,Q⟶Pl,Ll,Sl,Ql as t⟶∞.Figure 6
The dynamics of control variables.Similarly, to show the long-term behavior of each class wheneverR0>1 and the side condition of Theorem 8 holds, we take into consideration values from Table 4. From these values of parameters, we calculated R0=12.3077>1 and P∗,L∗,S∗,Q∗=127,142.1,164.5,354.9. Figure 3 proves the statement of Theorem 8 as well as the global stability of SPE numerically.Finally, we used values from Table5 and simulated both with and without control problems and presented the effect of the control variables. The simulation (in Figure 4) depicts in a clear way the effect of the control variables: to decrease the size of smoker population and to increase the number of quit smokers.


## 9. Concluding Remarks

In this work, we formulated and analyzed a giving-up smoking model utilizing three main factors related to smoking habit: the bilinear incidence rate showing the spread of smoking habit within the population from the contact of light smokers and smokers, the harmonic mean type of incidence rate (for incorporating contacts between potential smokers and occasional smokers), and the relapse factor associated to smoking. After model formulation, we investigated the model for bounded and positive solution subject to appropriate initial conditions. The model was checked for possible fixed points, and three equilibria were derived under certain conditions, namely, the smoking-free equilibrium, light smoker present equilibrium, and smoking present equilibrium. It was shown that the smoking-free fixed point always exists and is both locally and globally asymptotically stable whenR0<1. The positive light smoker present equilibrium exists if we assume R0>1, and the same fixed point is locally as well as globally asymptotically stable if we impose additional conditions on the basic reproduction number. Further, we proved that the model has a smoking present equilibrium whenever additional conditions are satisfied by R0. Also, a criteria for the local analysis has been derived, and the global stability of smoking present equilibrium was presented both analytically and numerically. Keeping in mind a few control measures, we set a control problem and optimal control strategies were achieved with the help of Pontryagin’s maximum principle. The obtained analytical results were verified through simulations and effectiveness of the control program is presented by comparing with and without control curves.This work is indeed a very good contribution to the existing literature on smoking epidemics. However, this social evil has a dramatic impact on every society, and thus, one cannot ignore the memory effect while studying such and related problems. Keeping in view the importance of generalized fractional derivatives in this regard [10–30], the authors have a keen interest to formulate and analyze the fractional models for smoking epidemics in near future.



---

*Source: 3771137-2022-04-14.xml*