Sensitivity and mathematical model analysis on secondhand smoking tobacco

Abstract

In this paper, we are concerned with a mathematical model of secondhand smoker. The model is biologically meaningful and mathematically well posed. The reproductive number $R_{0}$ is determined from the model, and it measures the average number of secondary cases generated by a single primary case in a fully susceptible population. If $R_{0}<1,$ the smoking-free equilibrium point is stable, and if $R_{0}>1,$ endemic equilibrium point is unstable. We also provide numerical simulation to show stability of equilibrium points. In addition, sensitivity analysis of parameters involving in the dynamic system of the proposed model has been included. The parameters involving in reproductive number measure the relative change in $R_{0}$ when the value of the parameter changes.

Introduction

Secondhand smoke is defined as “the combination of smoke emitted from the burning end of a cigarette or other tobacco products and smoke exhaled by the smoker” [1]. Thus, secondhand smoke exposure consists of an unintentional inhalation of smoke that occurs close to people smoking and/or in indoor environments where tobacco was recently used.

Subjects may be exposed to secondhand smoke in multiple sites, such as the home [2], public places, cars, homes of relatives [3], and the workplace [4]. Exposure to secondhand smoke is a significant risk factor for a plethora of diseases and adverse health-related outcomes at a global scale [5] and is estimated to cause over 600,000 deaths annually, the majority of which are due to ischemic heart disease among adults [6].

There is a connection between secondhand smoke and an increased risk of stroke. Regular exposure to secondhand smoke, such as in restaurants, increases one’s chance of stroke by 50 % [7]. Secondhand smoke is as damaging to a fetus as if the mother were inhaling the smoke directly from a cigarette [8]. Just 30 min of exposure to secondhand smoke can cause heart damage similar to that of habitual smokers [9].

Due to these, secondhand smoke exposure is a topic of great concern for us because of its well-known adverse effects on human health [10]. Worldwide, 40% of children, 33% of male nonsmokers and 35% of female nonsmokers were exposed to secondhand smoke. The highest proportions exposed were estimated in Europe, the Western Pacific, and Southeast Asia, with more than 50% of population exposed. Proportion of people exposed was lowest in Africa [6]. The prevalence of secondhand smoke exposure among adolescents in Ethiopia is highest [11].

Different scholars discussed about secondhand smoke in relation to health problem. For instance,  [12] discussed the global health burden of secondhand smoke, and [13] emphasized a new and alarming consequence of smoking in indoor environments.

All studies mentioned above showed secondhand smoke is a big problem in the world as the result of leading peoples to diseases and death. There are a number of studies on mathematical model of smoking by taking different assumption. We direct refer the reader to the papers [14,15,16,17].

In this article, we will see the effect of secondhand smoke in the society. We will prepare a mathematical model and interpret the model graphically. The model, we will develop extends the work of Eduardo L. L. M and Maribel L. B in paper [18]. We consider death $(\nu P)$ of secondhand smoker by this infection in addition to natural death [19, 20], and we discuss sensitivity analysis of $R_{0}$. Different parameters will be used for the infection rate $\lambda PS$ from P to S and infection rate $\delta QS$ from Q to S. The infection rate $\delta QS$ from Q to S means the infection rate when relapse smoking, that is, restarting smoking after a period of abstinence. In [21], relapse rates within the first year of abstinence ranged from 60 to 90%, while 2 years of continuous cessation indicated a likelihood of 80% to maintain long-term abstinence. A research conducted in the USA has the following results [22].

• 1 The prevalence of relapse is 6.8%.

• 2 Prevalence and odds of relapse were higher among young people compared with elders.

• 3 Former smokers living in smoke-free homes had 60% lower odds of relapse compared with those living in homes that allowed smoking inside.

• 4 Odds of relapse were higher among never married, widowed, divorced and separated individuals, compared with the married group.

• 5 Continuous smoking cessation for 6 months or more significant decreased odds of relapse.

We organized this paper as follows. In Section 2, we prepare the proposed model. In Sections 3, 4 and 5, we discuss free equilibrium point, reproduction number, sensitivity analysis and endemic equilibrium point, respectively. In Section 6, we discuss stability of equilibrium solution. In Section 7, we see the numerical simulation for the equilibrium solution. Finally, conclusion, abbreviations and declaration are included.

