Curtailing smoking dynamics: A mathematical modeling approach
Available online at www.sciencedirect.com
Applied Mathematics and Computation 195 (2008) 475–499 www.elsevier.com/locate/amc
Curtailing smoking dynamics: A mathematical modeling approach O. Sharomi, A.B. Gumel
*
Department of Mathematics, University of Manitoba, Winnipeg, Manitoba, Canada R3T 2N2
Abstract This paper provides a rigorous mathematical study for assessing the dynamics of smoking and its public health impact in a community. A basic mathematical model, which is a slight refinement of the model presented in [F. Brauer, C. CastilloChavez. Mathematical Models in Population Biology and Epidemiology. Text in Applied Mathematics. Springer, 2000; G.C. Castillo, S.G. Jordan, A.H. Rodriguez. Mathematical models for the dynamics of tobacco use, recovery and relapse. Technical Report Series, BU-1505-M. Department of Biometrics, Cornell University. 2000], is designed first of all. It is based on subdividing the total population in the community into non-smokers, smokers and those smokers who quit smoking either temporarily or permanently. The theoretical analysis of the basic model reveals that the associated smoking-free equilibrium is globally-asymptotically stable whenever a certain threshold, known as the smokers-generation number, is less than unity, and unstable if this threshold is greater than unity. The public health implication of this result is that the number of smokers in the community will be effectively controlled (or eliminated) at steady-state if the threshold is made to be less than unity. Such a control is not feasible if the threshold exceeds unity (a global stability result for the smoking-present equilibrium is provided for a special case). The basic model is extended to account for variability in smoking frequency, by introducing two classes of mild and chain smokers as well as the development and the public health impact of smoking-related illnesses. The analysis and simulations of the extended model, using an arbitrary but reasonable set of parameter values, reveal that the number of smokers in the community will be significantly reduced (or eliminated) if chain smokers do not remain as chain smokers for longer than 1.5 years before reverting to the mild smoking class, regardless of the time spent by mild smokers in their (mild smoking) class. Similarly, if mild smokers practice their mild smoking habit for less than 1.5 years, the number of smokers in the community will be effectively controlled irrespective of the dynamics in the chain smoking class. 2007 Elsevier Inc. All rights reserved. Keywords: Tobacco smoking; Smoking-related illnesses; Equilibrium; Stability; Mild/chain smokers
1. Introduction Since the advent of tobacco in 6000 BC [29], tobacco smoking remains a major public health menace globally. It inflicts significant mortality, morbidity and continues to cost billions of dollars in health-care (see [27] *
Corresponding author. E-mail address: [email protected] (A.B. Gumel).
0096-3003/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.05.012
476
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
and the references therein). Recent figures from the World Health Organization (WHO) estimate that 100 million people died of tobacco-associated diseases (such as cancer, chronic lung disease, cardiovascular disease and stroke) last century [31]. Although smoking is widely acknowledged as the most important preventable cause of death, 4 million people die of smoking-related illnesses globally every year, and the number of new smokers continues to rise [30]. Owing to these facts, and the astronomical public health burden associated with smoking, it is of public health significance to qualitatively study the dynamics of smoking in a community, aimed at determining realistic methods for combatting this habit. This would allow for the theoretical assessment of tobacco smoking control and protection from second-hand smoking etc. [25]. Although mathematical modelling has been extensively used to address questions of public health importance, dating back to the pioneering works of Bernoulli (on modelling the dynamics of smallpox) in 1760 [3], Kermack and McKendrick [19–21] and those reported in more recent literature (such as in [1,2,4,6–9,11– 15,18,24,26,28] and the references therein), not much has been done in terms of the mathematical modeling of human social behaviour. In particular, with the exception of the basic model in [4,10], the authors are not aware of any mathematical study (talkless of the rigorous qualitative analyses of the models) for fully examining, and assessing, the impact of smoking in a human population. The aim of this paper is to provide a qualitative study of the dynamics of smoking, vis-a-vis its impact on public health. A basic compartmental model, which subdivides a given population into a number of mutuallyexclusive sub-populations (depending on their smoking status), is designed and qualitatively analysed. This model is then extended to incorporate variability in smoking habits (namely mild or chain smoking) as well as the impact of smoking-related illnesses in the community. 2. Model formulation of basic model The basic model is based on monitoring the dynamics of the sub-populations of potential smokers (nonsmokers); P(t), smokers; S(t), smokers who temporarily quit smoking; Qt(t), and smokers who permanently quit smoking; Qp(t). Thus, the total population at time t is given by N(t) = P(t) + S(t) + Qt(t) + Qp(t). Non-smokers (potential smokers) (P(t)) The class of non-smokers (potential smokers) is increased by the recruitment of individuals into the nonsmoking class (at a rate lN). It is assumed that non-smokers can acquire smoking habits (and become smokers) via effective ‘‘contacts’’ with smokers (at a rate b = cq, where c is the average number of contacts per unit time; and q is the probability of becoming a smoker for a member of the P class following contact with a smoker). In other words, it is assumed that the acquisition of smoking habit is analogous to acquiring disease infection; the mode ‘‘contact’’ a potential smoker makes, the higher the likelihood of such an individual acquiring smoking habit. Although not considered in this study, other factors such as tobacco advertising and mass consumption of tobacco in the community can influence this process (and, hence, should be included in the function governing the generation of new smokers in the community). Finally, potential smokers suffer natural death (at a rate l), so that the rate of change of the population of potential smokers is given by dP bPS ¼ lN lP : dt N
ð1Þ
Smokers (S(t)) The population of smokers is increased when non-smokers begin to smoke (regardless of their frequency of smoking per unit time) and when smokers who temporarily quit smoking revert to smoking (at a rate a). The population is decreased when smokers temporarily quit smoking (at a rate c) and natural death (at the rate l), giving dS bPS ¼ þ aQt ðl þ cÞS: dt N
ð2Þ
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
477
Smokers who temporary quit smoking (Qt(t)) This is generated by a fraction, 1 r, of smokers who temporarily quit smoking (at the rate c). It is decreased by natural death (at the rate l) and reversion to smoking (at the rate a). Thus, dQt ¼ cð1 rÞS ðl þ aÞQt : dt
ð3Þ
Smokers who permanently quit smoking (Qp(t)) This is generated by the remaining fraction, r, of smokers who permanently quit smoking (at the rate c), and decreased by natural death (at the rate l), so that dQp ¼ crS lQp : dt
ð4Þ
Thus, the basic model for smoking dynamics is given by Eqs. (1)–(4), (see flowchart diagram in Fig. 1a). It is an extension of the model in [4,10], by incorporating a class (and dynamics) of smokers who temporarily quit smoking (Qt(t)). Note that this model does not include smoking-related mortality.
2.1. Basic properties of the model Adding Eqs. (1)–(4) gives dN ¼ 0: dt
ð5Þ
Thus, the total population, N, is constant; and equals N*. Since the models (1)–(4) monitors human population, it is plausible to assume that all its state variables and parameters are non-negative for all t P 0. Further, it can be shown that the region D ¼ fðP ; S; Qt ; Qp Þ 2 R4þ : P þ S þ Qt þ Qp ¼ N g; is positively-invariant. Thus, every solution of the models (1)–(4), with initial conditions in D remains there for t > 0. Therefore, the x limit sets of solutions of (1)–(4), in D are contained in D. Furthermore, in D, the
Fig. 1a. Flowchart of the basic models (1)–(4).
478
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
Fig. 1b. Flowchart of the extended model (13).
usual existence, uniqueness and continuation results hold for the system, so that the system (1)–(4), is wellposed mathematically and epidemiologically [18]. Hence, it is sufficient to consider the dynamics of the flow generated by models (1)–(4), in D. 2.2. Stability of smoking-free equilibrium (SFE) 2.2.1. Local stability The models (1)–(4) has an SFE, obtained by setting the right-hand sides of the equations in (1)–(4) to zero, given by E0 : ðP ; S ; Qt ; Qp Þ ¼ ðN ; 0; 0; 0Þ;
with N ¼ P þ S þ Qt þ Qp ¼ P :
ð6Þ
The linear stability of E0 can be established using the next generation method [11,28] (see Appendix A for a detailed description of the next generation method). Note, first of all, that system (1)–(4) satisfies the axioms A1–A5 in Appendix A. Thus, using the notation in [28], the next generation method entails writing the right hand sides of the equations of the rate of change of the smoking-present variables (dS/dt and dQt/dt) in terms of two matrices F and V, where F is a matrix consisting of the smokers-generation terms (i.e., all terms with b in them in the equations for (dS/dt and dQt/dt) and V is an M-matrix consisting of the remaining transition terms in the two equations (it should be recalled that a matrix A is an M-matrix if and only if every off-diagonal entry of A is non-positive and the diagonal entries are all non-negative). That is, Eqs. (2) and (3) are re-written as (note that P* = N*) ! dSðtÞ b 0 K 1 a SðtÞ dt ¼ ; dQt ðtÞ K 2 K 3 0 0 Qt ðtÞ dt SðtÞ ¼ ðF V Þ ; ð7Þ Qt ðtÞ
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
479
with, K 1 ¼ l þ c;
K 2 ¼ cð1 rÞ;
K 3 ¼ l þ a:
That is, for the models (1)–(4), the next generation matrices F and V are given by b 0 K 1 a F ¼ and V ¼ : 0 0 K 2 K 3 Using the next generation method, the local stability of the SFE, given by (6), is based on whether or not q(FV1) < 1, where q is the spectral radius. If q(FV1) < 1, then all eigenvalues of the linearized version of the models (1)–(4), at the SFE, given by ! dSðtÞ SðtÞ dt ¼ ðF V Þ ; ð8Þ dQt ðtÞ Qt ðtÞ dt
have negative real parts, so that the SFE is locally-asymptotically stable (LAS). For q(FV1) > 1, at least one of the eigenvalues of the linearization has a positive real part; thus, the SFE is unstable in this case. Letting Rs ¼ qðFV 1 Þ, it is easy to show that Rs ¼
bK 3 bðl þ aÞ ¼ : K 1 K 3 aK 2 lðl þ aÞ þ cðar þ lÞ
ð9Þ
Consequently, using Theorem 2 of [28] (see also Theorem 8 of Appendix A), the following result is established. Lemma 1. For the system (1)–(4), the SFE, E0 , is LAS if Rs < 1, and unstable if Rs > 1. The quantity Rs is the smokers generation number. It measures the average number of new smokers generated by a single smoker in a population of potential smokers. The local stability result in Lemma 1 implies that for Rs < 1, the total number of smokers in the population can be reduced to zero if the initial sizes of the sub-populations of the model are in the basin of attraction of E0 . That is, a small influx of smokers into the community will not generate large number of smokers if Rs < 1 (but will do so if Rs > 1). To ensure that the effective control (or elimination) of the number of smokers in the community at steady-state is independent of the initial sizes of the sub-populations of the model, it is imperative to show that the SFE is globally-asymptotically stable (GAS). This is done below. 2.2.2. Global stability of SFE We claim the following result: Theorem 1. The SFE is GAS in D if Rs < 1. Proof. First, it should be noted that P 6 N in D for all t P 0. Consider the following Lyapunov function: F ¼ ðl þ aÞS þ aQt ; with Lyapunov derivative (where a dot represents differentiation with respect to t) _ ¼ ðl þ aÞS_ þ aQ_ t F bPS ¼ ðl þ aÞ ðl þ cÞS þ aQt þ a½cð1 rÞS ðl þ aÞQt N bPSðl þ aÞ ðl þ aÞðl þ cÞS þ ðl þ aÞaQt þ acð1 rÞS aðl þ aÞQt N bPSðl þ aÞ ðl þ aÞðl þ cÞS þ acð1 rÞS ¼ N bP ðl þ aÞ ½ðl þ aÞðl þ cÞ acð1 rÞ S ¼ N ¼
480
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
6 fbðl þ aÞ ½lðl þ aÞ þ cðar þ lÞgS for P 6 N bðl þ aÞ 1 S ¼ ½lðl þ aÞ þ cðar þ lÞ lðl þ aÞ þ cðar þ lÞ ¼ ½lðl þ aÞ þ cðar þ lÞðRs 1ÞS < 0
for Rs < 1:
_ < 0 for Rs < 1; with F _ ¼ 0 only if Since all the model parameters are non-negative, it follows that F S = 0. Hence, F is a Lyapunov function on D. Therefore, by the LaSalle’s Invariance Principle [16], every solution to the equations of the models (1)–(4), with initial conditions in D, approaches E0 as t ! 1. h It is worth stating that this theorem can also be proved using other approaches such as the use of comparison theorem (see Appendix B). 2.3. Smoking-present equilibrium (SPE) 2.3.1. Existence of smoking-present equilibrium The existence of positive smoking-present equilibria (steady-states with the smoking-present components, S** and Q ¼ bS , where ** represents the comt , positive) of the models (1)–(4), is explored as follows. Let G N ponent of the SPE at steady-state (with N ¼ P þ S þ Qt þ Qp ). Then, the state variables of the models (1)–(4), can be expressed at an arbitrary SPE, denoted by E1 ¼ ðP ; S ; Q t ; QP Þ, as: P ¼ Q t
lN ; G þ l
S ¼
K 3 lG N ; ðG þ lÞ½lðl þ aÞ þ cðar þ lÞ
K 2 lG N ; ¼ ðG þ lÞ½lðl þ aÞ þ cðar þ lÞ
Q p
crK 3 G N : ¼ ðG þ lÞ½lðl þ aÞ þ cðar þ lÞ
ð10Þ
Substituting the equations in (10) into the expression for G** above, and simplifying, gives the following fixed-point problem: G ¼
bG K 3 l : l½lðl þ aÞ þ cðar þ lÞ þ G ðlK 3 þ lK 2 þ crK 3 Þ
ð11Þ
The positive smoking-present equilibria of (1)–(4) can then be obtained by solving for G** in (11) and substituting the result into (10). Clearly, G** = 0 is a fixed-point of (11), which corresponds to the SFE ðE0 Þ. Since we are only interested in the possible smoking-present equilibria of (1)–(4), we consider only the case where G** 5 0. For G** 5 0, Eq. (11) can be simplified to: a0 G a1 ¼ 0; K 3 ðlþcrÞþlK 2 lðlþaÞþcðarþlÞ
ð12Þ
¼ 1 and a1 ¼ lðRs 1Þ: where a0 ¼ It follows from (12) that G ¼ aa10 . Further, a1 > 0 for Rs > 1. Thus, the linear Eq. (12) has a unique positive solution when Rs > 1 (i.e., G** > 0 iff Rs > 1). It should be noted that if Rs < 1, then a1 < 0. Hence, for Rs < 1, the quantity G** < 0 (which is biologically meaningless). Consequently, the SPE ðE1 Þ does not exist in this case. Furthermore, for Rs ¼ 1, the coefficient a1 = 0; making G** = 0 (corresponding to the SFE). These results are summarized below. Theorem 2. The model (1)–(4), has a unique SPE, E1 , if Rs > 1, and no SPE if Rs 6 1. 2.3.2. Local stability of smoking-present equilibrium The local stability of the unique SPE is now explored. We claim the following: Theorem 3. The unique SPE of the model (1)–(4), is LAS if Rs > 1.
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
481
Proof. Evaluating the Jacobian of (1)–(4), at E1 gives 2 3 G 0 0 G l 1 6 G G a 0 7 6 7 1 K1 J ðE1 Þ ¼ 6 7; 4 K 3 0 5 0 K2 0 K4 0 l where, bP and K 4 ¼ cr; N from which it follows that the eigenvalues of J ðE1 Þ are the roots of the polynomial G1 ¼
b0 k4 þ b1 k3 þ b2 k2 þ b3 k þ b4 ¼ 0; with, b0 ¼ 1; aK 2 ; b1 ¼ K 3 þ 2l þ G þ K3 aK 2 b2 ¼ K 3 G þ 2l K 3 þ G þ þ cG þ l2 ; K3 aK 2 þ 2lK 3 G þ l2 K 3 þ cð2l þ raÞG ; b3 ¼ l2 G þ K3 b4 ¼ l2 K 3 G þ lcðl þ raÞG : Thus, the coefficients bi ði ¼ 0; 1; 2; 3; 4Þ are all positive if Rs > 1 (note that G** > 0 iff Rs > 1). Finally, it can be shown that b1 b2 b3 b23 b21 b4 ¼
ðf1 þ f2 þ f3 Þf4 > 0; K 33
where, f1 ¼ la2 K 22 þ l2 aK 2 K 3 þ 2G laK 2 K 3 þ 2laK 2 K 23 þ 2lK 33 G ; f2 ¼ lK 43 þ G l2 K 23 þ lK 23 G 2 þ l2 K 33 þ aK 2 K 23 G þ K 33 G2 ; f3 ¼ K 43 G þ acK 2 K 3 G þ cK 23 G fG þ ½l þ að1 rÞg; f4 ¼ craG K 3 þ 2l2 aK 2 þ 2l2 K 3 G þ 2l3 K 3 þ 2lK 23 G þ 2l2 K 23 þ 2clG K 3 : Hence, by the Routh–Hurwitz criterion, all eigenvalues of J ðE1 Þ have negative real parts whenever Rs > 1. Thus, E1 is LAS whenever it exists. h It is worth mentioning that for the special case r = 0, Qp ! 0 as t ! 1, giving a 3D limiting model (in P, S and Qt). This (reduced) model can further be simplified into a 2D model using the relation N = S + P + Qt; and the corresponding SPE can be shown to be globally-asymptotically stable (using, for instance, the Dulac criterion). For the full models, (1)–(4), although no proof for global stability is presented, numerous simulations of the model suggest that E1 is globally-asymptotically stable whenever Rs > 1 (as illustrated in Fig. 2b, d and f). Consequently, we offer the following conjecture. Conjecture 1. The unique smoking-present equilibrium of the models (1)–(4), is GAS in D n D0 , with D0 ¼ fðP ; S; Qt ; Qp Þ 2 D : S ¼ Qt ¼ 0g, whenever Rs > 1. 2.4. Numerical simulations of the basic model To illustrate the various theoretical results presented thus far, the basic models (1)–(4), is simulated using the parameter values/ranges in Table 1. These parameters are chosen arbitrarily, albeit plausibly, owing to the absence of available reliable data in the literature.
482
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499 1200
1000 900
1000
800
800
Smokers S(t)
Smokers S(t)
700 600 500 400
600
400
300 200
200
100 0 0
10
20
30
40
50
0 0
60
50
100
500
500
450
450
400
400
350 300 250 200 150
350 300 250 200 150
100
100
50
50
0 0
5
10
15
20
25
0 0
30
10
20
30
40
50
60
70
80
Time (years)
1200
1200
1000
1000
Total smokers s(t) + Qt(t)
Total smokers s(t) + Qt(t)
Time (years)
800
600
400
800
600
400
200
200
0 0
150
Time (years)
Temporary smokers Qt(t)
Temporary Quitters Qt(t)
Time (years)
10
20
30
Time (years)
40
50
60
0 0
20
40
60
80
100
120
140
Time (years)
Fig. 2. Time series plot of the basic models (1)–(4), with different initial conditions. Parameter as in Table 1. (a) S(t) for Rs < 1; (b) S(t) for Rs > 1; (c) Qt(t) for Rs < 1; (d) Qt(t) for Rs > 1; (e) S(t) + Qt(t) for Rs < 1; (f) S(t) + Qt(t) for Rs > 1.
Table 2 depicts the variables of the model at steady-state as a function of b and Rs . It is clear from this table that the number of smokers at steady-state increases with increasing Rs > 1. It also shows a marked decrease in the number of potential smokers as Rs increases. Further, in line with Theorems 1 and 3 and Conjecture 1,
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
483
Table 1 Description of parameters for the basic models (1)–(4) Parameter
Description
Estimated value/range
Nominal value
1/l a
Average duration of smoking activity Rate at which temporary quitters Revert back to smoking Rate of quitting smoking Fraction of smokers who permanently quit Contact rate
20 years
20
c r b
a 2 [0, 1) c 2 [0, 1) r 2 ð0; 1Þ 0.3 per year
3 0.3 0.5 0.3
Table 2 Effect of b on Rs and P ; S; Qt ; Qp at steady-state using the basic models (1)–(4) b
Rs
P
S
Qt
Qp
Stable equilibrium
0.003 0.03 0.02 0.01 0.1 0.2 0.3 0.5 1 2 3 5
0.0118 0.1184 0.0789 0.0394 0.3948 0.7896 1.1844 1.9741 3.9482 7.8964 11.8446 19.7411
1600 1600 1600 1600 1600 1600 1351 811 405 203 135 81
0 0 0 0 0 0 49 156 236 276 289 300
0 0 0 0 0 0 3 10 16 18 19 20
0 0 0 0 0 0 197 623 943 1103 1157 1199
E0 E0 E0 E0 E0 E0 E1 E1 E1 E1 E1 E1
this table shows that E0 is always stable whenever Rs 6 1 and E1 is stable if Rs > 1 (see also Fig. 2 for simulations using various initial conditions). This table further shows a marked decrease in the population of potential smokers as Rs increases. The combined effect of the rate of reversion to smoking (a) and quitting smoking temporarily (c) on Rs is investigated using various values of r. The results, depicted in Fig. 3, show a decrease in Rs as r increases. In particular, for r = 0.5 (Fig. 3), almost all of the Rs contours are less than unity, signifying that this value of r (and other fixed parameters) will lead to the elimination of smoking from the community (in line with Theorem 1). It is worth noting that this result is further supported by carrying out a sensitivity analysis on Rs (i.e., differentiating Rs partially with respect to r), giving oRs barðl þ aÞ ¼ < 0; 2 or ½lðl þ aÞ þ cðl þ raÞ from which it follows that Rs is a decreasing function of r. Fig. 3 further shows that Rs increases with increasing a. Here, too, sensitivity analysis on Rs shows oRs blcð1 rÞ ¼ > 0: 2 oa ½lðl þ aÞ þ cðl þ raÞ It is further evident from Fig. 3 that increasing the quitting rate, c, decreases Rs ; which is in line with the calculation below: oRs bðl þ aÞðl þ arÞ ¼ < 0: oc ½lðl þ aÞ þ cðl þ raÞ2 3. Extension of the basic model Although the basic models (1)–(4), allows for the assessment of the dynamics of smoking in a closed population, it does not incorporate some of the other important aspects (and consequences) associated with
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
2
2 3
Number of days before quitting temporarily(1/γ)
4
1.8 1.6 3
4
1
2
1.4 1.2 3
1
2
0.8 1
2
0.6 1
0.4 1
0.2 1000
2000
3000
4000
5000
6000
7000
8000
2.5
1.8 1.6 1.4
1.5
2
1.2
1 5 0.
1.5
1
1
0.8 1
0.5
0.6 0.5
0.4 0.5
0.2
9000 10000
1000
Number of days before returning back to smoking(1/α)
1.6
1.8
1.4
1.2
1.4
1.2 1
1.4 1.2
0.8
1
1.6
0.8
0.6
1.2 1
0.8
0.6
1
0.4
0.8
0.8
0.6
0.4
0.6
0.6
0.4
0.4
0.2
0.4 0.2
0.2
0.2 1000
2000
3000
4000
3000
4000
5000
6000
7000
8000
9000 10000
2
1.4
1.8 1.6
2000
Number of days before returning back to smoking(1/α)
Number of days before quitting temporarily(1/γ)
Number of days before quitting temporarily(1/γ)
2
1.5
2
1
5
Number of days before quitting temporarily(1/γ)
2
1
484
5000
6000
7000
8000
1.2
0.8
1
1.8 1.6
1.2
1
0.8 0.6
1.4 1
0.8
1.2
0.6 0.8
1 0.8
0.6 0.4
0.6
0.4
0.6
0.4
0.2
0.4
0.2 0.2
0.2
9000 10000
1000
Number of days before returning back to smoking(1/α)
2000
3000
4000
5000
6000
7000
8000
9000 10000
Number of days before returning back to smoking(1/α)
2
Number of days before quitting temporarily(1/γ)
1 0.8
1.8 1 0.8
1.6
0.6
1.4
0.8
0.6
1.2
0.6
1 0.4
0.8
0.4 0.4
0.6 0.2
0.4
0.2 0.2
0.2 1000
2000
3000
4000
5000
6000
7000
8000
9000 10000
Number of days before returning back to smoking(1/α)
Fig. 3. Contour plots of the basic models (1)–(4), depicting Rs as a function of 1c and 1a. Parameters as in Table 1 with various values of r. (a) r = 0.1; (b) r = 0.1; (c) r = 0.1; (d) r = 0.1; (e) r = 0.1.
tobacco smoking. One of such is the degree, or frequency, of smoking, where some smokers do so in a mild manner (categorized as mild smokers) while others may smoke more frequently (defined as chain smokers) per
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
485
unit time. Further, the transition dynamics (back-and-forth movement) between these subclasses, of varying smoking frequencies, also needs to be incorporated into the model. Another key aspect lacking in the basic model is the public health consequence of smoking. In other words, since millions of people die of smoking-related illnesses every year [17], it is necessary to include such fact in a population model of smoking dynamics. To account for the aforementioned drawbacks, the basic model is now extended to include two classes of smokers – a mild smoking class (S1(t)) and a chain smoking class (S2(t)) – and a class of smoking-related illnesses (C(t)). The extended model is given by: dP bP ðS 1 þ wS 2 Þ ¼P lP ; dt N dS 1 ð1 hÞbP ðS 1 þ wS 2 Þ ðg1 þ l þ c1 þ sÞS 1 þ g2 S 2 þ a1 Qt ; ¼ N dt dS 2 hbP ðS 1 þ wS 2 Þ ðg2 þ l þ c2 þ snÞS 2 þ g1 S 1 þ a2 Qt ; ¼ N dt dQt ¼ c1 ð1 r1 ÞS 1 þ c2 ð1 r2 ÞS 2 ðl þ a1 þ a2 ÞQt ; dt dQp ¼ c1 r1 S 1 þ c2 r2 S 2 lQp ; dt dC ¼ sS 1 þ snS 2 lC dC; dt
ð13Þ
where, now, the total population is N(t) = P(t) + S1(t) + S2(t) + Qt(t) + Qp(t) + C(t). In (13), a proportion, 1 h (with 0 < h < 1), of new smokers are assumed to be mild smokers, and the remaining fraction, h, are chain smokers. It is further assumed that chain smokers have a higher probability of generating more new smokers (by a factor w P 1) relative to the mild smokers. This may be true due to a possible higher social influence the chain smokers may have (amongst peers in the smoking environment) over mild smokers, coupled with the possibility that their high smoking frequency may tempt people around them (who are probably passive smokers anyway) to become smokers. Mild smokers become chain smokers (at a rate g1) and the reverse process occurs at a rate g2. While mild smokers develop smoking-related illnesses at a rate s, those with chain smoking habit are assumed to do so at an increased rate, sn (where n > 1). That is, it is assumed that those in the chain smoking class have a higher probability of developing smoking-related illnesses in comparison to those in the mild smoking class. Individuals with such illnesses die at an increased mortality rate, d. A fraction, 1 ri (i = 1, 2), of individuals in the smoking classes S1 and S2 temporarily quit smoking, at rates c1 and c2, respectively; whilst the remaining fraction, ri, quit smoking permanently. Individuals who temporarily quit smoking revert to mild smoking habit at a rate a1, and to chain smoking habit at a rate a2. Finally, individuals in all classes suffer a natural death at a rate l. A schematic description of the extended model (13) is given in Fig. 1b. 3.1. Basic properties of the extended model Here, too, it is plausible to assume that all the state variables and parameters of the extended model (13) are non-negative. Consider the biologically-feasible region M ¼ fðP ; S 1 ; S 2 ; Qt ; Qp ; CÞ 2 R6þ : P þ S 1 þ S 2 þ Qt þ Qp þ C 6 P=lg: The rate of change of the total population, obtained by adding all the equations in the extended model (13), is given by dN ¼ P lN dC: dt
ð14Þ
It is worth noting here that, unlike in the case of the basic models (1)–(4), the total population N(t) is not constant (except if d = 0). Since the right hand side of (14) is bounded by P lN, a standard comparison theorem can be used to show that N ðtÞ 6 N ð0Þelt þ Pl ð1 elt Þ. In particular, N ðtÞ 6 Pl if N ð0Þ 6 Pl. Thus,
486
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
M is positively-invariant. Hence, it is sufficient to consider the dynamics of the flow generated by (13) in M. It is easy to see, by comparison theorem, that lim inft ! 1N(t) 6 P/l. Thus, the omega limit sets of all solutions in R6þ are contained in M. That is, solutions in M remain in M for all time, and those outside M (but in R6þ ) are eventually attracted to M. 3.2. Stability of the SFE of the Extended Model 3.2.1. Local stability The extended model (13) has a smoking-free equilibrium, given by X0 ¼ ðP ; S 1 ; S 2 ; Qt ; Qp ; C Þ ¼ ðP=l; 0; 0; 0; 0; 0Þ: Here, now, the next generation matrices, F and V, are given 1 0 0 J1 ð1 hÞb ð1 hÞwb 0 0 C B B hwb 0 0C B hb B g C and V ¼ B 1 F ¼B B C B J 0 0 0 0A @ @ 3 0
0
0 0
s
by
1
g2
a1
J2
a2
J4
J5
C 0C C; 0C A
ns
0
J8
0
where, J 1 ¼ g1 þ l þ c1 þ s; J 4 ¼ c2 ð1 r2 Þ;
J 2 ¼ g2 þ l þ c2 þ sn;
J 5 ¼ l þ a1 þ a2 ;
J 3 ¼ c1 ð1 r1 Þ;
J 6 ¼ c1 r1 ;
J 7 ¼ c 2 r2 ;
J 8 ¼ l þ d:
It follows that, using the next generation approach, the smokers generation number of the extended model (13), denoted by Rsd ¼ qðFV 1 Þ; is Rsd ¼
bðD1 þ D2 Þ ; D3
with, D1 ¼ h½g2 J 5 þ a1 J 4 þ wðJ 1 J 5 a1 J 3 Þ; D2 ¼ ð1 hÞ½J 2 J 5 a2 J 4 þ wða2 J 3 þ g1 J 5 Þ; D3 ¼ J 1 J 2 J 5 g2 a2 J 3 a1 J 2 J 3 g1 g2 J 5 g1 a1 J 4 a2 J 1 J 4 : It is shown (see Appendix C) that D1 ; D2 ; and D3 are all positive, so that Rsd > 0: Thus, the following result is established by Theorem 2 of [28]. Lemma 2. The SFE, X0 , of the extended model (13) is LAS whenever Rsd < 1, and unstable if Rsd > 1. The quantity Rsd measures the average number of new smokers generated by a single smoker in a population of completely susceptible (potential smokers) individuals. If this number (Rsd ) is less than unity, then a small influx of smokers will not generate large number (‘‘outbreaks’’) of smokers. If Rsd > 1, on the other hand, such influx would result in large number of new smokers. It is easy to see that when g1 ¼ g2 ¼ c2 ¼ a2 ¼ s ¼ n ¼ h ¼ w ¼ d ¼ r2 ¼ 0, and c1 ¼ c; a1 ¼ a; r1 ¼ r, then Rsd ¼ Rs (i.e., the extended model (13) reduces to the basic models (1)–(4), in this case). A global stability proof for X0 is given below. 3.2.2. Global stability Theorem 4. The SFE of the extended model (13) is GAS in M if Rsd < 1. Proof. Consider the Lyapunov function: L ¼ l1 S 1 þ l 2 S 2 þ l 3 Q t ; where,
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
487
l1 ¼ J 2 J 5 a2 J 4 þ wða2 J 3 þ g1 J 5 Þ; l2 ¼ g2 J 5 þ a1 J 4 þ wðJ 1 J 5 a1 J 3 Þ; l3 ¼ a1 ðwg1 þ J 2 Þ þ a2 ðg2 þ wJ 1 Þ: Thus,
_ ¼ l1 ð1 hÞbP ðS 1 þ wS 2 Þ J 1 S 1 þ g2 S 2 þ a1 Qt L N hbP ðS 1 þ wS 2 Þ J 2 S 2 þ g1 S 1 þ a2 Qt þ l3 ðJ 3 S 1 þ J 4 S 2 J 5 Qt Þ þ l2 N S 1 þ wS 2 ¼ ½bP ðD1 þ D2 Þ ND3 N bP ðD1 þ D2 Þ ¼ D3 ðS 1 þ wS 2 Þ 1 ND3 bðD1 þ D2 Þ 6 D3 ðS 1 þ wS 2 Þ 1 since P 6 N inM D3 ¼ D3 ðS 1 þ wS 2 ÞðRsd 1Þ < 0
for Rsd < 1:
Hence, using similar argument as in the proof of Theorem 1, the SFE of the extended model, X0 , is globally-asymptotically stable whenever Rsd < 1. h An alternative proof for Theorem 4, using comparison theorem, is given in Appendix D. 3.3. SPE of the extended model 3.3.1. Existence of smoking-present equilibrium Using the method described in Section 2.3, with G ¼
bðS 1 þ wS 2 Þ ; N
ð15Þ
it follows that the extended model has a unique SPE, denoted by X1 ¼ ðP ; S 1 ; S 2 ; Qt ; Qp ; C Þ; where,
P ¼ Q t
P ; G þl
S 1 ¼
PG H 2 ; ¼ D3 ðG þ lÞ
PG H 0 ; D3 ðG þ lÞ
Q p
S 2 ¼
PG H 1 ; D3 ðG þ lÞ
PG ðH 11 þ H 12 Þ ; ¼ D3 lðG þ lÞ
C
PG sðH 21 þ H 22 Þ ¼ D3 J 8 ðG þ lÞ
ð16Þ
with, H 0 ¼ ½ðJ 5 J 2 J 4 a2 Þð1 hÞ þ ða1 J 4 þ J 5 g2 Þh; H 1 ¼ ½ðg1 J 5 þ a2 J 3 Þð1 hÞ þ ðJ 5 J 1 J 3 a1 Þh; H 2 ¼ ½ðJ 2 J 3 þ g1 J 4 Þð1 hÞ þ ðJ 1 J 4 þ g2 J 3 Þh; H 11 ¼ ½ðJ 4 a1 þ J 5 g2 ÞJ 6 h þ ðg1 J 5 þ J 3 a2 ÞJ 7 ð1 hÞ; H 12 ¼ ½ðJ 5 J 1 J 3 a1 ÞJ 7 h þ ðJ 5 J 2 J 4 a2 ÞJ 6 ð1 hÞ; H 21 ¼ ½ðJ 5 g2 þ a1 J 4 Þh þ ðg1 J 5 þ J 3 a2 Þnð1 hÞ; H 22 ¼ ½ðJ 5 J 1 J 3 a1 Þnh þ ðJ 5 J 2 J 4 a2 Þð1 hÞ: Using (16) in (15), and simplifying, gives bG D3
G ¼
1þ
G D3
½H 0 þ wH 1
22 Þ 12 ½H 0 þ H 1 þ H 2 þ H 11 þH þ sðH 21JþH l 8
:
ð17Þ
488
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
As before, solving for G** from (17) and substituting the result into (16) gives the corresponding SPE of the model (13). Here, too, G** = 0 corresponds to the smoking-free equilibrium, X0 . For G** 5 0 (i.e., we are looking for an SPE of (13)), (17) can be simplified to c0 G c1 ¼ 0; where, c0 ¼
ð18Þ
1 H 11 þ H 12 sðH 21 þ H 22 Þ þ H0 þ H1 þ H2 þ and D3 J8 l
c1 ¼ Rsd 1:
From (18), G ¼ cc10 . It follows that since the terms in c0 are all positive, then the coefficients c0 and c1 are positive for Rsd > 1. Thus, the extended model (13) has a unique SPE whenever Rsd > 1. It is further evident that if Rsd < 1, then the coefficient c1 < 0, so that G** < 0 (which is biologically meaningless). Here, too, G** = 0 whenever Rsd ¼ 1. Thus, the SPE ðX1 Þ does not exist whenever Rsd 6 1. These results are summarized below. Theorem 5. The extended model (13) has a unique positive smoking-present equilibrium if Rsd > 1, and no positive smoking-present equilibrium if Rsd 6 1. 3.3.2. Local stability of the SPE of the extended model We claim the following result. Theorem 6. The unique smoking-present equilibrium of the extended model (13) is LAS whenever Rsd > 1 and is close to 1. Proof. Linearizing the extended model (13) around its unique SPE is laborious (and not really tractable) owing to its high dimensionality. Consequently, an alternative approach is considered. It is based on the use of the centre manifold theory [5], as described in [7, (Theorem 4.1)]. To apply this method, the following simplification and change of variables are made first of all. Let P ¼ x1 ; S 1 ¼ x2 ; S 2 ¼ x3 ; Qt ¼ x4 ; Qp ¼ x5 , and C = x6, so that N ¼ x1 þ x2 þ x3 þ x4 þ x5 þ x6 . Further, by using vector notation X ¼ ðx1 ; x2 ; x3 ; x4 ; x5 ; x6 ÞT , T ¼ ðf1 ; f2 ; f3 ; f4 ; f5 ; f6 Þ , as follows: the extended model (13) can be written in the form dX dt dx1 dt dx2 dt dx3 dt dx4 dt dx5 dt dx6 dt
bx1 ðx2 þ wx3 Þ lx1 ; x1 þ x2 þ x3 þ x4 þ x5 þ x6 bð1 hÞx1 ðx2 þ wx3 Þ ¼ f2 ¼ J 1 x2 þ g1 x3 þ a1 x4 ; x1 þ x2 þ x3 þ x4 þ x5 þ x6 bhx1 ðx2 þ wx3 Þ ¼ f3 ¼ J 2 x3 þ g1 x2 þ a2 x4 ; x1 þ x2 þ x3 þ x4 þ x5 þ x6 ¼ f1 ¼ P
ð19Þ
¼ f4 ¼ J 3 x2 þ J 4 x3 J 5 x4 ; ¼ f5 ¼ J 6 x2 þ J 7 x3 lx5 ; ¼ f6 ¼ sx2 þ snx3 J 8 x6 :
The Jacobian of the system (19) at the SFE, X0 , is 0 l b bw B 0 ð1 hÞb J ð1 hÞbw þ g1 B 1 B B 0 hb þ g1 hbw J 2 B J ðX0 Þ ¼ B B 0 J4 J3 B B 0 J6 J7 @ 0 s sn
given by 0
0
a1
0
a2
0
J 5
0
0
l
0
0
0
1
C 0 C C 0 C C C; 0 C C 0 C A J 8
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
489
from which it can also be shown (as before) that Rsd ¼
bðD1 þ D2 Þ : D3
Consider the case when Rsd ¼ 1. Suppose, further, that b = b* is chosen as a bifurcation parameter. Solving for b from Rsd ¼ 1 gives b ¼ b ¼
D3 : D1 þ D2
The right eigenvectors of J(b*), the Jacobian of (19) evaluated at b = b*, are given by w ¼ ½w1 ; w2 ; w3 ; w4 ; w5 ; w6 T , where b w2 b ww4 ; w2 ¼ w2 > 0; w3 ¼ w3 > 0; l J 3 w2 þ J 4 w3 J 6 w2 þ J 7 w3 sw2 þ snw3 ; w6 ¼ ; w5 ¼ : w4 ¼ J5 l J8
w1 ¼
Further, ðJ b Þ has a left eigenvector v ¼ ½v1 ; v2 ; v3 ; v4 ; v5 ; v6 ; where v1 ¼ 0;
v2 ¼ v2 > 0;
v3 ¼ v3 > 0;
v4 ¼
a1 v 2 þ a2 v 3 ; J5
v5 ¼ 0;
v6 ¼ 0:
Note that the transformed system (19), with b = b*, has at least one non-hyperbolic equilibrium point (i.e., the linearized system has at least one eigenvalue with zero real part). Hence, the center manifold theory can be used to analyse the dynamics of (19), with b = b*. In particular, the following theorem (Theorem 4.1 in [7]) will be used to show the LAS of the SPE of (19) (which is the same as the SPE of the original system 13). The theorem in [7] is reproduced below for convenience. Theorem 7 (Castillo-Chavez and Song [7]). Consider the following general system of ordinary differential equations with a parameter / dx ¼ f ðx; /Þ; dt
f : Rn R ! R
and
f 2 C2 ðRn RÞ;
where 0 is an equilibrium point of the system (i.e., f(0, /) 0 for all /) and assume ofi ð0; 0ÞÞ is the linearization matrix of the system around the equilibrium 0 with / evaluA1: A ¼ Dx f ð0; 0Þ ¼ ðox j ated at 0. Zero is a simple eigenvalues of A and other eigenvalues of A have negative real parts; A2: Matrix A has a right eigenvectors w and a left eigenvector v corresponding to the zero eigenvalue.
Let fk be the kth component of f and a¼
n X
vk wi wj
k;i;j¼1
b¼
n X k;i¼1
vk wi
o2 fk ð0; 0Þ; oxi oxj
o2 fk ð0; 0Þ: oxi o/
The local dynamics of the system around 0 is totally by a and b. (i) a > 0; b > 0. When / < 0 with j/j 1; 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0 < / 1, 0 is locally asymptotically stable equilibrium; (ii) a < 0; b < 0. When / < 0 with j/j 1; 0 is unstable; when 0 < / 1, 0 is locally asymptotically stable equilibrium, and there exists a positive unstable equilibrium; (iii) a > 0; b < 0. When / < 0 with j/j 1; 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0 < / 1, 0 is unstable, and a positive unstable equilibrium appears;
490
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
(iv) a < 0; b > 0. When / changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable. To apply the above theorem, the following computations are necessary. Computations of a and b For the system (19), the associated non-zero partial derivatives of the right hand side functions (fi) are given by o2 f2 ð1 hÞð1 þ wÞb l o2 f3 hð1 þ wÞb l ; ; ¼ ¼ P P ox2 ox3 ox2 ox3 o2 f2 o2 f2 o2 f2 o2 f2 ð1 hÞb l ; ¼ 2 ¼ ¼ ¼ P ox2 ox4 ox2 ox2 ox5 ox2 ox6 o2 f2 o2 f2 o2 f2 o2 f2 ð1 hÞb wl ; ¼ 2 ¼ ¼ ¼ P ox3 ox4 ox3 ox3 ox5 ox3 ox6 o2 f3 o2 f3 o2 f3 o2 f3 hb l ; ¼ 2 ¼ ¼ ¼ P ox2 ox4 ox2 ox2 ox5 ox2 ox6 o2 f3 o2 f3 o2 f3 o2 f3 hb wl ; ¼ 2 ¼ ¼ ¼ P ox3 ox4 ox3 ox3 ox5 ox3 ox6 o2 f 2 ¼ 1 h; ox2 ob
o2 f 2 ¼ ð1 hÞw; ox3 ob
o2 f3 ¼ h; ox2 ob
o2 f3 ¼ hw: ox3 ob
Thus, it follows from the above expressions that a¼
6 X
vk wi wj
k;i;j¼1
¼
6 X
v2 wi wj
i;j¼1
o2 f k ð0; 0Þ oxi oxj
6 X o2 f2 o2 f3 þ v3 wi wj oxi oxj i;j¼1 oxi oxj
2b lðw2 þ w3 þ w4 þ w5 þ w6 Þðw2 þ ww3 Þðv2 ð1 hÞ þ v3 hÞ < 0; ¼ P and b¼
6 X
vk w i
o2 fk ð0; 0Þ oxi o/
v2 wi
6 X o2 f2 o2 f3 þ v3 w i oxi o/ i¼1 oxi o/
k;i¼1
¼
6 X i¼1
¼ ½v2 ð1 hÞ þ v3 hðw2 þ ww3 Þ > 0: Thus, a < 0 and b > 0. Hence, by Theorem 7(iv) above, the unique SPE of the extended model, which exists whenever Rsd > 1 (Theorem 6), is LAS whenever Rsd > 1 and b* < b with b close to b*. h It is worth mentioning that a similar approach is used in [22] to establish the local stability of endemic equilibrium. Here, too, extensive numerical simulations (see Fig. 4b, d and f) suggest the following conjecture: Conjecture 2. The unique smoking equilibrium of model 13 is globally-asymptotically stable in M n M0 , with M0 ¼ fðP ; S 1 ; S 2 ; Qt ; Qp ; CÞ 2 M : S 1 ¼ S 2 ¼ Qt ¼ 0g, whenever Rsd > 1. Overall, the qualitative analyses in Sections 2 and 3 show that both the basic models, (1)–(4), and the extended model (13) exhibit the same dynamical behaviour (each has a globally-stable SFE whenever its associated smoking-generation number is less than unity; and has a unique LAS SPE whenever this quantity
800
800
700
700
600
600
Mild Smokers S1(t)
Mild Smokers S1(t)
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
500 400 300
500 400 300
200
200
100
100
0
0
2
4
6
8
10
12
14
0 0
16
2
4
6
500
500
450
450
400
400
350
350
300 250 200 150
0 10
12
14
0 0
16
2
4
6
500
500
450
450
400
400
350 300 250 200 150
8
16
10
12
14
16
200 150
50
0
Time(years)
14
250
100
6
12
300
50
4
10
350
100
2
8
Time(years)
Temporary Quitters Qt(t)
Temporary Quitters Qt(t)
Time(years)
0
16
150
50
8
14
200
100
6
12
250
50
4
10
300
100
2
8
Time(years)
Chain Smokers S2(t)
Chain Smokers S2(t)
Time(years)
0
491
10
12
14
16
0 0
2
4
6
8
Time(years)
Fig. 4. Time series plot of the extended model (13) with different initial conditions. Parameter as in Table 3. (a) S1(t) for Rsd < 1; (b) S1(t) for Rsd > 1; (c) S2(t) for Rsd < 1; (d) S2(t) for Rsd > 1; (e) Qt(t) for Rsd < 1; (f) Qt(t) for Rsd > 1.
exceeds unity). In other words, adding the compartments (and dynamics) of the mild and chain smokers, as well as the class of smoking-related illnesses, to the basic model does not alter its dynamical features. The public health implication of these results is that the number of smokers in a community, together with the
492
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
associated smoking-related mortality, can be effectively controlled if the smoking-generation threshold is brought down to a value less than unity. On the other hand, if this threshold exceeds unity, then the number of smokers and smoking-related mortality will rise. 3.4. Numerical simulations of the extended model The extended model (13) is simulated using the parameter values in Table 3. As evident from Table 4, the solution profiles converge to X0 whenever Rsd < 1, and to X1 for Rsd > 1 (in line with Theorems 4 and 6). In addition to showing increasing number of smokers and smoking-related mortality (at steady-state) for increasing Rsd values, this table further shows that large values of the effective contact rate, b, could lead to drastic reductions in the entire population (including the near-complete absence of potential smokers). Time series plots for cases where Rsd < 1 and Rsd > 1 are depicted in Fig. 4, further supporting Theorems 4 and 6. Simulations for cumulative mortality and mortality per unit time are also carried out. Fig. 5a shows the cumulative death resulting from the illnesses for Rsd < 1; suggesting that about 38% of the total population would have acquired smoking habits at steady-state (which is attained in about six years time). Fig. 5b and c depict the time series plots of the mortality per unit time for the case Rsd < 1 and Rsd > 1, respectively, showing a peak mortality of approximately 230 (for Fig. 5b and 260 (Fig. 5c) people. These figures show reductions in smoking-related mortality with decreasing values of Rsd . Table 3 Description of parameters for the extended model (13) Parameter
Description
Estimated value/range
Nominal value
P 1/l a1 a2 c1 c2 r1 r2 g1 g2 b s n d w h
Recruitment rate into the potential smokers class Average duration of smoking activity Rate at which temporary quitters revert back to S1 class Rate at which temporary quitters revert back to S2 class Quitting rate for mild smokers Quitting rate for chain smokers Fraction of quitters from S1 class that permanently quit smoking Fraction of quitters from S2 class that permanently quit smoking Rate at which mild smokers become chain smokers Rate at which chain smokers become mild smokers Contact rate (b = cq) Rate at which mild smokers develop smoking-related illnesses Modification parameter Smoking-related mortality rate Modification parameter Fraction of new smokers that begin chain smoking
30 per year 20 years a1 2 (0, 1) a2 2 (0, 1) c1 2 (0, 1) c2 2 (0, 1) r1 2 ð0; 1Þ r2 2 ð0; 1Þ g1 2 ð0; 1Þ g2 2 ð0; 1Þ b 2 ð0; 1Þ s 2 ð0; 1Þ n>1 d 2 ð0; 1Þ wP1 h 2 ð0; 1Þ
30 20 0.3 0.3 0.4 0.5 0.4 0.5 2 3 0.2 1 3 9 9 0.6
Table 4 Effect of b on Rsd and P ; S 1 ; S 2 ; Qt ; Qp ; C at steady-state using model (13) b
Rsd
P
S1
S2
Qt
Qp
C
Stable equilibrium
0.003 0.03 0.02 0.1 0.2 0.3 0.5 1 2 3 5 7
0.0061 0.0605 0.0403 0.2017 0.4034 0.6051 1.0085 2.0169 4.0338 6.0507 10.084 14.1183
600 600 600 600 600 600 571 86 32 20 11 8
0 0 0 0 0 0 0 8 8 9 9 9
0 0 0 0 0 0 0 5 5 6 6 6
0 0 0 0 0 0 0 5 5 5 5 5
0 0 0 0 0 0 3 49 54 55 56 57
0 0 0 0 0 0 1 21 24 24 24 25
X0 X0 X0 X0 X0 X0 X1 X1 X1 X1 X1 X1
250
250
200
200
150
150
Mortality δC
Cummulative Death
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
100
50
0 0
493
100
50
2
4
6
8
10
12
14
0 0
16
2
4
6
Time(years)
8
10
12
14
16
Time(years)
300
250
Mortality δC
200
150
100
50
0
0
2
4
6
8
10
12
14
16
Time(years)
Fig. 5. Simulations of the extended model (13). (a) Cumulative mortality; (b) total death for Rsd < 1; (c) total death for Rsd < 1.
In order to ascertain the impact of the relative sojourn (duration) in the two smoking classes on the overall smoking dynamics, more simulations are carried out for a 3 year time period as follows. In Fig. 6, contours of Rsd are plotted as a function of number of days spent in the S1 class ðc11 Þ and the S2 class ðc12 Þ. Based on the parameter values used in these simulations, this figure shows that if chain smokers do not spend more than 1.5 years in this class, the number of smokers in the community (at steady state) will be effectively controlled or eliminated (since Rsd < 1 in this case) regardless of the time spent by smokers in the mild smoking class (S1). Similarly, if mild smokers practice their mild smoking habit for less than 1.5 years, then the number of smokers in the community (at steady state) will be greatly reduced (or eliminated) from the population irrespective of the dynamics in the chain smoking class. Furthermore, the number of smokers in the community (at steady state) will persist if individuals in either classes remain there for a duration longer than 1.5 years (here Rsd > 1). Of course, it should be emphasized that these quantitative estimates are dependent on the parameter values chosen, and, owing to expected uncertainties in these parameters, the results should be taken only as a general guide. The impact of the transition from mild to chain smoking class is also monitored by plotting Rsd contours as a function of the number of days a smoker of class i spends in class j (i 5 j). Fig. 7a shows a marked decrease in Rsd for longer duration in the mild smoking class by individuals in the chain smoking class. In other words, the number of smokers will be greatly reduced if chain smokers revert to mild smoking, and stay in the mild smoking class for a long duration. On the other hand, Rsd increases if mild smokers move to the chain
494
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
0.9
0.8
0.7 3
1
2.5
2 2
1/γ
0.9
0.8
0.7
1
1.5 1
1
0.9 0.
8
0.7
0.5
0.9
0.8
0.8
07
0.5
1
1.5
2
1/γ
2.5
3
1
Fig. 6. Contour plots of the extended model (13) depicting Rsd as a function of
1 c1
and c12 . Parameters as in Table 3.
smoking class and stay there for long duration (Fig. 7b). Thus, it is more beneficial (in terms of public health) for chain smokers to acquire mild smoking habits for extended duration, than mild smokers becoming chain smokers for a long periods. 3.5. Conclusions A four-dimensional compartmental deterministic model is designed and used to monitor smoking dynamics in a population. This model was extended to incorporate the heterogeneity of smoking frequency and smoking-related illnesses. The models were rigorously analysed to gain insight into their dynamical features. The study shows the following: (i) Each of the two models considered in this study has a globally-stable smoking-free equilibrium if a certain threshold quantity, known as the smokers-generation number, is less than unity; indicating that the number of smokers in the community will be brought to zero if public health measures that make (and keep) the threshold to a value less than unity are carried out; (ii) For the basic model, the number of smokers can be reduced by reducing the contact rate (b), reversion rate (a) and increasing quitting rates (c) and the fraction (r) of temporary quitters who permanently giveup smoking; (iii) For the extended model, the number of smokers can similarly be reduced by reducing the contact rates (b and bw), transition rates, (g1 and g2), reversion rates (a1 and a2), and increasing the rate at which smokers temporarily quit smoking (c1 as well as the c2) and fraction of permanent quitters (r1 and r2); (iv) The transition from chain to the mild smoking class plays a crucial role in the smoking dynamics. More new smokers are generated if mild smokers adopt a chain smoking habit; and markedly reduced number of new smokers are recorded if chain smokers move to the mild smoking class; (v) If mild and chain smokers spend less than 1.5 years in their respective smoking classes, then the number of smokers in the community will be effectively controlled (or eliminated). The number of smokers will persist otherwise. Similarly, if chain smokers spend less than 1.5 years in this class, the number of smokers will be effectively controlled (as the same reason applies).
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
495
1.3
1.25
1.2
1.15
Rsd
1.1
1.05
1
0.95
0.9
0.85
0.8 0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
Duration in mild smoking class by a chain smoker: (η /(η + μ + γ + τ)) 2
1
1
1.4
1.35
1.3
1.25
R
sd
1.2
1.15
1.1
1.05
1
0.95
0.9 0
0.5
1
1.5
2
2.5
Duration in chain smoking class by a mild smoker: (η /(η + μ + γ + τξ)) 1
2
2
Fig. 7. Simulations of the extended model (13). (a) Plot of the smokers generation number Rsd as a function of the number of days (years) mild smokers spend in the chain smoking class S2; (b) Plot of the smokers generation number Rsd as a function of the number of days (years) chain smokers spend in the mild smoking class S1.
In summary, although the models considered here can further be refined to incorporate other smokingrelated phenomena, such as the impact of second-hand (passive) smoking, recruitment of smokers from external sources into the population (so that a fraction of recruited people, P, go directly into the S1 and S2 classes, etc), age structure, smoking advertising and mass consumption, this study shows that smoking and smoking-related illnesses can be effectively controlled in a community if public health measures that can bring the smokers generation number to a value less than unity are implemented.
496
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
Acknowledgements One of the authors (ABG) acknowledges, with thanks, the support in part of the Natural Science and Engineering Research Council (NSERC) and Mathematics of Information Technology and Complex Systems (MITACS) of Canada. Appendix A. Next generation method [28] The next generation method is used to establish the local asymptotic stability of the SFE [11,28]. The method was first introduced by Diekmann and Hesterbeek [11], and refined for epidemiological models by van den Driessche and Watmough [28]. The formulation in [28], which is based on disease transmission model, is reproduced below (note that this method can be implemented on the smoking models in this paper by essentially replacing the word ‘‘disease’’ with ‘‘smoking’’). Suppose the given disease transmission model, with nonnegative initial conditions, can be written in terms of the following system: x_ i ¼ fi ðxÞ ¼ Fi ðxÞ Vi ðxÞ;
i ¼ 1; . . . ; n;
ð20Þ
þ where Vi ¼ V and the function satisfy the five axioms below. First of all, X s ¼ fx P 0jxi ¼ i Vi 0; i ¼ 1; . . . ; mg with m < n is defined as the disease-free states (non-infected state variables of the model) of the model, where x ¼ ðx1 ; . . . ; xn Þt ; xi P 0 represents the number of individuals in each compartment of the model.
(A1) (A2) (A3) (A4) (A5)
if x P 0, then Fi ; Vþ i ; Vi P 0 for i ¼ 1; . . . ; m. if xi = 0, then Vi ¼ 0. In particular, if x 2 Xs then V i ¼ 0 for i ¼ 1; . . . ; m. Fi ¼ 0 if i > m. if x 2 Xs, then Fi ðxÞ ¼ 0 and Vþ i ðxÞ ¼ 0 for i ¼ 1; . . . ; m. If FðxÞ is set to zero, then all eigenvalues of Df(x0) have negative real parts.
Here, Fi ðxÞ represents the rate of appearance of new infections in compartment i; Vþ i ðxÞ represents the rate of transfer of individuals into compartment i by all other means, and V i ðxÞ represents the rate of transfer of individuals out of compartment i. It is assumed that these functions are at least twice continuously differentiable in each variable [28]. Lemma 3 (van den Driessche and Watmough [28]). If x is a DFE of 20 and fi(x) satisfy (A1)–(A5), then the derivatives DFðxÞ and DVðxÞ are partitioned as F 0 V 0 DFðxÞ ¼ ; ; DVðxÞ ¼ 0 0 Q3 Q4 where F and V are the m · m matrices defined by, h i h i oF i ðxÞ and V ¼ oV F ¼ oxji ðxÞ with 1 6 i; j 6 m: oxj Further, F is non-negative, V is a non-singular M-matrix and Q3 ; Q4 are matrices associated with the transition terms of the model, and all eigenvalues of Q4 have positive real parts. Finally, the following stability result follows. Theorem 8 (van den Driessche and Watmough [28]). Consider an arbitrary disease transmission model with fi(x) satisfying the axioms (A1)–(A5). If x is a DFE of the model, then x is LAS if R0 ¼ qðFV 1 Þ < 1, but unstable if R0 > 1. In the above, q represents the spectral radius.
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
497
Appendix B. Proof of Theorem 1 using comparison theorem Theorem B.1. The SFE of the basic models (1)–(4), given by E0 , is GAS in D if Rs < 1. Proof. The above result will now be proved using a comparison theorem. As in Section 2.2.1, the equations for the rate of change of the smoking components (S, Qt) of the model (1)–(4), are re-written in terms of ! dSðtÞ SðtÞ b 0 SðtÞ P dt ¼ ðF V Þ 1 ; dQt ðtÞ N Qt ðtÞ 0 0 Qt ðtÞ dt
where the matrices F and V are as defined in Section 2.2.1. Since P 6 N (for all t P 0) in D, it follows that ! dSðtÞ SðtÞ dt 6 ðF V Þ : ðB:1Þ dQt ðtÞ Qt ðtÞ dt
Using the fact that the eigenvalues of the matrix F V all have negative real parts (see local stability result in Lemma 1, where q(FV1) < 1 if Rs < 1, which is equivalent to F V having eigenvalues with negative real parts when Rs < 1 [28]), it follows that the linearized differential inequality system (B.1) is stable whenever Rs < 1. Consequently, ðSðtÞ; Qt ðtÞÞ ! ð0; 0Þ as t ! 1. It follows by comparison theorem (see, for instance, [23, p. 31] and [26], Theorem B.1; Appendix B) that ðSðtÞ; Qt ðtÞÞ ! ð0; 0Þ. Substituting S = Qt = 0 in the first and the last equations of model (1)–(4), gives P(t) ! P*, and Qp(t) ! 0 as t ! 1. Thus, ðP ðtÞ; SðtÞ; Qt ðtÞ; Qp ðtÞÞ ! ðP ; 0; 0; 0Þ as t ! 1 for Rs < 1, so that E0 is GAS if Rs < 1. h Appendix C. Positivity of D1, D2 and D3 Given, D1 ¼ h½g2 J 5 þ a1 J 4 þ wðJ 1 J 5 a1 J 3 Þ; D2 ¼ ð1 hÞ½J 2 J 5 a2 J 4 þ wða2 J 3 þ g1 J 5 Þ; it follows that J 1 J 5 a1 J 3 ¼ ðg1 þ l þ c1 þ sÞðl þ a1 þ a1 Þ a1 c1 ð1 r1 Þ ¼ l2 þ ðg1 þ a1 þ s þ a2 þ c1 Þl þ g1 a1 þ g1 a2 þ c1 a2 þ sa1 þ sa2 þ c1 a1 r1 > 0 and J 2 J 5 a2 J 4 ¼ ðg2 þ l þ c2 þ snÞðl þ a1 þ a1 Þ a2 c2 ð1 r2 Þ ¼ l2 þ ðg2 þ a1 þ sn þ a2 þ c2 Þl þ g2 a1 þ g2 a2 þ c2 a1 þ sna1 þ sna2 þ c2 a2 r2 > 0: Thus, D1 > 0 and D2 > 0. Furthermore, D 3 ¼ J 1 J 2 J 5 g2 a2 J 3 a 1 J 2 J 3 g 1 g2 J 5 g1 a1 J 4 a2 J 1 J 4 ¼ ðg1 þ l þ c1 þ sÞðg2 þ l þ c2 þ snÞðl þ a1 þ a1 Þ g2 a2 c1 ð1 r1 Þ a1 c1 ðg2 þ l þ c2 þ snÞð1 r1 Þ g1 g2 ðl þ a1 þ a2 Þ g1 a1 c2 ð1 r2 Þ a2 c2 ðg1 þ l þ c1 þ sÞð1 r2 Þ ¼ l3 þ ðc1 þ a2 þ c2 þ a1 þ sn þ s þ g1 þ g2 Þl2 þ ½g2 a2 þ c1 c2 þ c1 g2 þ c1 a2 þ g2 a1 þ c2 a2 r2 þ c2 a1 þ g1 c2 þ g1 a1 þ g1 a2 þ c1 a1 r1 þ ðc2 þ a2 þ g2 þ a1 Þs þ ðs2 þ fa2 þ c1 þ g1 þ a1 gsÞnl þ ðfa2 þ a1 gs2 þ fg1 a2 þ c1 a2 þ c1 a1 r1 þ g1 a1 gsÞn þ a2 c2 c1 r2 þ ðg2 a1 þ c2 a2 r2 þ g2 a2 þ c2 a1 Þs þ c1 g2 a2 r1 þ a1 c1 g2 r1 þ a1 c1 c2 r1 þ a2 c2 g1 r2 þ g1 c2 a1 r2 > 0 Thus, D3 is also positive.
498
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
Appendix D. Proof of Theorem 4 using comparison theorem Theorem 4. The SFE of the extended model (13) is GAS if Rsd < 1. Proof. Here, too, the equations for the rate of change of the smoking components (S 1 ; S 2 ; Qt ; C) of the model 13 can be written in terms of 0 dS 1 ðtÞ 1 0 1 1 0 10 S 1 ðtÞ ð1 hÞb ð1 hÞwb 0 0 S 1 ðtÞ dt B dS 2 ðtÞ C B C C CB B C B B S 2 ðtÞ C B hb hwb 0 0 CB S 2 ðtÞ C B dt C C 1 P B C CB B dQ ðtÞ C ¼ ðF V ÞB B Q ðtÞ C B CB Q ðtÞ C; B t C N 0 0 0 0 @ A @ A @ A t t @ dt A dCðtÞ CðtÞ CðtÞ 0 0 0 0 dt
where the matrices F and V are as defined in Section 3.2.1. Since P 6 N (for all t P 0) in M, it follows that 0 dS ðtÞ 1 1 0 1 S 1 ðtÞ B dSdtðtÞ C B S ðtÞ C B 2 C B dt C B 2 C ðD:1Þ C: B dQ ðtÞ C 6 ðF V ÞB @ Qt ðtÞ A B t C @ dt A dCðtÞ CðtÞ dt
Using the fact that the eigenvalues of the matrix F V all have negative real parts (see local stability result in Lemma 2, where q(FV1) < 1 if Rsd < 1, which is equivalent to F V having eigenvalues with negative real parts when Rsd < 1 [28]), it follows that the linearized differential inequality system (D.1) is stable whenever Rsd < 1. Consequently, ðS 1 ðtÞ; S 2 ðtÞ; Qt ðtÞ; CðtÞÞ ! ð0; 0; 0; 0Þ as t ! 1. It follows, by comparison theorem, that ðS 1 ðtÞ; S 2 ðtÞ; Qt ðtÞ; CðtÞÞ ! ð0; 0; 0; 0Þ. Substituting S1 = S2 = Qt = C = 0 in the first and penultimate equations of the model (13) gives P(t) ! P* and Qp(t) ! 0 as t ! 1. Thus, ðP ðtÞ; S 1 ðtÞ; S 2 ðtÞ; Qt ðtÞ; Qp ðtÞ; CðtÞÞ ! ðP ; 0; 0; 0; 0; 0Þ as t ! 1 for Rsd < 1, so that X0 is GAS if Rsd < 1. h References [1] R.M. Anderson, R.M. May (Eds.), Population Biology of Infectious Diseases, Springer-Verlag, Berlin, Heidelberg, New York, 1982. [2] R.M. Anderson, R.M. May, Infectious Diseases of Humans, Oxford University Press, London, 1991. [3] D. Bernoulli, Essai d’une nouvelle analyse de la mortalite´ cause´e par la petite ve´role, et des avantages de l’inoculation pour la pre´venir. Histoire de l’Acade´mie Royale des Sciences. me´moire de mathe´matiques et de physique, Paris (1760). pp. 1–45. [4] F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Text in Applied Mathematics, Springer, 2000. [5] J. Carr, Applications Centre Manifold Theory, Springer-Verlag, New York, 1981. [6] C. Castillo-Chavez, K. Cooke, W. Huang, S.A. Levin, Results on the dynamics for models for the sexual transmission of the human immunodeficiency virus, Appl. Math. Lett. 2 (1989) 327–331. [7] C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. 1 (2) (2004) 361–404. [8] C. Castillo-Chavez, S. Blower, P. van den Driessche, D. Kirschner, A.A. Yakubu (Eds.), Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Introduction, The IMA Volumes in Mathematics and its Applications, vol. 125, Springer, New York, 2002. [9] C. Castillo-Chavez, H. Wenzhang, L. Jia, Competitive exclusion and coexistence of multiple strains in an SIS STD model, SIAM J. Appl. Math. 5 (1999) 1790–1811. [10] G.C. Castillo, S.G. Jordan, and A.H. Rodriguez. Mathematical models for the dynamics of tobacco use, recovery and relapse. Technical Report Series, BU-1505-M. Department of Biometrics, Cornell University, 2000. [11] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, J. Math. Biol. 28 (1990) 503–522. [12] E.H. Elbasha, A.B. Gumel, Theoretical assessment of public health impact of imperfect Prophylactic HIV-1 vaccines with therapeutic benefits, Bull. Math. Biol. 68 (2005) 577–614. [13] A.B. Gumel, C.C. McCluskey, J. Watmough, An SVEIR model for assessing potential impact of an imperfect anti-SARS vaccine, Math. Biosci. Eng. 3 (3) (2006) 485–512. [14] A.B. Gumel, McCluskey, C. Connell, P. van den Driessche, Mathematical study of a staged-progression HIV model with imperfect vaccine. Bull. Math. Biol., 2006, 10.1007/s11538-006-9095-7.
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
499
[15] A.B. Gumel, X.W. Zhang, P.N. Shivakumar, M.L. Garba, B.M. Sahai, A new mathematical model for assessing therapeutic strategies for HIV infection, J. Theor. Med. 4 (2) (2002) 147–155. [16] J.K. Hale, Ordinary Differential Equations, John Wiley and Sons, New York, 1969. [17] Health Canada . New directions for tobacco control in Canada: a strategy. Steering Committee of the National Strategy to reduce Tobacco Use in Canada, Ottawa, Ontario. Available from: http://dsp-psd.pwgsc.gc.ca/Collection/H39-505-1999E.pdf, 2005. [18] H.M. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599–963. [19] W.O. Kermack, A.G. McKendrick, A contribution to the mathematical theory of epidemics, Proc. Royal Soc. London. 115 (1927) 700–721. [20] W.O. Kermack, A.G. McKendrick, Contributions to the mathematical theory of epidemics, part II, Proc. Royal Soc. London 138 (1932) 55–83. [21] W.O. Kermack, A.G. McKendrick, Contributions to the mathematical theory of epidemics, part III, Proc. Royal Soc. London 141 (1933) 94–112. [22] M. Kgosimore, E.M. Lungu, The effects of vaccination and treatment on the spread of HIV/AIDS, J Biol Syst. 12 (4) (2004) 399–417. [23] V. Lakshmikantham, S. Leela, A.A. Martynyuk, Stability Analysis of Nonlinear Systems, Marcel Dekker Inc., New York, 1989. [24] A. Perelson, P. Nelson, Mathematical analysis of HIV-1 dynamics in vivo, SIAM Rev. 41 (1999) 3–44. [25] T.L. Scott, S. Robert, V.A. Kirsh, Beliefs about tobacco industry (mal)practices and youth smoking behaviour: insight for future tobacco control campaigns (Canada), Cancer Causes Control 17 (2006) 705–711. [26] R.J. Smith, S.M. Blower, Could disease-modifying HIV vaccine cause population-level pervasity? Lancet Infect. Dis. 4 (2004) 636– 639. [27] US Department of Health and Human Services. The health consequences of smoking: a report of the surgeon general. US Department of Health and Human Services, Public Health Service, Centers for Disease Control and Prevention, National Center for Chronic Disease Prevention and Health Promotion, Office of Smoking and Health, Atlanta, GA. Available from: http:// www.cdc.gov/Tobacco/sgr/sgr_2004/consumerpiece/SGR2004_Whatitmeanstoyou.pdf, 2004. [28] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48. [29] J .Walton, History of smoking. http://www.forestonline.org/output/page34.asp, 2003. [30] World Health Organization. Global smoking statistics. http://quitsmoking.about.com/cs/antismoking/a/statistics.htm, 2002. [31] World Health Organization. World cancer report. http://www.who.int/mediacentre/news/releases/2003/pr27/en, 2003.
Applied Mathematics and Computation 195 (2008) 475–499 www.elsevier.com/locate/amc
Curtailing smoking dynamics: A mathematical modeling approach O. Sharomi, A.B. Gumel
*
Department of Mathematics, University of Manitoba, Winnipeg, Manitoba, Canada R3T 2N2
Abstract This paper provides a rigorous mathematical study for assessing the dynamics of smoking and its public health impact in a community. A basic mathematical model, which is a slight refinement of the model presented in [F. Brauer, C. CastilloChavez. Mathematical Models in Population Biology and Epidemiology. Text in Applied Mathematics. Springer, 2000; G.C. Castillo, S.G. Jordan, A.H. Rodriguez. Mathematical models for the dynamics of tobacco use, recovery and relapse. Technical Report Series, BU-1505-M. Department of Biometrics, Cornell University. 2000], is designed first of all. It is based on subdividing the total population in the community into non-smokers, smokers and those smokers who quit smoking either temporarily or permanently. The theoretical analysis of the basic model reveals that the associated smoking-free equilibrium is globally-asymptotically stable whenever a certain threshold, known as the smokers-generation number, is less than unity, and unstable if this threshold is greater than unity. The public health implication of this result is that the number of smokers in the community will be effectively controlled (or eliminated) at steady-state if the threshold is made to be less than unity. Such a control is not feasible if the threshold exceeds unity (a global stability result for the smoking-present equilibrium is provided for a special case). The basic model is extended to account for variability in smoking frequency, by introducing two classes of mild and chain smokers as well as the development and the public health impact of smoking-related illnesses. The analysis and simulations of the extended model, using an arbitrary but reasonable set of parameter values, reveal that the number of smokers in the community will be significantly reduced (or eliminated) if chain smokers do not remain as chain smokers for longer than 1.5 years before reverting to the mild smoking class, regardless of the time spent by mild smokers in their (mild smoking) class. Similarly, if mild smokers practice their mild smoking habit for less than 1.5 years, the number of smokers in the community will be effectively controlled irrespective of the dynamics in the chain smoking class. 2007 Elsevier Inc. All rights reserved. Keywords: Tobacco smoking; Smoking-related illnesses; Equilibrium; Stability; Mild/chain smokers
1. Introduction Since the advent of tobacco in 6000 BC [29], tobacco smoking remains a major public health menace globally. It inflicts significant mortality, morbidity and continues to cost billions of dollars in health-care (see [27] *
Corresponding author. E-mail address: [email protected] (A.B. Gumel).
0096-3003/$ - see front matter 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.05.012
476
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
and the references therein). Recent figures from the World Health Organization (WHO) estimate that 100 million people died of tobacco-associated diseases (such as cancer, chronic lung disease, cardiovascular disease and stroke) last century [31]. Although smoking is widely acknowledged as the most important preventable cause of death, 4 million people die of smoking-related illnesses globally every year, and the number of new smokers continues to rise [30]. Owing to these facts, and the astronomical public health burden associated with smoking, it is of public health significance to qualitatively study the dynamics of smoking in a community, aimed at determining realistic methods for combatting this habit. This would allow for the theoretical assessment of tobacco smoking control and protection from second-hand smoking etc. [25]. Although mathematical modelling has been extensively used to address questions of public health importance, dating back to the pioneering works of Bernoulli (on modelling the dynamics of smallpox) in 1760 [3], Kermack and McKendrick [19–21] and those reported in more recent literature (such as in [1,2,4,6–9,11– 15,18,24,26,28] and the references therein), not much has been done in terms of the mathematical modeling of human social behaviour. In particular, with the exception of the basic model in [4,10], the authors are not aware of any mathematical study (talkless of the rigorous qualitative analyses of the models) for fully examining, and assessing, the impact of smoking in a human population. The aim of this paper is to provide a qualitative study of the dynamics of smoking, vis-a-vis its impact on public health. A basic compartmental model, which subdivides a given population into a number of mutuallyexclusive sub-populations (depending on their smoking status), is designed and qualitatively analysed. This model is then extended to incorporate variability in smoking habits (namely mild or chain smoking) as well as the impact of smoking-related illnesses in the community. 2. Model formulation of basic model The basic model is based on monitoring the dynamics of the sub-populations of potential smokers (nonsmokers); P(t), smokers; S(t), smokers who temporarily quit smoking; Qt(t), and smokers who permanently quit smoking; Qp(t). Thus, the total population at time t is given by N(t) = P(t) + S(t) + Qt(t) + Qp(t). Non-smokers (potential smokers) (P(t)) The class of non-smokers (potential smokers) is increased by the recruitment of individuals into the nonsmoking class (at a rate lN). It is assumed that non-smokers can acquire smoking habits (and become smokers) via effective ‘‘contacts’’ with smokers (at a rate b = cq, where c is the average number of contacts per unit time; and q is the probability of becoming a smoker for a member of the P class following contact with a smoker). In other words, it is assumed that the acquisition of smoking habit is analogous to acquiring disease infection; the mode ‘‘contact’’ a potential smoker makes, the higher the likelihood of such an individual acquiring smoking habit. Although not considered in this study, other factors such as tobacco advertising and mass consumption of tobacco in the community can influence this process (and, hence, should be included in the function governing the generation of new smokers in the community). Finally, potential smokers suffer natural death (at a rate l), so that the rate of change of the population of potential smokers is given by dP bPS ¼ lN lP : dt N
ð1Þ
Smokers (S(t)) The population of smokers is increased when non-smokers begin to smoke (regardless of their frequency of smoking per unit time) and when smokers who temporarily quit smoking revert to smoking (at a rate a). The population is decreased when smokers temporarily quit smoking (at a rate c) and natural death (at the rate l), giving dS bPS ¼ þ aQt ðl þ cÞS: dt N
ð2Þ
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
477
Smokers who temporary quit smoking (Qt(t)) This is generated by a fraction, 1 r, of smokers who temporarily quit smoking (at the rate c). It is decreased by natural death (at the rate l) and reversion to smoking (at the rate a). Thus, dQt ¼ cð1 rÞS ðl þ aÞQt : dt
ð3Þ
Smokers who permanently quit smoking (Qp(t)) This is generated by the remaining fraction, r, of smokers who permanently quit smoking (at the rate c), and decreased by natural death (at the rate l), so that dQp ¼ crS lQp : dt
ð4Þ
Thus, the basic model for smoking dynamics is given by Eqs. (1)–(4), (see flowchart diagram in Fig. 1a). It is an extension of the model in [4,10], by incorporating a class (and dynamics) of smokers who temporarily quit smoking (Qt(t)). Note that this model does not include smoking-related mortality.
2.1. Basic properties of the model Adding Eqs. (1)–(4) gives dN ¼ 0: dt
ð5Þ
Thus, the total population, N, is constant; and equals N*. Since the models (1)–(4) monitors human population, it is plausible to assume that all its state variables and parameters are non-negative for all t P 0. Further, it can be shown that the region D ¼ fðP ; S; Qt ; Qp Þ 2 R4þ : P þ S þ Qt þ Qp ¼ N g; is positively-invariant. Thus, every solution of the models (1)–(4), with initial conditions in D remains there for t > 0. Therefore, the x limit sets of solutions of (1)–(4), in D are contained in D. Furthermore, in D, the
Fig. 1a. Flowchart of the basic models (1)–(4).
478
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
Fig. 1b. Flowchart of the extended model (13).
usual existence, uniqueness and continuation results hold for the system, so that the system (1)–(4), is wellposed mathematically and epidemiologically [18]. Hence, it is sufficient to consider the dynamics of the flow generated by models (1)–(4), in D. 2.2. Stability of smoking-free equilibrium (SFE) 2.2.1. Local stability The models (1)–(4) has an SFE, obtained by setting the right-hand sides of the equations in (1)–(4) to zero, given by E0 : ðP ; S ; Qt ; Qp Þ ¼ ðN ; 0; 0; 0Þ;
with N ¼ P þ S þ Qt þ Qp ¼ P :
ð6Þ
The linear stability of E0 can be established using the next generation method [11,28] (see Appendix A for a detailed description of the next generation method). Note, first of all, that system (1)–(4) satisfies the axioms A1–A5 in Appendix A. Thus, using the notation in [28], the next generation method entails writing the right hand sides of the equations of the rate of change of the smoking-present variables (dS/dt and dQt/dt) in terms of two matrices F and V, where F is a matrix consisting of the smokers-generation terms (i.e., all terms with b in them in the equations for (dS/dt and dQt/dt) and V is an M-matrix consisting of the remaining transition terms in the two equations (it should be recalled that a matrix A is an M-matrix if and only if every off-diagonal entry of A is non-positive and the diagonal entries are all non-negative). That is, Eqs. (2) and (3) are re-written as (note that P* = N*) ! dSðtÞ b 0 K 1 a SðtÞ dt ¼ ; dQt ðtÞ K 2 K 3 0 0 Qt ðtÞ dt SðtÞ ¼ ðF V Þ ; ð7Þ Qt ðtÞ
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
479
with, K 1 ¼ l þ c;
K 2 ¼ cð1 rÞ;
K 3 ¼ l þ a:
That is, for the models (1)–(4), the next generation matrices F and V are given by b 0 K 1 a F ¼ and V ¼ : 0 0 K 2 K 3 Using the next generation method, the local stability of the SFE, given by (6), is based on whether or not q(FV1) < 1, where q is the spectral radius. If q(FV1) < 1, then all eigenvalues of the linearized version of the models (1)–(4), at the SFE, given by ! dSðtÞ SðtÞ dt ¼ ðF V Þ ; ð8Þ dQt ðtÞ Qt ðtÞ dt
have negative real parts, so that the SFE is locally-asymptotically stable (LAS). For q(FV1) > 1, at least one of the eigenvalues of the linearization has a positive real part; thus, the SFE is unstable in this case. Letting Rs ¼ qðFV 1 Þ, it is easy to show that Rs ¼
bK 3 bðl þ aÞ ¼ : K 1 K 3 aK 2 lðl þ aÞ þ cðar þ lÞ
ð9Þ
Consequently, using Theorem 2 of [28] (see also Theorem 8 of Appendix A), the following result is established. Lemma 1. For the system (1)–(4), the SFE, E0 , is LAS if Rs < 1, and unstable if Rs > 1. The quantity Rs is the smokers generation number. It measures the average number of new smokers generated by a single smoker in a population of potential smokers. The local stability result in Lemma 1 implies that for Rs < 1, the total number of smokers in the population can be reduced to zero if the initial sizes of the sub-populations of the model are in the basin of attraction of E0 . That is, a small influx of smokers into the community will not generate large number of smokers if Rs < 1 (but will do so if Rs > 1). To ensure that the effective control (or elimination) of the number of smokers in the community at steady-state is independent of the initial sizes of the sub-populations of the model, it is imperative to show that the SFE is globally-asymptotically stable (GAS). This is done below. 2.2.2. Global stability of SFE We claim the following result: Theorem 1. The SFE is GAS in D if Rs < 1. Proof. First, it should be noted that P 6 N in D for all t P 0. Consider the following Lyapunov function: F ¼ ðl þ aÞS þ aQt ; with Lyapunov derivative (where a dot represents differentiation with respect to t) _ ¼ ðl þ aÞS_ þ aQ_ t F bPS ¼ ðl þ aÞ ðl þ cÞS þ aQt þ a½cð1 rÞS ðl þ aÞQt N bPSðl þ aÞ ðl þ aÞðl þ cÞS þ ðl þ aÞaQt þ acð1 rÞS aðl þ aÞQt N bPSðl þ aÞ ðl þ aÞðl þ cÞS þ acð1 rÞS ¼ N bP ðl þ aÞ ½ðl þ aÞðl þ cÞ acð1 rÞ S ¼ N ¼
480
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
6 fbðl þ aÞ ½lðl þ aÞ þ cðar þ lÞgS for P 6 N bðl þ aÞ 1 S ¼ ½lðl þ aÞ þ cðar þ lÞ lðl þ aÞ þ cðar þ lÞ ¼ ½lðl þ aÞ þ cðar þ lÞðRs 1ÞS < 0
for Rs < 1:
_ < 0 for Rs < 1; with F _ ¼ 0 only if Since all the model parameters are non-negative, it follows that F S = 0. Hence, F is a Lyapunov function on D. Therefore, by the LaSalle’s Invariance Principle [16], every solution to the equations of the models (1)–(4), with initial conditions in D, approaches E0 as t ! 1. h It is worth stating that this theorem can also be proved using other approaches such as the use of comparison theorem (see Appendix B). 2.3. Smoking-present equilibrium (SPE) 2.3.1. Existence of smoking-present equilibrium The existence of positive smoking-present equilibria (steady-states with the smoking-present components, S** and Q ¼ bS , where ** represents the comt , positive) of the models (1)–(4), is explored as follows. Let G N ponent of the SPE at steady-state (with N ¼ P þ S þ Qt þ Qp ). Then, the state variables of the models (1)–(4), can be expressed at an arbitrary SPE, denoted by E1 ¼ ðP ; S ; Q t ; QP Þ, as: P ¼ Q t
lN ; G þ l
S ¼
K 3 lG N ; ðG þ lÞ½lðl þ aÞ þ cðar þ lÞ
K 2 lG N ; ¼ ðG þ lÞ½lðl þ aÞ þ cðar þ lÞ
Q p
crK 3 G N : ¼ ðG þ lÞ½lðl þ aÞ þ cðar þ lÞ
ð10Þ
Substituting the equations in (10) into the expression for G** above, and simplifying, gives the following fixed-point problem: G ¼
bG K 3 l : l½lðl þ aÞ þ cðar þ lÞ þ G ðlK 3 þ lK 2 þ crK 3 Þ
ð11Þ
The positive smoking-present equilibria of (1)–(4) can then be obtained by solving for G** in (11) and substituting the result into (10). Clearly, G** = 0 is a fixed-point of (11), which corresponds to the SFE ðE0 Þ. Since we are only interested in the possible smoking-present equilibria of (1)–(4), we consider only the case where G** 5 0. For G** 5 0, Eq. (11) can be simplified to: a0 G a1 ¼ 0; K 3 ðlþcrÞþlK 2 lðlþaÞþcðarþlÞ
ð12Þ
¼ 1 and a1 ¼ lðRs 1Þ: where a0 ¼ It follows from (12) that G ¼ aa10 . Further, a1 > 0 for Rs > 1. Thus, the linear Eq. (12) has a unique positive solution when Rs > 1 (i.e., G** > 0 iff Rs > 1). It should be noted that if Rs < 1, then a1 < 0. Hence, for Rs < 1, the quantity G** < 0 (which is biologically meaningless). Consequently, the SPE ðE1 Þ does not exist in this case. Furthermore, for Rs ¼ 1, the coefficient a1 = 0; making G** = 0 (corresponding to the SFE). These results are summarized below. Theorem 2. The model (1)–(4), has a unique SPE, E1 , if Rs > 1, and no SPE if Rs 6 1. 2.3.2. Local stability of smoking-present equilibrium The local stability of the unique SPE is now explored. We claim the following: Theorem 3. The unique SPE of the model (1)–(4), is LAS if Rs > 1.
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
481
Proof. Evaluating the Jacobian of (1)–(4), at E1 gives 2 3 G 0 0 G l 1 6 G G a 0 7 6 7 1 K1 J ðE1 Þ ¼ 6 7; 4 K 3 0 5 0 K2 0 K4 0 l where, bP and K 4 ¼ cr; N from which it follows that the eigenvalues of J ðE1 Þ are the roots of the polynomial G1 ¼
b0 k4 þ b1 k3 þ b2 k2 þ b3 k þ b4 ¼ 0; with, b0 ¼ 1; aK 2 ; b1 ¼ K 3 þ 2l þ G þ K3 aK 2 b2 ¼ K 3 G þ 2l K 3 þ G þ þ cG þ l2 ; K3 aK 2 þ 2lK 3 G þ l2 K 3 þ cð2l þ raÞG ; b3 ¼ l2 G þ K3 b4 ¼ l2 K 3 G þ lcðl þ raÞG : Thus, the coefficients bi ði ¼ 0; 1; 2; 3; 4Þ are all positive if Rs > 1 (note that G** > 0 iff Rs > 1). Finally, it can be shown that b1 b2 b3 b23 b21 b4 ¼
ðf1 þ f2 þ f3 Þf4 > 0; K 33
where, f1 ¼ la2 K 22 þ l2 aK 2 K 3 þ 2G laK 2 K 3 þ 2laK 2 K 23 þ 2lK 33 G ; f2 ¼ lK 43 þ G l2 K 23 þ lK 23 G 2 þ l2 K 33 þ aK 2 K 23 G þ K 33 G2 ; f3 ¼ K 43 G þ acK 2 K 3 G þ cK 23 G fG þ ½l þ að1 rÞg; f4 ¼ craG K 3 þ 2l2 aK 2 þ 2l2 K 3 G þ 2l3 K 3 þ 2lK 23 G þ 2l2 K 23 þ 2clG K 3 : Hence, by the Routh–Hurwitz criterion, all eigenvalues of J ðE1 Þ have negative real parts whenever Rs > 1. Thus, E1 is LAS whenever it exists. h It is worth mentioning that for the special case r = 0, Qp ! 0 as t ! 1, giving a 3D limiting model (in P, S and Qt). This (reduced) model can further be simplified into a 2D model using the relation N = S + P + Qt; and the corresponding SPE can be shown to be globally-asymptotically stable (using, for instance, the Dulac criterion). For the full models, (1)–(4), although no proof for global stability is presented, numerous simulations of the model suggest that E1 is globally-asymptotically stable whenever Rs > 1 (as illustrated in Fig. 2b, d and f). Consequently, we offer the following conjecture. Conjecture 1. The unique smoking-present equilibrium of the models (1)–(4), is GAS in D n D0 , with D0 ¼ fðP ; S; Qt ; Qp Þ 2 D : S ¼ Qt ¼ 0g, whenever Rs > 1. 2.4. Numerical simulations of the basic model To illustrate the various theoretical results presented thus far, the basic models (1)–(4), is simulated using the parameter values/ranges in Table 1. These parameters are chosen arbitrarily, albeit plausibly, owing to the absence of available reliable data in the literature.
482
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499 1200
1000 900
1000
800
800
Smokers S(t)
Smokers S(t)
700 600 500 400
600
400
300 200
200
100 0 0
10
20
30
40
50
0 0
60
50
100
500
500
450
450
400
400
350 300 250 200 150
350 300 250 200 150
100
100
50
50
0 0
5
10
15
20
25
0 0
30
10
20
30
40
50
60
70
80
Time (years)
1200
1200
1000
1000
Total smokers s(t) + Qt(t)
Total smokers s(t) + Qt(t)
Time (years)
800
600
400
800
600
400
200
200
0 0
150
Time (years)
Temporary smokers Qt(t)
Temporary Quitters Qt(t)
Time (years)
10
20
30
Time (years)
40
50
60
0 0
20
40
60
80
100
120
140
Time (years)
Fig. 2. Time series plot of the basic models (1)–(4), with different initial conditions. Parameter as in Table 1. (a) S(t) for Rs < 1; (b) S(t) for Rs > 1; (c) Qt(t) for Rs < 1; (d) Qt(t) for Rs > 1; (e) S(t) + Qt(t) for Rs < 1; (f) S(t) + Qt(t) for Rs > 1.
Table 2 depicts the variables of the model at steady-state as a function of b and Rs . It is clear from this table that the number of smokers at steady-state increases with increasing Rs > 1. It also shows a marked decrease in the number of potential smokers as Rs increases. Further, in line with Theorems 1 and 3 and Conjecture 1,
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
483
Table 1 Description of parameters for the basic models (1)–(4) Parameter
Description
Estimated value/range
Nominal value
1/l a
Average duration of smoking activity Rate at which temporary quitters Revert back to smoking Rate of quitting smoking Fraction of smokers who permanently quit Contact rate
20 years
20
c r b
a 2 [0, 1) c 2 [0, 1) r 2 ð0; 1Þ 0.3 per year
3 0.3 0.5 0.3
Table 2 Effect of b on Rs and P ; S; Qt ; Qp at steady-state using the basic models (1)–(4) b
Rs
P
S
Qt
Qp
Stable equilibrium
0.003 0.03 0.02 0.01 0.1 0.2 0.3 0.5 1 2 3 5
0.0118 0.1184 0.0789 0.0394 0.3948 0.7896 1.1844 1.9741 3.9482 7.8964 11.8446 19.7411
1600 1600 1600 1600 1600 1600 1351 811 405 203 135 81
0 0 0 0 0 0 49 156 236 276 289 300
0 0 0 0 0 0 3 10 16 18 19 20
0 0 0 0 0 0 197 623 943 1103 1157 1199
E0 E0 E0 E0 E0 E0 E1 E1 E1 E1 E1 E1
this table shows that E0 is always stable whenever Rs 6 1 and E1 is stable if Rs > 1 (see also Fig. 2 for simulations using various initial conditions). This table further shows a marked decrease in the population of potential smokers as Rs increases. The combined effect of the rate of reversion to smoking (a) and quitting smoking temporarily (c) on Rs is investigated using various values of r. The results, depicted in Fig. 3, show a decrease in Rs as r increases. In particular, for r = 0.5 (Fig. 3), almost all of the Rs contours are less than unity, signifying that this value of r (and other fixed parameters) will lead to the elimination of smoking from the community (in line with Theorem 1). It is worth noting that this result is further supported by carrying out a sensitivity analysis on Rs (i.e., differentiating Rs partially with respect to r), giving oRs barðl þ aÞ ¼ < 0; 2 or ½lðl þ aÞ þ cðl þ raÞ from which it follows that Rs is a decreasing function of r. Fig. 3 further shows that Rs increases with increasing a. Here, too, sensitivity analysis on Rs shows oRs blcð1 rÞ ¼ > 0: 2 oa ½lðl þ aÞ þ cðl þ raÞ It is further evident from Fig. 3 that increasing the quitting rate, c, decreases Rs ; which is in line with the calculation below: oRs bðl þ aÞðl þ arÞ ¼ < 0: oc ½lðl þ aÞ þ cðl þ raÞ2 3. Extension of the basic model Although the basic models (1)–(4), allows for the assessment of the dynamics of smoking in a closed population, it does not incorporate some of the other important aspects (and consequences) associated with
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
2
2 3
Number of days before quitting temporarily(1/γ)
4
1.8 1.6 3
4
1
2
1.4 1.2 3
1
2
0.8 1
2
0.6 1
0.4 1
0.2 1000
2000
3000
4000
5000
6000
7000
8000
2.5
1.8 1.6 1.4
1.5
2
1.2
1 5 0.
1.5
1
1
0.8 1
0.5
0.6 0.5
0.4 0.5
0.2
9000 10000
1000
Number of days before returning back to smoking(1/α)
1.6
1.8
1.4
1.2
1.4
1.2 1
1.4 1.2
0.8
1
1.6
0.8
0.6
1.2 1
0.8
0.6
1
0.4
0.8
0.8
0.6
0.4
0.6
0.6
0.4
0.4
0.2
0.4 0.2
0.2
0.2 1000
2000
3000
4000
3000
4000
5000
6000
7000
8000
9000 10000
2
1.4
1.8 1.6
2000
Number of days before returning back to smoking(1/α)
Number of days before quitting temporarily(1/γ)
Number of days before quitting temporarily(1/γ)
2
1.5
2
1
5
Number of days before quitting temporarily(1/γ)
2
1
484
5000
6000
7000
8000
1.2
0.8
1
1.8 1.6
1.2
1
0.8 0.6
1.4 1
0.8
1.2
0.6 0.8
1 0.8
0.6 0.4
0.6
0.4
0.6
0.4
0.2
0.4
0.2 0.2
0.2
9000 10000
1000
Number of days before returning back to smoking(1/α)
2000
3000
4000
5000
6000
7000
8000
9000 10000
Number of days before returning back to smoking(1/α)
2
Number of days before quitting temporarily(1/γ)
1 0.8
1.8 1 0.8
1.6
0.6
1.4
0.8
0.6
1.2
0.6
1 0.4
0.8
0.4 0.4
0.6 0.2
0.4
0.2 0.2
0.2 1000
2000
3000
4000
5000
6000
7000
8000
9000 10000
Number of days before returning back to smoking(1/α)
Fig. 3. Contour plots of the basic models (1)–(4), depicting Rs as a function of 1c and 1a. Parameters as in Table 1 with various values of r. (a) r = 0.1; (b) r = 0.1; (c) r = 0.1; (d) r = 0.1; (e) r = 0.1.
tobacco smoking. One of such is the degree, or frequency, of smoking, where some smokers do so in a mild manner (categorized as mild smokers) while others may smoke more frequently (defined as chain smokers) per
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
485
unit time. Further, the transition dynamics (back-and-forth movement) between these subclasses, of varying smoking frequencies, also needs to be incorporated into the model. Another key aspect lacking in the basic model is the public health consequence of smoking. In other words, since millions of people die of smoking-related illnesses every year [17], it is necessary to include such fact in a population model of smoking dynamics. To account for the aforementioned drawbacks, the basic model is now extended to include two classes of smokers – a mild smoking class (S1(t)) and a chain smoking class (S2(t)) – and a class of smoking-related illnesses (C(t)). The extended model is given by: dP bP ðS 1 þ wS 2 Þ ¼P lP ; dt N dS 1 ð1 hÞbP ðS 1 þ wS 2 Þ ðg1 þ l þ c1 þ sÞS 1 þ g2 S 2 þ a1 Qt ; ¼ N dt dS 2 hbP ðS 1 þ wS 2 Þ ðg2 þ l þ c2 þ snÞS 2 þ g1 S 1 þ a2 Qt ; ¼ N dt dQt ¼ c1 ð1 r1 ÞS 1 þ c2 ð1 r2 ÞS 2 ðl þ a1 þ a2 ÞQt ; dt dQp ¼ c1 r1 S 1 þ c2 r2 S 2 lQp ; dt dC ¼ sS 1 þ snS 2 lC dC; dt
ð13Þ
where, now, the total population is N(t) = P(t) + S1(t) + S2(t) + Qt(t) + Qp(t) + C(t). In (13), a proportion, 1 h (with 0 < h < 1), of new smokers are assumed to be mild smokers, and the remaining fraction, h, are chain smokers. It is further assumed that chain smokers have a higher probability of generating more new smokers (by a factor w P 1) relative to the mild smokers. This may be true due to a possible higher social influence the chain smokers may have (amongst peers in the smoking environment) over mild smokers, coupled with the possibility that their high smoking frequency may tempt people around them (who are probably passive smokers anyway) to become smokers. Mild smokers become chain smokers (at a rate g1) and the reverse process occurs at a rate g2. While mild smokers develop smoking-related illnesses at a rate s, those with chain smoking habit are assumed to do so at an increased rate, sn (where n > 1). That is, it is assumed that those in the chain smoking class have a higher probability of developing smoking-related illnesses in comparison to those in the mild smoking class. Individuals with such illnesses die at an increased mortality rate, d. A fraction, 1 ri (i = 1, 2), of individuals in the smoking classes S1 and S2 temporarily quit smoking, at rates c1 and c2, respectively; whilst the remaining fraction, ri, quit smoking permanently. Individuals who temporarily quit smoking revert to mild smoking habit at a rate a1, and to chain smoking habit at a rate a2. Finally, individuals in all classes suffer a natural death at a rate l. A schematic description of the extended model (13) is given in Fig. 1b. 3.1. Basic properties of the extended model Here, too, it is plausible to assume that all the state variables and parameters of the extended model (13) are non-negative. Consider the biologically-feasible region M ¼ fðP ; S 1 ; S 2 ; Qt ; Qp ; CÞ 2 R6þ : P þ S 1 þ S 2 þ Qt þ Qp þ C 6 P=lg: The rate of change of the total population, obtained by adding all the equations in the extended model (13), is given by dN ¼ P lN dC: dt
ð14Þ
It is worth noting here that, unlike in the case of the basic models (1)–(4), the total population N(t) is not constant (except if d = 0). Since the right hand side of (14) is bounded by P lN, a standard comparison theorem can be used to show that N ðtÞ 6 N ð0Þelt þ Pl ð1 elt Þ. In particular, N ðtÞ 6 Pl if N ð0Þ 6 Pl. Thus,
486
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
M is positively-invariant. Hence, it is sufficient to consider the dynamics of the flow generated by (13) in M. It is easy to see, by comparison theorem, that lim inft ! 1N(t) 6 P/l. Thus, the omega limit sets of all solutions in R6þ are contained in M. That is, solutions in M remain in M for all time, and those outside M (but in R6þ ) are eventually attracted to M. 3.2. Stability of the SFE of the Extended Model 3.2.1. Local stability The extended model (13) has a smoking-free equilibrium, given by X0 ¼ ðP ; S 1 ; S 2 ; Qt ; Qp ; C Þ ¼ ðP=l; 0; 0; 0; 0; 0Þ: Here, now, the next generation matrices, F and V, are given 1 0 0 J1 ð1 hÞb ð1 hÞwb 0 0 C B B hwb 0 0C B hb B g C and V ¼ B 1 F ¼B B C B J 0 0 0 0A @ @ 3 0
0
0 0
s
by
1
g2
a1
J2
a2
J4
J5
C 0C C; 0C A
ns
0
J8
0
where, J 1 ¼ g1 þ l þ c1 þ s; J 4 ¼ c2 ð1 r2 Þ;
J 2 ¼ g2 þ l þ c2 þ sn;
J 5 ¼ l þ a1 þ a2 ;
J 3 ¼ c1 ð1 r1 Þ;
J 6 ¼ c1 r1 ;
J 7 ¼ c 2 r2 ;
J 8 ¼ l þ d:
It follows that, using the next generation approach, the smokers generation number of the extended model (13), denoted by Rsd ¼ qðFV 1 Þ; is Rsd ¼
bðD1 þ D2 Þ ; D3
with, D1 ¼ h½g2 J 5 þ a1 J 4 þ wðJ 1 J 5 a1 J 3 Þ; D2 ¼ ð1 hÞ½J 2 J 5 a2 J 4 þ wða2 J 3 þ g1 J 5 Þ; D3 ¼ J 1 J 2 J 5 g2 a2 J 3 a1 J 2 J 3 g1 g2 J 5 g1 a1 J 4 a2 J 1 J 4 : It is shown (see Appendix C) that D1 ; D2 ; and D3 are all positive, so that Rsd > 0: Thus, the following result is established by Theorem 2 of [28]. Lemma 2. The SFE, X0 , of the extended model (13) is LAS whenever Rsd < 1, and unstable if Rsd > 1. The quantity Rsd measures the average number of new smokers generated by a single smoker in a population of completely susceptible (potential smokers) individuals. If this number (Rsd ) is less than unity, then a small influx of smokers will not generate large number (‘‘outbreaks’’) of smokers. If Rsd > 1, on the other hand, such influx would result in large number of new smokers. It is easy to see that when g1 ¼ g2 ¼ c2 ¼ a2 ¼ s ¼ n ¼ h ¼ w ¼ d ¼ r2 ¼ 0, and c1 ¼ c; a1 ¼ a; r1 ¼ r, then Rsd ¼ Rs (i.e., the extended model (13) reduces to the basic models (1)–(4), in this case). A global stability proof for X0 is given below. 3.2.2. Global stability Theorem 4. The SFE of the extended model (13) is GAS in M if Rsd < 1. Proof. Consider the Lyapunov function: L ¼ l1 S 1 þ l 2 S 2 þ l 3 Q t ; where,
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
487
l1 ¼ J 2 J 5 a2 J 4 þ wða2 J 3 þ g1 J 5 Þ; l2 ¼ g2 J 5 þ a1 J 4 þ wðJ 1 J 5 a1 J 3 Þ; l3 ¼ a1 ðwg1 þ J 2 Þ þ a2 ðg2 þ wJ 1 Þ: Thus,
_ ¼ l1 ð1 hÞbP ðS 1 þ wS 2 Þ J 1 S 1 þ g2 S 2 þ a1 Qt L N hbP ðS 1 þ wS 2 Þ J 2 S 2 þ g1 S 1 þ a2 Qt þ l3 ðJ 3 S 1 þ J 4 S 2 J 5 Qt Þ þ l2 N S 1 þ wS 2 ¼ ½bP ðD1 þ D2 Þ ND3 N bP ðD1 þ D2 Þ ¼ D3 ðS 1 þ wS 2 Þ 1 ND3 bðD1 þ D2 Þ 6 D3 ðS 1 þ wS 2 Þ 1 since P 6 N inM D3 ¼ D3 ðS 1 þ wS 2 ÞðRsd 1Þ < 0
for Rsd < 1:
Hence, using similar argument as in the proof of Theorem 1, the SFE of the extended model, X0 , is globally-asymptotically stable whenever Rsd < 1. h An alternative proof for Theorem 4, using comparison theorem, is given in Appendix D. 3.3. SPE of the extended model 3.3.1. Existence of smoking-present equilibrium Using the method described in Section 2.3, with G ¼
bðS 1 þ wS 2 Þ ; N
ð15Þ
it follows that the extended model has a unique SPE, denoted by X1 ¼ ðP ; S 1 ; S 2 ; Qt ; Qp ; C Þ; where,
P ¼ Q t
P ; G þl
S 1 ¼
PG H 2 ; ¼ D3 ðG þ lÞ
PG H 0 ; D3 ðG þ lÞ
Q p
S 2 ¼
PG H 1 ; D3 ðG þ lÞ
PG ðH 11 þ H 12 Þ ; ¼ D3 lðG þ lÞ
C
PG sðH 21 þ H 22 Þ ¼ D3 J 8 ðG þ lÞ
ð16Þ
with, H 0 ¼ ½ðJ 5 J 2 J 4 a2 Þð1 hÞ þ ða1 J 4 þ J 5 g2 Þh; H 1 ¼ ½ðg1 J 5 þ a2 J 3 Þð1 hÞ þ ðJ 5 J 1 J 3 a1 Þh; H 2 ¼ ½ðJ 2 J 3 þ g1 J 4 Þð1 hÞ þ ðJ 1 J 4 þ g2 J 3 Þh; H 11 ¼ ½ðJ 4 a1 þ J 5 g2 ÞJ 6 h þ ðg1 J 5 þ J 3 a2 ÞJ 7 ð1 hÞ; H 12 ¼ ½ðJ 5 J 1 J 3 a1 ÞJ 7 h þ ðJ 5 J 2 J 4 a2 ÞJ 6 ð1 hÞ; H 21 ¼ ½ðJ 5 g2 þ a1 J 4 Þh þ ðg1 J 5 þ J 3 a2 Þnð1 hÞ; H 22 ¼ ½ðJ 5 J 1 J 3 a1 Þnh þ ðJ 5 J 2 J 4 a2 Þð1 hÞ: Using (16) in (15), and simplifying, gives bG D3
G ¼
1þ
G D3
½H 0 þ wH 1
22 Þ 12 ½H 0 þ H 1 þ H 2 þ H 11 þH þ sðH 21JþH l 8
:
ð17Þ
488
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
As before, solving for G** from (17) and substituting the result into (16) gives the corresponding SPE of the model (13). Here, too, G** = 0 corresponds to the smoking-free equilibrium, X0 . For G** 5 0 (i.e., we are looking for an SPE of (13)), (17) can be simplified to c0 G c1 ¼ 0; where, c0 ¼
ð18Þ
1 H 11 þ H 12 sðH 21 þ H 22 Þ þ H0 þ H1 þ H2 þ and D3 J8 l
c1 ¼ Rsd 1:
From (18), G ¼ cc10 . It follows that since the terms in c0 are all positive, then the coefficients c0 and c1 are positive for Rsd > 1. Thus, the extended model (13) has a unique SPE whenever Rsd > 1. It is further evident that if Rsd < 1, then the coefficient c1 < 0, so that G** < 0 (which is biologically meaningless). Here, too, G** = 0 whenever Rsd ¼ 1. Thus, the SPE ðX1 Þ does not exist whenever Rsd 6 1. These results are summarized below. Theorem 5. The extended model (13) has a unique positive smoking-present equilibrium if Rsd > 1, and no positive smoking-present equilibrium if Rsd 6 1. 3.3.2. Local stability of the SPE of the extended model We claim the following result. Theorem 6. The unique smoking-present equilibrium of the extended model (13) is LAS whenever Rsd > 1 and is close to 1. Proof. Linearizing the extended model (13) around its unique SPE is laborious (and not really tractable) owing to its high dimensionality. Consequently, an alternative approach is considered. It is based on the use of the centre manifold theory [5], as described in [7, (Theorem 4.1)]. To apply this method, the following simplification and change of variables are made first of all. Let P ¼ x1 ; S 1 ¼ x2 ; S 2 ¼ x3 ; Qt ¼ x4 ; Qp ¼ x5 , and C = x6, so that N ¼ x1 þ x2 þ x3 þ x4 þ x5 þ x6 . Further, by using vector notation X ¼ ðx1 ; x2 ; x3 ; x4 ; x5 ; x6 ÞT , T ¼ ðf1 ; f2 ; f3 ; f4 ; f5 ; f6 Þ , as follows: the extended model (13) can be written in the form dX dt dx1 dt dx2 dt dx3 dt dx4 dt dx5 dt dx6 dt
bx1 ðx2 þ wx3 Þ lx1 ; x1 þ x2 þ x3 þ x4 þ x5 þ x6 bð1 hÞx1 ðx2 þ wx3 Þ ¼ f2 ¼ J 1 x2 þ g1 x3 þ a1 x4 ; x1 þ x2 þ x3 þ x4 þ x5 þ x6 bhx1 ðx2 þ wx3 Þ ¼ f3 ¼ J 2 x3 þ g1 x2 þ a2 x4 ; x1 þ x2 þ x3 þ x4 þ x5 þ x6 ¼ f1 ¼ P
ð19Þ
¼ f4 ¼ J 3 x2 þ J 4 x3 J 5 x4 ; ¼ f5 ¼ J 6 x2 þ J 7 x3 lx5 ; ¼ f6 ¼ sx2 þ snx3 J 8 x6 :
The Jacobian of the system (19) at the SFE, X0 , is 0 l b bw B 0 ð1 hÞb J ð1 hÞbw þ g1 B 1 B B 0 hb þ g1 hbw J 2 B J ðX0 Þ ¼ B B 0 J4 J3 B B 0 J6 J7 @ 0 s sn
given by 0
0
a1
0
a2
0
J 5
0
0
l
0
0
0
1
C 0 C C 0 C C C; 0 C C 0 C A J 8
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
489
from which it can also be shown (as before) that Rsd ¼
bðD1 þ D2 Þ : D3
Consider the case when Rsd ¼ 1. Suppose, further, that b = b* is chosen as a bifurcation parameter. Solving for b from Rsd ¼ 1 gives b ¼ b ¼
D3 : D1 þ D2
The right eigenvectors of J(b*), the Jacobian of (19) evaluated at b = b*, are given by w ¼ ½w1 ; w2 ; w3 ; w4 ; w5 ; w6 T , where b w2 b ww4 ; w2 ¼ w2 > 0; w3 ¼ w3 > 0; l J 3 w2 þ J 4 w3 J 6 w2 þ J 7 w3 sw2 þ snw3 ; w6 ¼ ; w5 ¼ : w4 ¼ J5 l J8
w1 ¼
Further, ðJ b Þ has a left eigenvector v ¼ ½v1 ; v2 ; v3 ; v4 ; v5 ; v6 ; where v1 ¼ 0;
v2 ¼ v2 > 0;
v3 ¼ v3 > 0;
v4 ¼
a1 v 2 þ a2 v 3 ; J5
v5 ¼ 0;
v6 ¼ 0:
Note that the transformed system (19), with b = b*, has at least one non-hyperbolic equilibrium point (i.e., the linearized system has at least one eigenvalue with zero real part). Hence, the center manifold theory can be used to analyse the dynamics of (19), with b = b*. In particular, the following theorem (Theorem 4.1 in [7]) will be used to show the LAS of the SPE of (19) (which is the same as the SPE of the original system 13). The theorem in [7] is reproduced below for convenience. Theorem 7 (Castillo-Chavez and Song [7]). Consider the following general system of ordinary differential equations with a parameter / dx ¼ f ðx; /Þ; dt
f : Rn R ! R
and
f 2 C2 ðRn RÞ;
where 0 is an equilibrium point of the system (i.e., f(0, /) 0 for all /) and assume ofi ð0; 0ÞÞ is the linearization matrix of the system around the equilibrium 0 with / evaluA1: A ¼ Dx f ð0; 0Þ ¼ ðox j ated at 0. Zero is a simple eigenvalues of A and other eigenvalues of A have negative real parts; A2: Matrix A has a right eigenvectors w and a left eigenvector v corresponding to the zero eigenvalue.
Let fk be the kth component of f and a¼
n X
vk wi wj
k;i;j¼1
b¼
n X k;i¼1
vk wi
o2 fk ð0; 0Þ; oxi oxj
o2 fk ð0; 0Þ: oxi o/
The local dynamics of the system around 0 is totally by a and b. (i) a > 0; b > 0. When / < 0 with j/j 1; 0 is locally asymptotically stable and there exists a positive unstable equilibrium; when 0 < / 1, 0 is locally asymptotically stable equilibrium; (ii) a < 0; b < 0. When / < 0 with j/j 1; 0 is unstable; when 0 < / 1, 0 is locally asymptotically stable equilibrium, and there exists a positive unstable equilibrium; (iii) a > 0; b < 0. When / < 0 with j/j 1; 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0 < / 1, 0 is unstable, and a positive unstable equilibrium appears;
490
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
(iv) a < 0; b > 0. When / changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable. To apply the above theorem, the following computations are necessary. Computations of a and b For the system (19), the associated non-zero partial derivatives of the right hand side functions (fi) are given by o2 f2 ð1 hÞð1 þ wÞb l o2 f3 hð1 þ wÞb l ; ; ¼ ¼ P P ox2 ox3 ox2 ox3 o2 f2 o2 f2 o2 f2 o2 f2 ð1 hÞb l ; ¼ 2 ¼ ¼ ¼ P ox2 ox4 ox2 ox2 ox5 ox2 ox6 o2 f2 o2 f2 o2 f2 o2 f2 ð1 hÞb wl ; ¼ 2 ¼ ¼ ¼ P ox3 ox4 ox3 ox3 ox5 ox3 ox6 o2 f3 o2 f3 o2 f3 o2 f3 hb l ; ¼ 2 ¼ ¼ ¼ P ox2 ox4 ox2 ox2 ox5 ox2 ox6 o2 f3 o2 f3 o2 f3 o2 f3 hb wl ; ¼ 2 ¼ ¼ ¼ P ox3 ox4 ox3 ox3 ox5 ox3 ox6 o2 f 2 ¼ 1 h; ox2 ob
o2 f 2 ¼ ð1 hÞw; ox3 ob
o2 f3 ¼ h; ox2 ob
o2 f3 ¼ hw: ox3 ob
Thus, it follows from the above expressions that a¼
6 X
vk wi wj
k;i;j¼1
¼
6 X
v2 wi wj
i;j¼1
o2 f k ð0; 0Þ oxi oxj
6 X o2 f2 o2 f3 þ v3 wi wj oxi oxj i;j¼1 oxi oxj
2b lðw2 þ w3 þ w4 þ w5 þ w6 Þðw2 þ ww3 Þðv2 ð1 hÞ þ v3 hÞ < 0; ¼ P and b¼
6 X
vk w i
o2 fk ð0; 0Þ oxi o/
v2 wi
6 X o2 f2 o2 f3 þ v3 w i oxi o/ i¼1 oxi o/
k;i¼1
¼
6 X i¼1
¼ ½v2 ð1 hÞ þ v3 hðw2 þ ww3 Þ > 0: Thus, a < 0 and b > 0. Hence, by Theorem 7(iv) above, the unique SPE of the extended model, which exists whenever Rsd > 1 (Theorem 6), is LAS whenever Rsd > 1 and b* < b with b close to b*. h It is worth mentioning that a similar approach is used in [22] to establish the local stability of endemic equilibrium. Here, too, extensive numerical simulations (see Fig. 4b, d and f) suggest the following conjecture: Conjecture 2. The unique smoking equilibrium of model 13 is globally-asymptotically stable in M n M0 , with M0 ¼ fðP ; S 1 ; S 2 ; Qt ; Qp ; CÞ 2 M : S 1 ¼ S 2 ¼ Qt ¼ 0g, whenever Rsd > 1. Overall, the qualitative analyses in Sections 2 and 3 show that both the basic models, (1)–(4), and the extended model (13) exhibit the same dynamical behaviour (each has a globally-stable SFE whenever its associated smoking-generation number is less than unity; and has a unique LAS SPE whenever this quantity
800
800
700
700
600
600
Mild Smokers S1(t)
Mild Smokers S1(t)
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
500 400 300
500 400 300
200
200
100
100
0
0
2
4
6
8
10
12
14
0 0
16
2
4
6
500
500
450
450
400
400
350
350
300 250 200 150
0 10
12
14
0 0
16
2
4
6
500
500
450
450
400
400
350 300 250 200 150
8
16
10
12
14
16
200 150
50
0
Time(years)
14
250
100
6
12
300
50
4
10
350
100
2
8
Time(years)
Temporary Quitters Qt(t)
Temporary Quitters Qt(t)
Time(years)
0
16
150
50
8
14
200
100
6
12
250
50
4
10
300
100
2
8
Time(years)
Chain Smokers S2(t)
Chain Smokers S2(t)
Time(years)
0
491
10
12
14
16
0 0
2
4
6
8
Time(years)
Fig. 4. Time series plot of the extended model (13) with different initial conditions. Parameter as in Table 3. (a) S1(t) for Rsd < 1; (b) S1(t) for Rsd > 1; (c) S2(t) for Rsd < 1; (d) S2(t) for Rsd > 1; (e) Qt(t) for Rsd < 1; (f) Qt(t) for Rsd > 1.
exceeds unity). In other words, adding the compartments (and dynamics) of the mild and chain smokers, as well as the class of smoking-related illnesses, to the basic model does not alter its dynamical features. The public health implication of these results is that the number of smokers in a community, together with the
492
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
associated smoking-related mortality, can be effectively controlled if the smoking-generation threshold is brought down to a value less than unity. On the other hand, if this threshold exceeds unity, then the number of smokers and smoking-related mortality will rise. 3.4. Numerical simulations of the extended model The extended model (13) is simulated using the parameter values in Table 3. As evident from Table 4, the solution profiles converge to X0 whenever Rsd < 1, and to X1 for Rsd > 1 (in line with Theorems 4 and 6). In addition to showing increasing number of smokers and smoking-related mortality (at steady-state) for increasing Rsd values, this table further shows that large values of the effective contact rate, b, could lead to drastic reductions in the entire population (including the near-complete absence of potential smokers). Time series plots for cases where Rsd < 1 and Rsd > 1 are depicted in Fig. 4, further supporting Theorems 4 and 6. Simulations for cumulative mortality and mortality per unit time are also carried out. Fig. 5a shows the cumulative death resulting from the illnesses for Rsd < 1; suggesting that about 38% of the total population would have acquired smoking habits at steady-state (which is attained in about six years time). Fig. 5b and c depict the time series plots of the mortality per unit time for the case Rsd < 1 and Rsd > 1, respectively, showing a peak mortality of approximately 230 (for Fig. 5b and 260 (Fig. 5c) people. These figures show reductions in smoking-related mortality with decreasing values of Rsd . Table 3 Description of parameters for the extended model (13) Parameter
Description
Estimated value/range
Nominal value
P 1/l a1 a2 c1 c2 r1 r2 g1 g2 b s n d w h
Recruitment rate into the potential smokers class Average duration of smoking activity Rate at which temporary quitters revert back to S1 class Rate at which temporary quitters revert back to S2 class Quitting rate for mild smokers Quitting rate for chain smokers Fraction of quitters from S1 class that permanently quit smoking Fraction of quitters from S2 class that permanently quit smoking Rate at which mild smokers become chain smokers Rate at which chain smokers become mild smokers Contact rate (b = cq) Rate at which mild smokers develop smoking-related illnesses Modification parameter Smoking-related mortality rate Modification parameter Fraction of new smokers that begin chain smoking
30 per year 20 years a1 2 (0, 1) a2 2 (0, 1) c1 2 (0, 1) c2 2 (0, 1) r1 2 ð0; 1Þ r2 2 ð0; 1Þ g1 2 ð0; 1Þ g2 2 ð0; 1Þ b 2 ð0; 1Þ s 2 ð0; 1Þ n>1 d 2 ð0; 1Þ wP1 h 2 ð0; 1Þ
30 20 0.3 0.3 0.4 0.5 0.4 0.5 2 3 0.2 1 3 9 9 0.6
Table 4 Effect of b on Rsd and P ; S 1 ; S 2 ; Qt ; Qp ; C at steady-state using model (13) b
Rsd
P
S1
S2
Qt
Qp
C
Stable equilibrium
0.003 0.03 0.02 0.1 0.2 0.3 0.5 1 2 3 5 7
0.0061 0.0605 0.0403 0.2017 0.4034 0.6051 1.0085 2.0169 4.0338 6.0507 10.084 14.1183
600 600 600 600 600 600 571 86 32 20 11 8
0 0 0 0 0 0 0 8 8 9 9 9
0 0 0 0 0 0 0 5 5 6 6 6
0 0 0 0 0 0 0 5 5 5 5 5
0 0 0 0 0 0 3 49 54 55 56 57
0 0 0 0 0 0 1 21 24 24 24 25
X0 X0 X0 X0 X0 X0 X1 X1 X1 X1 X1 X1
250
250
200
200
150
150
Mortality δC
Cummulative Death
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
100
50
0 0
493
100
50
2
4
6
8
10
12
14
0 0
16
2
4
6
Time(years)
8
10
12
14
16
Time(years)
300
250
Mortality δC
200
150
100
50
0
0
2
4
6
8
10
12
14
16
Time(years)
Fig. 5. Simulations of the extended model (13). (a) Cumulative mortality; (b) total death for Rsd < 1; (c) total death for Rsd < 1.
In order to ascertain the impact of the relative sojourn (duration) in the two smoking classes on the overall smoking dynamics, more simulations are carried out for a 3 year time period as follows. In Fig. 6, contours of Rsd are plotted as a function of number of days spent in the S1 class ðc11 Þ and the S2 class ðc12 Þ. Based on the parameter values used in these simulations, this figure shows that if chain smokers do not spend more than 1.5 years in this class, the number of smokers in the community (at steady state) will be effectively controlled or eliminated (since Rsd < 1 in this case) regardless of the time spent by smokers in the mild smoking class (S1). Similarly, if mild smokers practice their mild smoking habit for less than 1.5 years, then the number of smokers in the community (at steady state) will be greatly reduced (or eliminated) from the population irrespective of the dynamics in the chain smoking class. Furthermore, the number of smokers in the community (at steady state) will persist if individuals in either classes remain there for a duration longer than 1.5 years (here Rsd > 1). Of course, it should be emphasized that these quantitative estimates are dependent on the parameter values chosen, and, owing to expected uncertainties in these parameters, the results should be taken only as a general guide. The impact of the transition from mild to chain smoking class is also monitored by plotting Rsd contours as a function of the number of days a smoker of class i spends in class j (i 5 j). Fig. 7a shows a marked decrease in Rsd for longer duration in the mild smoking class by individuals in the chain smoking class. In other words, the number of smokers will be greatly reduced if chain smokers revert to mild smoking, and stay in the mild smoking class for a long duration. On the other hand, Rsd increases if mild smokers move to the chain
494
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
0.9
0.8
0.7 3
1
2.5
2 2
1/γ
0.9
0.8
0.7
1
1.5 1
1
0.9 0.
8
0.7
0.5
0.9
0.8
0.8
07
0.5
1
1.5
2
1/γ
2.5
3
1
Fig. 6. Contour plots of the extended model (13) depicting Rsd as a function of
1 c1
and c12 . Parameters as in Table 3.
smoking class and stay there for long duration (Fig. 7b). Thus, it is more beneficial (in terms of public health) for chain smokers to acquire mild smoking habits for extended duration, than mild smokers becoming chain smokers for a long periods. 3.5. Conclusions A four-dimensional compartmental deterministic model is designed and used to monitor smoking dynamics in a population. This model was extended to incorporate the heterogeneity of smoking frequency and smoking-related illnesses. The models were rigorously analysed to gain insight into their dynamical features. The study shows the following: (i) Each of the two models considered in this study has a globally-stable smoking-free equilibrium if a certain threshold quantity, known as the smokers-generation number, is less than unity; indicating that the number of smokers in the community will be brought to zero if public health measures that make (and keep) the threshold to a value less than unity are carried out; (ii) For the basic model, the number of smokers can be reduced by reducing the contact rate (b), reversion rate (a) and increasing quitting rates (c) and the fraction (r) of temporary quitters who permanently giveup smoking; (iii) For the extended model, the number of smokers can similarly be reduced by reducing the contact rates (b and bw), transition rates, (g1 and g2), reversion rates (a1 and a2), and increasing the rate at which smokers temporarily quit smoking (c1 as well as the c2) and fraction of permanent quitters (r1 and r2); (iv) The transition from chain to the mild smoking class plays a crucial role in the smoking dynamics. More new smokers are generated if mild smokers adopt a chain smoking habit; and markedly reduced number of new smokers are recorded if chain smokers move to the mild smoking class; (v) If mild and chain smokers spend less than 1.5 years in their respective smoking classes, then the number of smokers in the community will be effectively controlled (or eliminated). The number of smokers will persist otherwise. Similarly, if chain smokers spend less than 1.5 years in this class, the number of smokers will be effectively controlled (as the same reason applies).
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
495
1.3
1.25
1.2
1.15
Rsd
1.1
1.05
1
0.95
0.9
0.85
0.8 0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
Duration in mild smoking class by a chain smoker: (η /(η + μ + γ + τ)) 2
1
1
1.4
1.35
1.3
1.25
R
sd
1.2
1.15
1.1
1.05
1
0.95
0.9 0
0.5
1
1.5
2
2.5
Duration in chain smoking class by a mild smoker: (η /(η + μ + γ + τξ)) 1
2
2
Fig. 7. Simulations of the extended model (13). (a) Plot of the smokers generation number Rsd as a function of the number of days (years) mild smokers spend in the chain smoking class S2; (b) Plot of the smokers generation number Rsd as a function of the number of days (years) chain smokers spend in the mild smoking class S1.
In summary, although the models considered here can further be refined to incorporate other smokingrelated phenomena, such as the impact of second-hand (passive) smoking, recruitment of smokers from external sources into the population (so that a fraction of recruited people, P, go directly into the S1 and S2 classes, etc), age structure, smoking advertising and mass consumption, this study shows that smoking and smoking-related illnesses can be effectively controlled in a community if public health measures that can bring the smokers generation number to a value less than unity are implemented.
496
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
Acknowledgements One of the authors (ABG) acknowledges, with thanks, the support in part of the Natural Science and Engineering Research Council (NSERC) and Mathematics of Information Technology and Complex Systems (MITACS) of Canada. Appendix A. Next generation method [28] The next generation method is used to establish the local asymptotic stability of the SFE [11,28]. The method was first introduced by Diekmann and Hesterbeek [11], and refined for epidemiological models by van den Driessche and Watmough [28]. The formulation in [28], which is based on disease transmission model, is reproduced below (note that this method can be implemented on the smoking models in this paper by essentially replacing the word ‘‘disease’’ with ‘‘smoking’’). Suppose the given disease transmission model, with nonnegative initial conditions, can be written in terms of the following system: x_ i ¼ fi ðxÞ ¼ Fi ðxÞ Vi ðxÞ;
i ¼ 1; . . . ; n;
ð20Þ
þ where Vi ¼ V and the function satisfy the five axioms below. First of all, X s ¼ fx P 0jxi ¼ i Vi 0; i ¼ 1; . . . ; mg with m < n is defined as the disease-free states (non-infected state variables of the model) of the model, where x ¼ ðx1 ; . . . ; xn Þt ; xi P 0 represents the number of individuals in each compartment of the model.
(A1) (A2) (A3) (A4) (A5)
if x P 0, then Fi ; Vþ i ; Vi P 0 for i ¼ 1; . . . ; m. if xi = 0, then Vi ¼ 0. In particular, if x 2 Xs then V i ¼ 0 for i ¼ 1; . . . ; m. Fi ¼ 0 if i > m. if x 2 Xs, then Fi ðxÞ ¼ 0 and Vþ i ðxÞ ¼ 0 for i ¼ 1; . . . ; m. If FðxÞ is set to zero, then all eigenvalues of Df(x0) have negative real parts.
Here, Fi ðxÞ represents the rate of appearance of new infections in compartment i; Vþ i ðxÞ represents the rate of transfer of individuals into compartment i by all other means, and V i ðxÞ represents the rate of transfer of individuals out of compartment i. It is assumed that these functions are at least twice continuously differentiable in each variable [28]. Lemma 3 (van den Driessche and Watmough [28]). If x is a DFE of 20 and fi(x) satisfy (A1)–(A5), then the derivatives DFðxÞ and DVðxÞ are partitioned as F 0 V 0 DFðxÞ ¼ ; ; DVðxÞ ¼ 0 0 Q3 Q4 where F and V are the m · m matrices defined by, h i h i oF i ðxÞ and V ¼ oV F ¼ oxji ðxÞ with 1 6 i; j 6 m: oxj Further, F is non-negative, V is a non-singular M-matrix and Q3 ; Q4 are matrices associated with the transition terms of the model, and all eigenvalues of Q4 have positive real parts. Finally, the following stability result follows. Theorem 8 (van den Driessche and Watmough [28]). Consider an arbitrary disease transmission model with fi(x) satisfying the axioms (A1)–(A5). If x is a DFE of the model, then x is LAS if R0 ¼ qðFV 1 Þ < 1, but unstable if R0 > 1. In the above, q represents the spectral radius.
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
497
Appendix B. Proof of Theorem 1 using comparison theorem Theorem B.1. The SFE of the basic models (1)–(4), given by E0 , is GAS in D if Rs < 1. Proof. The above result will now be proved using a comparison theorem. As in Section 2.2.1, the equations for the rate of change of the smoking components (S, Qt) of the model (1)–(4), are re-written in terms of ! dSðtÞ SðtÞ b 0 SðtÞ P dt ¼ ðF V Þ 1 ; dQt ðtÞ N Qt ðtÞ 0 0 Qt ðtÞ dt
where the matrices F and V are as defined in Section 2.2.1. Since P 6 N (for all t P 0) in D, it follows that ! dSðtÞ SðtÞ dt 6 ðF V Þ : ðB:1Þ dQt ðtÞ Qt ðtÞ dt
Using the fact that the eigenvalues of the matrix F V all have negative real parts (see local stability result in Lemma 1, where q(FV1) < 1 if Rs < 1, which is equivalent to F V having eigenvalues with negative real parts when Rs < 1 [28]), it follows that the linearized differential inequality system (B.1) is stable whenever Rs < 1. Consequently, ðSðtÞ; Qt ðtÞÞ ! ð0; 0Þ as t ! 1. It follows by comparison theorem (see, for instance, [23, p. 31] and [26], Theorem B.1; Appendix B) that ðSðtÞ; Qt ðtÞÞ ! ð0; 0Þ. Substituting S = Qt = 0 in the first and the last equations of model (1)–(4), gives P(t) ! P*, and Qp(t) ! 0 as t ! 1. Thus, ðP ðtÞ; SðtÞ; Qt ðtÞ; Qp ðtÞÞ ! ðP ; 0; 0; 0Þ as t ! 1 for Rs < 1, so that E0 is GAS if Rs < 1. h Appendix C. Positivity of D1, D2 and D3 Given, D1 ¼ h½g2 J 5 þ a1 J 4 þ wðJ 1 J 5 a1 J 3 Þ; D2 ¼ ð1 hÞ½J 2 J 5 a2 J 4 þ wða2 J 3 þ g1 J 5 Þ; it follows that J 1 J 5 a1 J 3 ¼ ðg1 þ l þ c1 þ sÞðl þ a1 þ a1 Þ a1 c1 ð1 r1 Þ ¼ l2 þ ðg1 þ a1 þ s þ a2 þ c1 Þl þ g1 a1 þ g1 a2 þ c1 a2 þ sa1 þ sa2 þ c1 a1 r1 > 0 and J 2 J 5 a2 J 4 ¼ ðg2 þ l þ c2 þ snÞðl þ a1 þ a1 Þ a2 c2 ð1 r2 Þ ¼ l2 þ ðg2 þ a1 þ sn þ a2 þ c2 Þl þ g2 a1 þ g2 a2 þ c2 a1 þ sna1 þ sna2 þ c2 a2 r2 > 0: Thus, D1 > 0 and D2 > 0. Furthermore, D 3 ¼ J 1 J 2 J 5 g2 a2 J 3 a 1 J 2 J 3 g 1 g2 J 5 g1 a1 J 4 a2 J 1 J 4 ¼ ðg1 þ l þ c1 þ sÞðg2 þ l þ c2 þ snÞðl þ a1 þ a1 Þ g2 a2 c1 ð1 r1 Þ a1 c1 ðg2 þ l þ c2 þ snÞð1 r1 Þ g1 g2 ðl þ a1 þ a2 Þ g1 a1 c2 ð1 r2 Þ a2 c2 ðg1 þ l þ c1 þ sÞð1 r2 Þ ¼ l3 þ ðc1 þ a2 þ c2 þ a1 þ sn þ s þ g1 þ g2 Þl2 þ ½g2 a2 þ c1 c2 þ c1 g2 þ c1 a2 þ g2 a1 þ c2 a2 r2 þ c2 a1 þ g1 c2 þ g1 a1 þ g1 a2 þ c1 a1 r1 þ ðc2 þ a2 þ g2 þ a1 Þs þ ðs2 þ fa2 þ c1 þ g1 þ a1 gsÞnl þ ðfa2 þ a1 gs2 þ fg1 a2 þ c1 a2 þ c1 a1 r1 þ g1 a1 gsÞn þ a2 c2 c1 r2 þ ðg2 a1 þ c2 a2 r2 þ g2 a2 þ c2 a1 Þs þ c1 g2 a2 r1 þ a1 c1 g2 r1 þ a1 c1 c2 r1 þ a2 c2 g1 r2 þ g1 c2 a1 r2 > 0 Thus, D3 is also positive.
498
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
Appendix D. Proof of Theorem 4 using comparison theorem Theorem 4. The SFE of the extended model (13) is GAS if Rsd < 1. Proof. Here, too, the equations for the rate of change of the smoking components (S 1 ; S 2 ; Qt ; C) of the model 13 can be written in terms of 0 dS 1 ðtÞ 1 0 1 1 0 10 S 1 ðtÞ ð1 hÞb ð1 hÞwb 0 0 S 1 ðtÞ dt B dS 2 ðtÞ C B C C CB B C B B S 2 ðtÞ C B hb hwb 0 0 CB S 2 ðtÞ C B dt C C 1 P B C CB B dQ ðtÞ C ¼ ðF V ÞB B Q ðtÞ C B CB Q ðtÞ C; B t C N 0 0 0 0 @ A @ A @ A t t @ dt A dCðtÞ CðtÞ CðtÞ 0 0 0 0 dt
where the matrices F and V are as defined in Section 3.2.1. Since P 6 N (for all t P 0) in M, it follows that 0 dS ðtÞ 1 1 0 1 S 1 ðtÞ B dSdtðtÞ C B S ðtÞ C B 2 C B dt C B 2 C ðD:1Þ C: B dQ ðtÞ C 6 ðF V ÞB @ Qt ðtÞ A B t C @ dt A dCðtÞ CðtÞ dt
Using the fact that the eigenvalues of the matrix F V all have negative real parts (see local stability result in Lemma 2, where q(FV1) < 1 if Rsd < 1, which is equivalent to F V having eigenvalues with negative real parts when Rsd < 1 [28]), it follows that the linearized differential inequality system (D.1) is stable whenever Rsd < 1. Consequently, ðS 1 ðtÞ; S 2 ðtÞ; Qt ðtÞ; CðtÞÞ ! ð0; 0; 0; 0Þ as t ! 1. It follows, by comparison theorem, that ðS 1 ðtÞ; S 2 ðtÞ; Qt ðtÞ; CðtÞÞ ! ð0; 0; 0; 0Þ. Substituting S1 = S2 = Qt = C = 0 in the first and penultimate equations of the model (13) gives P(t) ! P* and Qp(t) ! 0 as t ! 1. Thus, ðP ðtÞ; S 1 ðtÞ; S 2 ðtÞ; Qt ðtÞ; Qp ðtÞ; CðtÞÞ ! ðP ; 0; 0; 0; 0; 0Þ as t ! 1 for Rsd < 1, so that X0 is GAS if Rsd < 1. h References [1] R.M. Anderson, R.M. May (Eds.), Population Biology of Infectious Diseases, Springer-Verlag, Berlin, Heidelberg, New York, 1982. [2] R.M. Anderson, R.M. May, Infectious Diseases of Humans, Oxford University Press, London, 1991. [3] D. Bernoulli, Essai d’une nouvelle analyse de la mortalite´ cause´e par la petite ve´role, et des avantages de l’inoculation pour la pre´venir. Histoire de l’Acade´mie Royale des Sciences. me´moire de mathe´matiques et de physique, Paris (1760). pp. 1–45. [4] F. Brauer, C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology, Text in Applied Mathematics, Springer, 2000. [5] J. Carr, Applications Centre Manifold Theory, Springer-Verlag, New York, 1981. [6] C. Castillo-Chavez, K. Cooke, W. Huang, S.A. Levin, Results on the dynamics for models for the sexual transmission of the human immunodeficiency virus, Appl. Math. Lett. 2 (1989) 327–331. [7] C. Castillo-Chavez, B. Song, Dynamical models of tuberculosis and their applications, Math. Biosci. Eng. 1 (2) (2004) 361–404. [8] C. Castillo-Chavez, S. Blower, P. van den Driessche, D. Kirschner, A.A. Yakubu (Eds.), Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Introduction, The IMA Volumes in Mathematics and its Applications, vol. 125, Springer, New York, 2002. [9] C. Castillo-Chavez, H. Wenzhang, L. Jia, Competitive exclusion and coexistence of multiple strains in an SIS STD model, SIAM J. Appl. Math. 5 (1999) 1790–1811. [10] G.C. Castillo, S.G. Jordan, and A.H. Rodriguez. Mathematical models for the dynamics of tobacco use, recovery and relapse. Technical Report Series, BU-1505-M. Department of Biometrics, Cornell University, 2000. [11] O. Diekmann, J.A.P. Heesterbeek, J.A.J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, J. Math. Biol. 28 (1990) 503–522. [12] E.H. Elbasha, A.B. Gumel, Theoretical assessment of public health impact of imperfect Prophylactic HIV-1 vaccines with therapeutic benefits, Bull. Math. Biol. 68 (2005) 577–614. [13] A.B. Gumel, C.C. McCluskey, J. Watmough, An SVEIR model for assessing potential impact of an imperfect anti-SARS vaccine, Math. Biosci. Eng. 3 (3) (2006) 485–512. [14] A.B. Gumel, McCluskey, C. Connell, P. van den Driessche, Mathematical study of a staged-progression HIV model with imperfect vaccine. Bull. Math. Biol., 2006, 10.1007/s11538-006-9095-7.
O. Sharomi, A.B. Gumel / Applied Mathematics and Computation 195 (2008) 475–499
499
[15] A.B. Gumel, X.W. Zhang, P.N. Shivakumar, M.L. Garba, B.M. Sahai, A new mathematical model for assessing therapeutic strategies for HIV infection, J. Theor. Med. 4 (2) (2002) 147–155. [16] J.K. Hale, Ordinary Differential Equations, John Wiley and Sons, New York, 1969. [17] Health Canada . New directions for tobacco control in Canada: a strategy. Steering Committee of the National Strategy to reduce Tobacco Use in Canada, Ottawa, Ontario. Available from: http://dsp-psd.pwgsc.gc.ca/Collection/H39-505-1999E.pdf, 2005. [18] H.M. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599–963. [19] W.O. Kermack, A.G. McKendrick, A contribution to the mathematical theory of epidemics, Proc. Royal Soc. London. 115 (1927) 700–721. [20] W.O. Kermack, A.G. McKendrick, Contributions to the mathematical theory of epidemics, part II, Proc. Royal Soc. London 138 (1932) 55–83. [21] W.O. Kermack, A.G. McKendrick, Contributions to the mathematical theory of epidemics, part III, Proc. Royal Soc. London 141 (1933) 94–112. [22] M. Kgosimore, E.M. Lungu, The effects of vaccination and treatment on the spread of HIV/AIDS, J Biol Syst. 12 (4) (2004) 399–417. [23] V. Lakshmikantham, S. Leela, A.A. Martynyuk, Stability Analysis of Nonlinear Systems, Marcel Dekker Inc., New York, 1989. [24] A. Perelson, P. Nelson, Mathematical analysis of HIV-1 dynamics in vivo, SIAM Rev. 41 (1999) 3–44. [25] T.L. Scott, S. Robert, V.A. Kirsh, Beliefs about tobacco industry (mal)practices and youth smoking behaviour: insight for future tobacco control campaigns (Canada), Cancer Causes Control 17 (2006) 705–711. [26] R.J. Smith, S.M. Blower, Could disease-modifying HIV vaccine cause population-level pervasity? Lancet Infect. Dis. 4 (2004) 636– 639. [27] US Department of Health and Human Services. The health consequences of smoking: a report of the surgeon general. US Department of Health and Human Services, Public Health Service, Centers for Disease Control and Prevention, National Center for Chronic Disease Prevention and Health Promotion, Office of Smoking and Health, Atlanta, GA. Available from: http:// www.cdc.gov/Tobacco/sgr/sgr_2004/consumerpiece/SGR2004_Whatitmeanstoyou.pdf, 2004. [28] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci. 180 (2002) 29–48. [29] J .Walton, History of smoking. http://www.forestonline.org/output/page34.asp, 2003. [30] World Health Organization. Global smoking statistics. http://quitsmoking.about.com/cs/antismoking/a/statistics.htm, 2002. [31] World Health Organization. World cancer report. http://www.who.int/mediacentre/news/releases/2003/pr27/en, 2003.