Formulation of the model

Model assumption

We have the following assumptions in preparing the model.

• 1 There are deaths as a result of smoking and secondhand smoking.

• 2 The number of healthy and completely recovered people is not directly considered in the system.

• 3 People can be in one of the three groups:

  • (a) Secondhand smokers or those at risk of others smoking, represented by P;

  • (b) People who have stopped smoking but are at risk due to their smoking habit before, represented by Q.

  • (c) People who are addicted to tobacco and now smoke it, represented by S.

• 4 There is a constant flow, $\alpha$, of healthy people who become at risk of active smokers to the state P.

• 5 People leave from the states P or S or Q under study due to factors such as living in a completely nonsmoking population or the death of the individual.

Description of variables and parameters

In the following table, we describe the variables and parameters to create the mathematical model that represents the dynamics of transmission of the habit of smoking (Table 1).

Table 1 Description of variables and parameters

The dynamic system

Fig. 1

Mathematical model of smoking

As we see in Fig. 1,

• 1 P increases by $\alpha$ per unit time

• 2 P decreases due to natural death ($\mu P$ people per unit time); the influence of factors that cause a person to move from the population at risk of being an active smoker to the population of healthy people ($\sigma P$ people per unit time); death as secondhand smoker ($\nu P$); and the impact of smoking ($\lambda PS$ people per unit time).

• 3 S increases by the impact of ($\lambda PS$ and $\delta QS$ people per unit time).

• 4 S decreases due to natural death ($\mu S$), death as a result of smoking tobacco ($\kappa S$), the influence of factors that cause a person to leave the population of active smokers and join the population of people who have stopped smoking.

• 5 Q increases due to the influence of factors that cause a person to leave the population of active smokers and join the population of people who have stopped smoking

• 6 Q decreases due to natural death ($\mu Q$), deaths as a result of smoker some days before ($\eta Q$) and the influence of factors that cause a person to move from the population to other areas and also medical treatment to the population of healthy people ($\xi Q$).

  We now describe the dynamic system as follows.

$$\begin{aligned} \dfrac{\hbox {d}P}{\hbox {d}t}=\, & {} \alpha -(\nu +\mu +\sigma +\lambda S)P \end{aligned}$$

(1)

$$\begin{aligned} \dfrac{\hbox {d}S}{\hbox {d}t}=\, & {} (\lambda P+\delta Q-(\mu +\kappa +\zeta )) S \end{aligned}$$

(2)

$$\begin{aligned} \dfrac{\hbox {d}Q}{t}=\, & {} \zeta S-(\delta S +\eta +\mu +\xi ) Q \end{aligned}$$

(3)

Positivity and boundedness of the solution

Theorem 1

If the initial population sizes of the model are positive, then the population sizes at any time are nonnegative. In other words, if $P(0)>0, S(0)>0$ and $Q(0)>0,$ then $P(t)>0, S(t)>0$ and $Q(t)>0$ for all t.

Proof

Equation (1) can be expressed as an inequality

$$\begin{aligned} \dfrac{\hbox {d}P}{P}\ge -(\nu +\mu +\sigma +\lambda S)\hbox {d}t. \end{aligned}$$

Integrating both sides from 0 to t, the solution is obtained as

$$\begin{aligned} P(t)\ge P(0)\hbox {e}^{-\int (\nu +\mu +\sigma +\lambda S)\hbox {d}t}. \end{aligned}$$

Since $P(0)>0, P(t)>0.$ In the same manner from Eqs. (2) and (3), we obtain

$$\begin{aligned} S(t)\ge S(0)\hbox {e}^{\int (\lambda P+\delta Q-(\mu +\kappa +\zeta ))\hbox {d}t}>0 \end{aligned}$$

and

$$\begin{aligned} Q(t)\ge Q(0)\hbox {e}^{-\int (\delta S +\eta +\mu +\xi )\hbox {d}t }>0 \end{aligned}$$

, respectively. The latter one is justified by $\zeta S(t)\ge 0$. $\square$

Theorem 2

All the solutions P(t), S(t) and Q(t) of system (1), (2) and (3) are bounded.

Proof

The population size N(t) is given by

$$\begin{aligned} N(t)=P(t)+S(t)+Q(t). \end{aligned}$$

Differentiating both sides with respect to t, we obtain

$$\begin{aligned} \dfrac{\hbox {d}N(t)}{\hbox {d}t}=\dfrac{\hbox {d}P(t)}{\hbox {d}t}+\dfrac{\hbox {d}S(t)}{\hbox {d}t}+\dfrac{\hbox {d}Q(t)}{\hbox {d}t}. \end{aligned}$$

(4)

Substituting (1), (2) and (3) in (4), we have

$$\begin{aligned} \dfrac{\hbox {d}N(t)}{\hbox {d}t}= & {} \alpha -(\nu +\mu +\sigma )P-(\mu +\kappa ) S-(\eta +\mu +\xi )Q \nonumber \\\le & {} \alpha -(P+S+Q)\mu \nonumber \\= & {} \alpha -N\mu . \end{aligned}$$

(5)

The solution of Eq. (5) is

$$\begin{aligned} N(t)\le \dfrac{\alpha }{\mu }+\left( N(0)-\dfrac{\alpha }{\mu }\right) \hbox {e}^{-t}. \end{aligned}$$

As $t\rightarrow \infty ,$ we have

$$\begin{aligned} N(t)\le \dfrac{\alpha }{\mu }. \end{aligned}$$

Therefore, the solution of the dynamic system (1), (2) and (3) is bounded. $\square$

So, let the total population size at time t be denoted by N(t). So, $N(t)=P(t)+Q(t)+S(t).$ We assume that N(t) is constant, and for the purpose of computation, we take P(t), Q(t) and S(t) as proportions of N(t) where $P(t)+Q(t)+S(t)=1.$ Thus, the region that we consider for our model is

$$\begin{aligned} \Omega =\{(P,Q,S)\in R^{3}:P+S+Q=1,\text {and}\; P,S,Q\ge 0\}. \end{aligned}$$

In this paper, P(t), Q(t) and S(t) are proportions of N(t).

Free equilibrium point

Smoking-free equilibrium points are steady-state solutions of a mathematical model indicating that there is no smoking (only presented secondhand smoker). Thus, setting $\dfrac{\hbox {d}P}{\hbox {d}t}=0,\dfrac{\hbox {d}S}{\hbox {d}t}=\,0$ and $\dfrac{\hbox {d}Q}{\hbox {d}t}=\,0,$ we obtain

$$\begin{aligned} \alpha -(\nu +\mu +\sigma +\lambda S)P= &\, 0 \end{aligned}$$

(6)

$$\begin{aligned} (\lambda P+\delta Q-(\mu +\kappa +\zeta )) S= &\, 0\end{aligned}$$

(7)

$$\begin{aligned} \zeta S-(\delta S +\eta +\mu +\xi ) Q= &\, 0 \end{aligned}$$

(8)

From (7), we get $S=0.$ Substituting this value in (6) and (8) and calculating P and Q, we have respectively $P=\dfrac{\alpha }{\nu +\mu +\sigma }$ and $Q=0.$ Thus, the free equilibrium point of the dynamic systems [(1), (2), (3)] is

$$\begin{aligned} (P,S,Q)=\left( \dfrac{\alpha }{\nu +\mu +\sigma },0,0\right) . \end{aligned}$$

Reproduction number

The basic reproduction number, $R_{0}$, of an infectious disease is the average number of secondary cases generated by a single primary case in a fully susceptible population [23]. $R_{0}$ is the most widely used epidemiological measurement of the transmission potential in a given population. Statistical estimation of $R_{0}$ has been performed for various infectious diseases [24, 25], aiming toward understanding the dynamics of transmission and evolution and designing effective public health intervention strategies. In particular, $R_{0}$ has been used for determining the minimum coverage of immunization, because the threshold condition to prevent a major epidemic in a randomly-mixing.

It is primarily used as a threshold parameter: If $R_{0} < 1,$ the disease will fade out of the population, but if $R_{0} > 1$, the disease will persist and become endemic to the population. Furthermore, the larger the magnitude of $R_{0}$, the faster the disease will spread and presumably the more difficult it will be to control.

Even though there are different methods in which $R_{0}$ can be calculated, we use the next-generation method. This is the most common method of calculating $R_{0}.$ In this method, we place appropriate terms from the infected class equations into the vectors F and V. Terms that describe appearances of new infections in each compartment belong in F, and other terms belong in V. The Jacobian matrices obtained by differentiating F and V with respect to the relevant subset of variables are computed and evaluated at a nontrivial disease-free equilibrium, resulting in the matrices DF and DV , respectively. So, the reproductive number $R_{0}$ is defined as the spectral radius of the matrix $DF(DV)^{-1},$ denoted by $\rho (DF(DV)^{-1}).$

Let $X=(P,S,Q).$ Then,

$$\begin{aligned} \dfrac{\hbox {d}X}{\hbox {d}t}=F(X)-V(X), \end{aligned}$$

where

$$\begin{aligned} F(X)=\left( \begin{array}{c} \alpha \\ \lambda PS +\delta QS\\ \zeta S \end{array}\right) \;\text { and}\; V(X)=\left( \begin{array}{c} (\nu +\mu +\sigma +\lambda S )P\\ (\mu +\kappa +\zeta ) S \\ (\delta S +\eta +\mu +\xi )Q \end{array}\right) \end{aligned}$$

Let $M_{0}=\left( \dfrac{\alpha }{\nu +\mu +\sigma },0,0\right) .$ Then,

$$\begin{aligned} DF(M_{0})=\left( \begin{array}{ccc} 0 &{}\quad 0 &{}\quad 0\\ 0 &{}\quad \dfrac{\lambda \alpha }{\mu +\nu +\sigma } &{}\quad 0\\ 0 &{}\quad \zeta &{}\quad 0 \end{array}\right) \end{aligned}$$

and

$$\begin{aligned} DV(M_{0})=\left( \begin{array}{ccc} \nu +\mu +\sigma &{}\quad \dfrac{\lambda \alpha }{\nu +\mu +\sigma } &{}\quad 0\\ 0 &{}\quad \mu +\kappa +\zeta &{}\quad 0\\ 0 &{}\quad 0 &{}\quad \eta +\mu +\xi \end{array}\right) , \end{aligned}$$

Here,

$$\begin{aligned} \hbox {det}(DV)=(\nu +\mu +\sigma )(\mu +\kappa +\zeta )(\eta +\mu +\xi ). \end{aligned}$$

Hence,

$$\begin{aligned} DV^{-1}=\dfrac{1}{\hbox {det}(DV)}\left( \begin{array}{ccc} A &{}\quad B &{}\quad 0\\ 0 &{}\quad C &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad D \end{array}\right) , \end{aligned}$$

where

$$\begin{aligned} A=\, & {} (\mu +\kappa +\zeta )(\eta +\mu +\xi ),\\ B=\, & {} -\dfrac{\lambda \alpha }{\nu +\mu +\sigma }(\eta +\mu +\xi )\\ C=\, & {} (\nu +\mu +\sigma )(\eta +\mu +\xi )\\ D=\, & {} (\nu +\mu +\sigma )(\mu +\kappa +\zeta ). \end{aligned}$$

And

$$\begin{aligned} DF(DV^{-1})= \,& {} \dfrac{1}{\hbox {det}(DV)}\left( \begin{array}{ccc} 0 &{}\quad 0 &{}\quad 0\\ 0 &{}\quad \dfrac{\lambda \alpha }{\mu +\nu +\sigma }C &{}\quad 0\\ 0 &{}\quad \zeta C &{}\quad 0\\ \end{array}\right) ,\nonumber \\ R_{0}=\, & {} \rho (DF(DV)^{-1})=\dfrac{\lambda \alpha C}{(\mu +\nu +\sigma )\hbox {det}(DV)}\nonumber \\=\, & {} \dfrac{\lambda \alpha (\nu +\mu +\sigma )(\eta +\mu +\xi )}{(\nu +\mu +\sigma )(\nu +\mu +\sigma )(\mu +\kappa +\zeta )(\eta +\mu +\xi )}\nonumber \\= \,& {} \dfrac{\lambda \alpha }{(\nu +\mu +\sigma )(\mu +\kappa +\zeta )} \end{aligned}$$

(9)

Here, $R_{0}$ represents the average number of smokers that a single smoker would produce during their period of infection in a population of exposed people.

Sensitivity analysis

Sensitivity analysis tells us how important each parameter is to interpret the model. Sensitivity allows us to measure the relative change in a variable when a parameter changes. The sensitivity of $R_{0}$ with respect to a parameter is the ratio of the relative change in the variable ($R_{0}$) to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives [26]. Thus, for a parameter p, we denote sensitivity of $R_{0}$ with respect to p by $\mathcal {S}^{R_{0}}_{p}$ and as in [27] defined by

$$\begin{aligned} \mathcal {S}^{R_{0}}_{p}=\dfrac{p}{R_{0}}\left[ \dfrac{\partial R_{0}}{\partial p}\right] . \end{aligned}$$

We calculate the sensitivity analysis of parameters involved in $R_{0}$ as follows.

$$\begin{aligned} \mathcal {S}^{R_{0}}_{\lambda }=\, & {} \dfrac{\lambda }{R_{0}}\left[ \dfrac{\alpha }{(\nu +\mu +\sigma )(\mu +\kappa +\zeta )}\right] >0 \end{aligned}$$

(10)

$$\begin{aligned} \mathcal {S}^{R_{0}}_{\alpha }=\, & {} \dfrac{\alpha }{R_{0}}\left[ \dfrac{\lambda }{(\nu +\mu +\sigma )(\mu +\kappa +\zeta )}\right] >0\end{aligned}$$

(11)

$$\begin{aligned} \mathcal {S}^{R_{0}}_{\nu }= \,& {} \dfrac{\nu }{R_{0}}\left[ \dfrac{-\lambda \alpha }{(\nu +\mu +\sigma )^{2}(\mu +\kappa +\zeta )}\right] <0\end{aligned}$$

(12)

$$\begin{aligned} \mathcal {S}^{R_{0}}_{\mu }=\, & {} \dfrac{\mu }{R_{0}}\left[ \dfrac{-(2\mu +\nu +\sigma +\kappa +\zeta )\lambda \alpha }{(\nu +\mu +\sigma )^{2}(\mu +\kappa +\zeta )^{2}}\right] <0\end{aligned}$$

(13)

$$\begin{aligned} \mathcal {S}^{R_{0}}_{\sigma }=\, & {} \dfrac{\sigma }{R_{0}}\left[ \dfrac{-\lambda \alpha }{(\nu +\mu +\sigma )^{2}(\mu +\kappa +\zeta )}\right] <0\end{aligned}$$

(14)

$$\begin{aligned} \mathcal {S}^{R_{0}}_{\kappa }=\, & {} \dfrac{\kappa }{R_{0}}\left[ \dfrac{-\lambda \alpha }{(\nu +\mu +\sigma )(\mu +\kappa +\zeta )^{2}}\right] <0\end{aligned}$$

(15)

$$\begin{aligned} \mathcal {S}^{R_{0}}_{\zeta }=\, & {} \dfrac{\zeta }{R_{0}}\left[ \dfrac{-\lambda \alpha }{(\nu +\mu +\sigma )(\mu +\kappa +\zeta )^{2}}\right] <0 \end{aligned}$$

(16)

We can see that $\mathcal {S}^{R_{0}}_{\lambda }$ and $\mathcal {S}^{R_{0}}_{\alpha }$ are positive. This indicates that $\lambda$ and $\alpha$ are directly proportional to $R_{0}.$ This means an increase (or decrease) in $\lambda$ and $\alpha$ will cause an increase (or decrease) in $R_{0}$ with the same proportion. So, the average number of secondary cases of infection increases (or decreases) in the community. It can also be seen that the remaining five parameters $\nu ,\mu ,\sigma ,\kappa$ and $\zeta$ are inversely proportional to $R_{0}.$ So, an increase (or decrease) in these parameters leads to minimizing (or maximizing) the endemic nature of smoking in the community. The relationship between $R_{0}$ and the above parameters can be described graphically as shown below (Fig. 2).

Fig. 2

Parameters versus reproductive number

Endemic equilibrium point

Endemic equilibrium is steady-state solutions when the smoking persists in the population. Assume that the smoking will persist in the population. Thus, in (7), if $S\ne 0,$ then

$$\begin{aligned} \lambda P+\delta Q=\mu +\kappa +\zeta . \end{aligned}$$

(17)

From (6) and (8), we get

$$\begin{aligned} P=\, & {} \dfrac{\alpha }{\nu +\mu +\sigma +\lambda S} \end{aligned}$$

(18)

$$\begin{aligned} Q=\, & {} \dfrac{\zeta S}{\eta +\mu +\xi +\delta S} \end{aligned}$$

(19)

Substituting (18) and (19) in (17), we obtain the quadratic equation (in a variable S)

$$\begin{aligned} (\lambda \delta \zeta -\lambda \delta L)S^{2}+(\delta ML(R_{0}-1)+\delta \zeta L-\lambda NM)S+LMN(R_{0}-1)=0, \end{aligned}$$

(20)

where $R_{0}$ is given as in (9), $L=\nu +\mu +\sigma ,M=\mu +\kappa +\zeta ,$ and $N=\eta +\mu +\xi .$ From Eq. (20), we can calculate

$$\begin{aligned} S=\dfrac{-V\pm \sqrt{V^{2}-4UW}}{2U}, \end{aligned}$$

(21)

where $U=\lambda \delta \zeta -\lambda \delta M, V=\delta ML(R_{0}-1)+\delta \zeta L-\lambda NM$ and $W=LMN(R_{0}-1).$ If $R_{0}>1,$, then

$$\begin{aligned} S=\dfrac{-V-\sqrt{V^{2}-4UW}}{2U} \end{aligned}$$

(22)

since $U<0.$ In this case, P and Q can be found by substituting Eq. (22) into Eqs. (18) and (19), respectively. If $R_{0}<1,$ then $-V>0$ and $W<0$, and thus either S is negative real number or complex number. In this case, S does not exist and hence P and Q. If $R_{0}=1,$ then $V=\delta \zeta L-\lambda N M$ and $W=0.$ Consequently,

$$\begin{aligned} S=\dfrac{-V-|V|}{2U}=\left\{ \begin{array}{ll} \text {negative,} &{} \text {if}\;V>0\\ \text {zero,} &{} \text {if}\;V\le 0 \end{array}\right. \end{aligned}$$

We write the above discussion as a theorem as follows.

Theorem 3

The endemic equilibrium point of the dynamic system (1), (2) and (3) exists if $R_{0}>1$ or $R_{0}=1$ and $\delta \zeta L>\lambda N M$ and does not exist if $R_{0}<1.$ For $R_{0}>1$ or $R_{0}=1$ and $\delta \zeta L>\lambda N M,$ the endemic equilibrium point is the ordered triple (P, S, Q), where S, P and Q are given as in (22), (18) and (19), respectively.

Stability

The physical stability of the equilibrium solution of the dynamic system (6), (7) and (8) is related to its eigenvalue.

Definition 4

The equilibrium solution of the dynamic system is stable if all the eigenvalues of the Jacobian matrix of the dynamic system evaluated at the equilibrium solution have negative real part, otherwise the solution is unstable.

The Jacobian matrix of the dynamic system is

$$\begin{aligned} \left( \begin{array}{ccc} -(\nu +\mu +\sigma +\lambda S) &{} -\lambda P &{} 0\\ \lambda S &{} \lambda P+\delta Q-(\mu +\kappa +\zeta ) &{} \delta S\\ 0 &{} \zeta -\delta Q &{} -(\delta S +\eta +\mu +\xi ) \end{array}\right) \end{aligned}$$

(23)

Free equilibrium point

At the free equilibrium point, the Jacobian matrix (23) becomes

$$\begin{aligned} \left( \begin{array}{ccc} -(\nu +\mu +\sigma ) &{} \dfrac{-\lambda \alpha }{\nu +\mu +\sigma } &{} 0\\ 0 &{} (\mu +\kappa +\zeta )(R_{0}-1) &{} 0\\ 0 &{} \zeta &{} -(\eta +\mu +\xi ) \end{array}\right) . \end{aligned}$$

The eigenvalues of this matrix are

$$\begin{aligned} \lambda _{1}=-(\nu +\mu +\sigma ),\lambda _{2}=(\mu +\kappa +\zeta )(R_{0}-1)\;\text {and}\; \lambda _{3}=-(\eta +\mu +\xi ). \end{aligned}$$

We see that $\lambda _{1}$ and $\lambda _{3}$ are negative. If $R_{0}<1,$ then $\lambda _{2}<0$ and hence the free equilibrium point is stable. If $R_{0}>1,$ then $\lambda _{2}>0$ and hence the free equilibrium point is unstable. From this discussion, we now have the following theorem.

Theorem 5

If $R_{0}<1,$ then the free equilibrium point is stable. If $R_{0}>1,$ then the free equilibrium point is unstable.

Endemic equilibrium point

The Jacobian matrix (23) of the dynamic system (6), (7) and (8) at the endemic equilibrium point is

$$\begin{aligned} \left( \begin{array}{ccc} -(L +\lambda S^{*}) &{} -\lambda P^{*} &{} 0\\ \lambda S^{*} &{} 0 &{} \delta S^{*}\\ 0 &{} \zeta -\delta Q^{*} &{} -(\delta S^{*} +N) \end{array}\right) , \end{aligned}$$

where $(P^{*},S^{*},Q^{*})=(P,S,Q),$ here P, S and Q are given as in (18), (19) and (22). Its characteristic polynomial is

$$\begin{aligned} \mathcal {P}(\gamma )=\gamma ^{3}+a_{1}\gamma ^{2}+a_{2}\gamma +a_{3}, \end{aligned}$$

(24)

where

$$\begin{aligned} a_{1}=\, & {} L+\lambda S^{*}+\delta S^{*}+N \end{aligned}$$

(25)

$$\begin{aligned} a_{2}=\, & {} \delta LS^{*}+LN+\delta \lambda (S^{*})^{2}+\lambda NS^{*}+\delta ^{2}S^{*}Q^{*}\nonumber \\&-\,\zeta \delta S^{*}+\lambda ^{2}S^{*}P^{*} \end{aligned}$$

(26)

$$\begin{aligned} a_{3}=\, & {} \delta ^{2}LS^{*}Q^{*}+\lambda \delta ^{2}(S^{*})^{2}Q^{*}- L\delta \zeta S^{*}-\lambda \zeta \delta (S^{*})^{2}+\lambda ^{2}\delta P^{*}(S^{*})^{2}\nonumber \\&+\lambda ^{2}NS^{*}P^{*} \end{aligned}$$

(27)

Here, our objective is to show all the roots of the characteristic polynomial (24) have negative real part. However, it is difficult to show this due to the complexity of $a_{1},a_{2}$ and $a_{3}.$ Later, we will illustrate the stability of the endemic equilibrium solution using numerical simulations, by imposing the value of parameters or variables. It is easy to see that if $a_{3}>0$ and $a_{1}a_{2}-a_{3}>0,$ then all the roots of (24) have negative real part. In this case, the endemic equilibrium point is stable by Routh–Hurwitz criteria. By this criteria, if $a_{3}<0$ or $a_{1}a_{2}-a_{3}<0,$ then the endemic equilibrium point is unstable.

Numerical simulation

The following initial values (IV) for P, S and Q have been taken at $t=0$.

	IV1	IV2	IV3	IV4
P(0)	0.75	0.65	0.55	0.45
S(0)	0.15	0.20	0.25	0.30
Q(0)	0.10	0.15	0.20	0.25

For numerical simulation, we use Runge–Kutta 4–5 methods and MATLAB 2018 software. It will be seen separately for free and endemic equilibrium point.

Free equilibrium point

The values of parameters used for simulation are $\alpha =0.5, \delta =0.5, \nu =0.02, \mu =0.0135, \sigma =0.2, \kappa =0.02, \zeta =0.3, \eta =0.02, \xi =0.1$, and the remaining parameter $\lambda$ can be chosen by considering $R_{0}.$

	$\lambda$
$R_{0}<1$	0.035
$R_{0}>1$	0.165

Since the free equilibrium point does not depend on $\lambda ,$ in both cases, we have

$$\begin{aligned} (P,S,Q)=(2.1413,0,0). \end{aligned}$$

We will construct the graphs in each case for $t\in [0,400]$ as follows.

• (i) $\lambda =0.035$

  As it can be seen in Fig. 3, for $R_{0}<1$ we conclude that the graphs of P, S and Q for different initial values approach $P=2.1413,S=0$ and $Q=0$, respectively, as $t\rightarrow \infty$ respectively. We recall that the point $(P,S,Q)=(2.1413,0,0)$ is smoking-free equilibrium point of the model. Hence, the smoking-free equilibrium point is stable.

• (ii) $\lambda =0.165$

  As it is described in Fig. 4, for $R_{0}>1$ we conclude that the graphs of P, S and Q for different initial values do not approach $P=2.1413,S=0$ and $Q=0$, respectively, as $t\rightarrow \infty$. We again recall that the point $(P,S,Q)=(2.1413,0,0)$ is smoking-free equilibrium point of the model. Thus, smoking-free equilibrium point is unstable.

Fig. 3

$R_{0}=0.2242<1$

Fig. 4

$R_{0}=1.0594>1$

Endemic equilibrium point

We take $\alpha =0.3, \nu =0.01, \mu =0.0035, \sigma =0.2, \lambda =0.2, \delta =0.05,\kappa =0.01,\zeta =0.2, \eta =0.02$ and $\xi =0.1.$ In this case, the endemic equilibrium point is

$$\begin{aligned} (P,S,Q)=(0.8846,0.7085,0.8916). \end{aligned}$$

We will construct the graphs for $t\in [0,400]$ as follows. In Fig. 5, we see that for $R_{0}>1$ we conclude that the graphs of P, S and Q for different initial values approach $P=0.8846,S=0.7085$ and $Q=0.8916$ as $t\rightarrow \infty$, respectively. Therefore, the endemic equilibrium point is stable.

Fig. 5

$R_{0}=1.3163>1$

We next consider $\alpha =0.3, \nu =0.01, \mu =0.035, \sigma =0.2, \lambda =0.2001,\delta =0.2,\kappa =0.01,\zeta =0.2, \eta =0.02$ and $\xi =0.1.$ Here, the endemic equilibrium point is

$$\begin{aligned} (P,S,Q)=(0.6125,1.2232,0.6121). \end{aligned}$$

We will construct the graphs for $t\in [0,600]$ as follows. In Fig. 6, we notice that for $R_{0}=1$ and $\delta \zeta L>\lambda NM$ we conclude that the graphs of P, S and Q for different initial values approach $P=0.6125,S=1.2232$ and $Q=0.6121$ as $t\rightarrow \infty$, respectively. So, here the endemic equilibrium point is stable.

Fig. 6

$R_{0}=1$ and $\delta \zeta L>\lambda NM$

Conclusion

• 1 The average number of new smokers from a single addicted smoker is given by

  $$\begin{aligned} \dfrac{\alpha \lambda }{(\nu +\mu +\sigma )(\mu +\kappa +\zeta )}. \end{aligned}$$

  Secondhand smoker persists in the population if $R_{0}>1$ or $R_{0}=1$ and $\delta \zeta L>\lambda N M.$ Secondhand smoker will be disappeared over time if $R_{0}<1.$

• 2 An increase (or decrease) in $\lambda$ leads an increase (or decrease) secondhand smoker in the community.

• 3 An increase (or decrease) in the parameters $\nu ,\mu ,\sigma ,\kappa$ and $\zeta$ leads to minimizing (or maximizing) the endemic nature of smoking in the community. So secondhand smoker will be minimized or maximized.

• 4 If the average number of smokers that a single addicted smoker produces is less than one, the population of secondhand smokers disappears over time.

• 5 If the average number of smokers that a single addicted smoker produces is greater than one, the population of secondhand smokers persists.