Optimal intervention strategy for prevention tuberculosis using a smoking-tuberculosis model

Author's Accepted Manuscript

Optimal Intervention Strategy for Prevention Tuberculosis using a Smoking-Tuberculosis Model Sunhwa Choi, Eunok Jung, Seok-Min Lee

www.elsevier.com/locate/yjtbi

PII: DOI: Reference:

S0022-5193(15)00257-X http://dx.doi.org/10.1016/j.jtbi.2015.05.022 YJTBI8204

To appear in:

Journal of Theoretical Biology

Received date: 11 November 2014 Revised date: 13 May 2015 Accepted date: 14 May 2015 Cite this article as: Sunhwa Choi, Eunok Jung, Seok-Min Lee, Optimal Intervention Strategy for Prevention Tuberculosis using a SmokingTuberculosis Model, Journal of Theoretical Biology, http://dx.doi.org/10.1016/j. jtbi.2015.05.022 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Optimal Intervention Strategy for Prevention Tuberculosis using a Smoking-Tuberculosis Model Sunhwa Choia , Eunok Junga , Seok-Min Leeb,∗ a Department

of Mathematics, Konkuk University, 120 Neungdong-ro, Gwngjin-gu, Seoul, 143-701, Republic of Korea b Department of Liberal Arts, College of Engineering, Hongik University, 94 Wausan-ro, Mapo-gu, Seoul, 121-791, Republic of Korea

Abstract In this paper, we developed a dynamic model of smoking-tuberculosis (TB) transmission in South Korea, and investigated the effects of control strategies on the number of incidence of TB using optimal control theory. Model parameters regarding TB and smoking are estimated through least-squares fitting to real data. We considered three TB controls (distancing, case-finding, and case-holding) and two smoking controls (distancing and quitting), in order to minimize the number of exposed and infectious individuals and the cost of control. Numerical simulations for the various control strategies highlight that implementing the smoking controls, not with TB controls, can effectively reduce the incidence of TB transmission. Keywords: Epidemic model, Optimal control theory, Parameter Estimation 1. Introduction Tuberculosis (TB) is a major global health threat, with 8.6 million new cases and 1.3 million deaths reported in 2012, along with one third of the total population carrying a latent infection [52]. Susceptible individuals, in airborne contact with active pulmonary TB patients, are infected with Mycobacterium tuberculosis, usually through the lungs. Tobacco smoking is also a significant cause of disease, disability, and death. Considering both the direct smoking and the effects of second-hand smoke, smoking kills nearly six million people every year. The World Health Organization (WHO) has warned that of the one billion smokers worldwide, up to half of those will die of tobacco-related diseases in the long run [54]. TB is a serious problem in South Korea, which has the highest TB incidence, prevalence, and mortality rates among Organization for Economic Cooperation and Development (OECD) nations [51]. The smoking rate in Korea is in the intermediate level of the OECD member countries, although it is very high among male adults ∗ Corresponding

author Email address: [email protected] (Seok-Min Lee )

Preprint submitted to Journal of Theoretical Biology

May 20, 2015

[29, 30, 31, 32, 33, 34, 35]. Thus, TB and tobacco smoking are very important national health issues in Korea. It has long been reported that there is a relation between tobacco use and pulmonary TB, and there is increasing evidence of this association. Many researchers conducted meta-analyses or cohort studies using large-size samples, which increase the credence of the results. Smoking has been positively associated with the development of TB infection and incident active TB [2, 7, 16, 18, 19, 21, 24, 26, 41, 45], TB relapse [5, 45], and TB-related mortality rates [15, 16, 18, 21, 22, 24, 36]. Bates et al. [7] used meta-analytic methods to quantify the relationship between smoking and TB infection, disease, and mortality, and concluded that smoking is a risk factor for TB infection and TB disease. The research of Jee et al. [18] on smoking and TB among South Koreans was by a massive cohort (1,294,504 people) over a long period (1992-2006). The study reported an association of smoking with incident TB, and TB recurrence in men, as well as a 60% increase in the risk of death. In [55], Yen et al. tracked down over 5,000 individuals who successfully completed TB treatment and were confirmed as cured through bacteriological testing. Of these, 1.5% were observed to have developed a recurrent case of TB, with regular tobacco smokers twice as likely to suffer a recurrence compared with individuals who never smoked or had quit. Indeed, some papers have reported that quitting smoking can reduce TB incidence (Lin et al. [23]) and TB mortality (Wen et al. [48]). Based on the above, one may conclude that efforts to reduce smoking will support the TB intervention program in controlling the incidence, relapse, and mortality rate of TB. Some researchers have suggested the cessation of smoking as a TB control program: for example, Basu et al. [6] pointed out that smoking effects delay the millennium development target for TB mortality, while others [37, 38] proposed that active support for smoking cessation should be included in TB treatment programs. Sereno [39], Slama et al. [42, 43, 44] and Enarson et al. [12] suggested detailed methods for incorporating the cessation of smoking into TB care. Our first objective is to develop a mathematical model of smoking TB transmission based on South Korean data. To set up this dynamic TB-smoking model, we first need to consider a mathematical model for smoking. A social habit such as smoking is usually modeled as an infectious disease in the literature, as smoking is transmitted via humanto-human interactions, particularly with peer pressure. Little work has been done on the mathematical modeling of tobacco smoking compared with the modeling of epidemic diseases. Since Castillo-Garsow et al. [10] proposed simple models describing the acquisition and cessation of smoking, the mathematical models developed in later studies [3, 4, 17, 40, 56, 57] have been based on SIR or SEIR model, while some use finer compartmental decompositions. Stability analyses and optimal control theory have been applied to smoking models using numerical simulations, although the connection with real data is not always clear. Van Voorn and Kooi [46] presented an eco-epidemiological dynamic model that combined smoking with resource considerations, using parameters based on statistics from the Netherlands. A number of authors have presented mathematical models that describe the effect of smoking on the transmission of infectious diseases or other social behaviors. Acevedo2

Estefania et al. [1] developed a model for smoking and lung cancer, while Bhunu and Mushayabasa [8] combined smoking and alcoholism. Bhunu et al. [9] reported a combined smoking and tuberculosis mathematical model, and conducted a stability analysis on the equilibria, determination of the reproductive number, and some numerical simulations to evaluate the effects of smoking on active TB cases. The mathematical modeling of TB transmission in South Korea was studied with real data by Whang et al. [49] and Choi and Jung [11]. According to their research, the TB model consists of four compartments, the susceptible, the exposed, the infectious, and the permanently latent (SEIL). This is a modified SEIR model, with main parameters data fitted to South Korean TB data and optimal treatment strategies proposed. Three control mechanisms (distancing, case finding, and case holding) were considered to minimize the number of latent and infectious individuals. In this study, we extend Choi and Jung’s TB model [11] to a smoking-TB model: a SEIL model for tuberculosis, with each compartment of TB transmission sub-divided into smoking and non-smoking (potentially smoking) groups. Model parameters were estimated or fitted to South Korean data. This is the first study to develop a dynamic Korea-specific smoking-TB model. It has been reported that smoking enhances TB transmission in several ways as stated previously, and the model reflects this assumption. We also apply optimal control theory to the developed mathematical model, and thus propose optimal intervention strategies. Five controls are considered: TB distancing control, TB case-finding control, TB case-holding control, smoking distancing control, and smoking cessation control. The purpose of these controls is to minimize the number of infectious and exposed individuals and the cost of implementing the control treatment. In particular, we focus on the effects of controls on smoking behavior to the TB transmission. Numerical results are given for four different optimal control strategies, three of which adopt only smoking-related controls. The distinctive feature of our study is that the smoking control strategies can reduce TB transmission. This paper is organized as follows. In Section 2, we introduce a smoking-TB model that considers the effects of smoking on the transmission of TB, and then describe the parameter estimation. A probabilistic sensitivity analysis for the parameters is also reported. In Section 3, we present the model along with its optimal controls and numerical studies for four strategies. Our conclusions and final remarks are given in Section 4. Mathematical details for optimal control theory can be found in the appendix. 2. Materials and Methods 2.1. Mathematical Model We extend Choi and Jung’s TB model [11] to develop a smoking-TB model for South Korea. The smoking effect is modeled by smoking transmission and quitting rates, as developed in Bhunu’s model [9]. In this section, the TB model is presented first, and then the coupled smoking-TB model will be described. For the TB model, individuals are classified into compartments as susceptible (S), exposed (E), active TB infectious (I), and low-risk latent (L), based on the pathological 3

bN



μ

S β T NI

 E k

μ

/

O c pr



α

/

σ

μ

I

(1− p )r

d

/ /

#  L μ

 Figure 1: Flowchart of the TB transmission dynamics describing model structure of TB progression of state classes of individuals.

characteristics of TB. The total population size is given by N = S + E + I + L. We denote by b and μ the natural birth rate and death rate respectively. It is assumed that individuals are mixed homogeneously. The susceptible individuals progress to the class E at rate β T NI when exposed to infectious individuals, where β T is the average infected susceptible number of TB transmission and NI is the probability to contact with an infectious individual in the whole population. Individuals in E progress to the active TB infectious class (I) at the rate k, and to the latent class (L) at the rate α. A proportion r of infectious individuals move back to the class E at the incomplete treatment rate p, while the rest, at treatment success rate (1 − p), progress to the latent class L. Infectious individuals die from TB at the rate d. The relapse rate σ represents the proportion of individuals transmitting from the class L to the class E. Figure 1 illustrates the dynamics of TB transmission. The smoking-TB model divides the population into two groups: potential smokers (non-smokers) NP and smokers NS . Each group is subdivided into the four classes of the SEIL TB model described above. Therefore, a total of eight compartments are considered: susceptible potential smokers SP , TB exposed potential smokers who are not infectious but have a high risk of acquiring TB infections EP , TB infectious potential smokers IP , low-risk permanently latent potential smokers L P , smoking susceptible individuals SS , TB exposed smokers ES , TB infectious smokers IS , and permanent latent smokers LS . Thus, the total population is given by N, where N = NP + NS = 4

μ

o

bN



o

μ

μ

o

α

o

SP

( I +η I ) βT P N S

o

βS

O

o b

IP

o

(1− p )r

"  LP μ

α

o

φβ T

 /

ES

q σ

pr



μ

SS

q



d

/

N β S NS

EP k

NS N

ωk

( I P + η IS ) N

μ

O

IS

(1− p)rθ

N β S NS

/ " 

q



/

b prθ



ξq

/

σ

μ εd

/ /

LS



μ

Figure 2: Smoking-TB transmission diagram. The model is given by the system of equations (1).

(SP + EP + IP + L P ) + (SS + ES + IS + LS ). We also denote S, E, I, and L, the TB transmission groups combining smokers and potential smokers, with S = SP + SS , E = EP + ES , I = IP + IS , and L = L P + LS . The flow of chains SP -EP -IP -L P and SS -ES -IS -LS follows the SEIL TB model with a slight amendment on parameters considering the affection from smoking. Transmission between compartments in different chains are supposed parallel. We assume that potential smokers acquire a smoking habit from their contact with smokers, in a similar manner to the transmission of infectious disease. All the newcomers into the system are supposed to be in SP . The transmission of smoking is considered class-wise: SP to SS , EP to ES , and L P to LS . We assume that a non-smoking active TB infectious individual tends not to acquire a smoking habit. Potential smokers who are exposed to smokers or to scenes of smoking in media, which are assumed proportional to the number of smokers, acquire a smoking habit at rate β S NNS , where β S is the average number of effective contacts. Smokers in classes SS , ES and LS quit smoking at a rate q and return to the potential smoker classes. Because infectious smokers with TB symptoms are assumed to quit more than smokers in other classes, individuals in IS quit smoking at rate ξq with ξ > 1. We also assume that smoking increases the infectiousness of TB, and η > 1 denotes this increase in TB transmission rate for TB infectious smokers. Hence, β T SNP ( IP + η IS ) is the number of newly infected nonsmoking individuals per unit time. Susceptible smok5

ers have more chance of being infected, TB exposed smokers tend to develop active TB faster, and smokers with active TB have a lower recovery rate. It has been reported that smoking increases the disease mortality rate. These effects of smoking on TB transmission are reflected in the transmission rates for smokers, as multipliers to the rates. The positive numbers φ > 1, ω > 1, θ < 1, and ε > 1 are additional factors applied to the rates β T SNS ( IP + η IS ), k, r, and d, respectively. The structure of the proposed model is shown in Figure 2, and is described by the following system of equations: ⎧ SP SP ⎪  ⎪ S NS − μSP + qSS , = bN − β ( I + η I ) − β ⎪ T P S S P ⎪ N N ⎪ ⎪ ⎪ ⎪ ⎨ E = β SP ( I + η I ) − β EP N − (α + k + μ) E + qE + prI + σL , T P P P P S S S S P N N ⎪ ⎪ IP = kEP − (μ + r + d) IP + ξqIS , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ LP = (1 − p)rIP + αEP − (σ + μ) L P − β S L P NS + qLS , N ⎧ S S P ⎪ ⎪ SS = β S NS − φβ T S ( IP + η IS ) − (q + μ)SS , ⎪ (1) ⎪ N N ⎪ ⎪ ⎪ ⎪ ⎨ E = φβ SS ( I + η I ) + β EP N − (α + ωk + μ + q) E + prθ I + σL , T P S S S S S S S N N ⎪ ⎪ IS = ωkES − (ξq + μ + εd + rθ ) IS , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ LS = (1 − p)rθ IS + αES + β S L P NS − (μ + q + σ ) LS , N where NP = SP + EP + IP + L P , NS = SS + ES + IS + LS , and N = NP + NS . The model parameters are summarized in Table 2 in subsection 2.2.2. 2.2. Parameter Estimation The smoking-TB model (1) is based on the TB model in [11], and inherits the parameters b, μ, β T (as β in [11]), r, and p. The smoking acquisition and cessation parameters β S and q are obtained from a pure smoking model that is fitted to real smoking data using the least-squares method. Some smoking-TB parameters are obtained from previous studies [9, 7, 18]. Once the other parameters have been determined, the parameters related to TB incidence α, k, and σ are fitted to TB incidence and TB relapse data using a least-squares fitting method. The Matlab routine lsqcurvefit is used for this data fitting process. 2.2.1. Smoking Parameters To estimate the smoking transmission rate β S and cessation rate q, we use data from the National Health and Nutrition Examination Survey, conducted by the Korea Centers for Disease Control and Prevention (KCDC) (Table 1) [20]. The data include the proportion of smokers (present smokers) and former smokers (people who quit smoking) in certain 6

sample populations. The survey was conducted eight times during the period 2001-2012. As the corresponding survey data describe information on adult smoking habits, the data are adjusted to include the infant and adolescent population, and are normalized to the 2005 population structure. Table 1: National Health and Nutrition Examination Survey, Korea Centers for Disease Control and Prevention (KCDC)

Year

2001

2005

2007

2008

2009

2010

2011

2012

ratio of smokers (%)

23.1

22.0

19.3

21.2

20.8

21.0

20.6

19.7

ratio of former smokers (%)

6.6

11.1

13.4

12.8

12.6

12.7

12.8

12.0

We develop a smoking-only model to estimate β S and q from the data. We include the number of smokers NS , number of former smokers NQ , and number of people who have never smoked NN . Then, the number of potential smokers is NP = NQ + NN , and the total population is N = NS + NQ + NN . The smoking model can be formulated as follows:  NN = bN −

βS NN NS − μNN , N

βS ( NN + NQ ) NS − (μ + q) NS , N β  = qNS − S NQ NS − μNQ . NQ N NS =

(2)

Starting from the population data of 2001, we can analyze the above system of differential equations numerically to determine values of β S and q that ensure the ratio of N smokers ( NNS ) and the ratio of former smokers ( NQ ) best fit the given data. Figure 3 depicts the fitted values (−) and the real data (∗). The best fitted values for the transmission rate of smoking and the quitting rate are β S = 0.0720 and q = 0.0495. 2.2.2. Determining the other parameters The smoking-TB model is a modification of the TB-only model in [11], from which it inherits certain TB parameters. We use the same values for the birth rate b, the natural death rate μ, the TB transmission coefficient β T (β in [11]), the treatment success rate 1 − p, the treatment rate for the infectious r, and TB-induced death rate d. The parameters α, k, and σ are given by fitting to the real data, because these values are not explicitly determined from direct observations and may be changed by our modifications. The TB incidence and relapse data of South Korea from the WHO [53] are used to determine α, k, and σ. The incidence and relapse numbers correspond to k( EP + ωES ) and σ( L P + LS ), respectively, in the smoking-TB model of (1). Figure 4 compares the

7

0.2 0.3

fitted data

0.18

0.28

0.16

0.26

0.14

0.24

0.12

NQ/N

0.22

NS/N

fitted data

0.2

0.1

0.18

0.08

0.16

0.06

0.14

0.04

0.12

0.02

0.1 2001

2003

2005

2007

2009

0 2001

2011 2012

2003

2005

Time (year)

2007

2009

2011 2012

Time (year)

Figure 3: Smoking rate (left) and quitting rate (right) data and fitted curves. The best fit is given by β S = 0.0720 and q = 0.0495.

incidence data with values of k( EP + ωES ) computed by the fitted variables, and compares the relapse data with σ( L P + LS ) over the simulation period [2001, 2012]. The fitted values are k = 0.0588, α = 0.4818, σ = 0.0001626. The parameters in our smokingTB model are summarized in Table 2. Figure 5 displays the fractions of the population of classes combined smoking-nonsmoking, S(t)/N (t), E(t)/N (t), I (t)/N (t), and L(t)/N (t). The fraction of susceptible individuals is around 0.6 over the simulated time. The fraction of exposed and infectious individuals increases, while the fraction of permanently latent individuals remains constant at around 0.4 for the simulation period. Note that the prevalence rate, ( E(t) + I (t) + L(t))/N (t), is around 0.4, which is in good agreement with the data in [52]. From this result, we can predict that the number of infectious individuals in South Korea is gradually increasing. 4

x 10

6000

4 5000

number of relapse

incidence of TB

3.5

3

4000

3000

2000

2.5 1000 fitted data 2 2001

2004

2007 Time (year)

2010

fitted data 0 2001

2012

2004

2007 Time (year)

2010

2012

Figure 4: The reported TB incidence (left) and relapse (right) data and their fitted curves (kE(t) and σL(t), respectively). The fitted values are k = 0.0588, α = 0.4818 and σ = 0.0001626.

8

1

0.015

0.01

0.6

E/N

S/N

0.8

0.4

0.005

0.2 0 2001

2004

2007

0 2001

2010 2012

2010 2012

1 0.8 L/N

1.5 I/N

2007 L/N

−3

x 10

1 0.5 0 2001

2004

0.6 0.4 0.2

2004 2007 Time (year)

0 2001

2010 2012

2004 2007 Time (year)

2010 2012

Figure 5: State variables with fitted parameters for the smoking-TB model from 2001 to 2012.

2.3. Bootstrap Bootstrapping is a general approach to statistical inference based on building a sample distribution for a statistic by resampling from the given data [14]. To ensure the reliability of the estimates of β S , q, k, α, and σ obtained through the least-squares fitting, we employ a parametric bootstrap technique. The proportion of smokers (present smokers) and former smokers (people who have stopped smoking) are resampled from normal distribution. Since the smoking rate and the ratio of former smokers are from samples with appropriate weight values, the distributions of the proportions follow modified binomial distributions, which can be approximated by the normal distributions due to the big sample size. The mean and standard deviation are calculated from the original data. Using the resampled data, we estimate β S and q. The new incidence numbers of active TB individuals and relapses are generated from Poisson distributions with means that give the best fit to the observed data. Values for k, α, and σ are then estimated using the generated data. The distribution of estimates obtained from 1000 generated datasets is shown in Figure 6 and the mean, standard deviation (S.D.), and 95% confidence interval (CI) are listed in Table 3. The results demonstrate the confidence in the parameter estimates.

9

Table 2: Model parameters and their interpretations. 

Symbol b

Population

μ

natural death rate

Smoking ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

frequency

treatment success rate treatment rate for latent treatment rate for infectious

σ



Smoking-TB

progression rate

r

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

Unit

References

0.0210

yr −1

[11], [49]

0.0126

yr −1

[11], [49]

10

transmission coefficient of TB

α

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

Value

birth rate

⎧ βT ⎪ ⎪ k ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ 1− p TB

Definition

relapse rate

[11], [49]

0.0588

yr −1

data fitting

0.8293

yr −1

[11], [49]

0.4818

yr −1

data fitting

0.9874

yr −1

[11], [49]

0.000 162 6

yr −1

data fitting

yr −1

[49]

d

TB induced death rate

0

βS

transmission coefficient of smoking

0.0720

q

quitting rate

0.0495

η

transmission increase factor

1.04

[9]

ϕ

transmission increase factor

1.73

[7]

ω

progression increase factor

2

[6]

θ

treatment decrease factor

0.75

[9]

ξ

increased smoking quitting factor

1.05

[9]

ε

increased TB induced factor

1.6

[18]

90

90

90

90

90

80

80

80

80

80

70

70

70

70

70

60

60

60

60

60

50

50

50

50

50

40

40

40

40

40

30

30

30

30

30

20

20

20

20

20

10

10

10

10

10

0

0

0.05

β

s

0.1

data fitting yr −1

0.15

0

0

0.05

q

0.1

0

0.05

0.06

k

0.07

0

0.4

0.45

0.5

α

0.55

0

data fitting

1.6 1.62 1.64 1.66

σ

−4

x 10

Figure 6: Bootstrap distributions with 1000 resampling simulations according to a normal distribution (left two) and Poisson distribution (right three).

10

Table 3: Parameter estimates and mean, standard deviation, and 95% confidence interval calculated from the simulation study.

Parameter

Estimate

Mean

S.D.

95% CI

βS

0.0720

0.0719

0.0051

(0.0619, 0.0819)

q

0.0495

0.0495

0.0032

(0.0432, 0.0558)

k

0.0588

0.0587

0.0024

(0.0540, 0.0634)

α

0.4818

0.4814

0.0212

(0.4401, 0.5224)

σ

0.0001626

0.0001626

0.0000009

(0.0001609, 0.0001643)

2.4. Sensitivity Analysis Sensitivity analyses quantify the level of uncertainty in any type of complex model. The objective is to identify critical model inputs (parameters and initial conditions) and quantify how input uncertainty impacts the model outcomes [25]. It is often difficult to accurately estimate the parameters involved in system (1). Because there is no data for parameters η, θ, and ξ, we examine how sensitive the state variables SP , EP , IP , L P , SS , ES , IS , and LS are to these parameters. The partial rank correlation coefficient (PRCC) measures the statistical influence of each of the three input parameters and model variables. Figure 7 shows the PRCC of all state variables with respect to (η, θ, ξ ) at the final simulation time, t = 2012 yr. Table 4 shows the PRCC results given by Latin hypercube sampling (LHS). As shown in Figure 7 and Table 4, the factors denoting increased transmission η and decreased treatment θ are highly correlated with the model variables. The number of infectious individuals (I = IP + IS ) is positively correlated with η, and negatively correlated with θ. Results of the sensitivity analysis suggest that all outcome measure, the state variables (SP , EP , IP , L P , SS , ES , IS , and LS ) are sensitive to changes in η and θ. That is, the two parameters, η and θ are important contributors to uncertainty in the model. 3. Optimal Treatment Strategies for Smoking-TB in South Korea In this section, we formulate an optimal control problem in order to minimize the number of individuals who are actively TB infectious (I = IP + IS ) and TB latent (E = EP + ES ), taking into account the cost of treatment and social efforts within a finite time period. We calibrate our model using the control functions ui (t) (i = 1, 2, 3, 4, 5), which quantify the intervention strategies or policies. To investigate the optimal treatment policy, our smoking-TB model is extended from 2014 to 2030. We assume the same parameter values determined from the period [2001, 2012] for the years after 2012. Before we investigate optimal control strategies, it is essential to examine the behavior of the number of infected individuals E + I in the smoking-TB model when no controls 11

Figure 7: Sensitivity analysis: general LHS scheme and PRCC performed on the model. The y-axis indicates PRCC values of the set (SP , EP , IP , L P , SS , ES , IS , LS ) for model parameters η, θ, and ξ at the final simulation time (t = 2012 yr ). Table 4: Model parameters used in sensitivity analysis and PRCC values of SP , EP , IP , L P , SS , ES , IS , and LS at t = 2012 for three perturbed parameters η, θ, ξ. The range (minimum/maximum) of the three perturbed parameters and their baseline value are shown on the right. Sample size n = 1,000

Parameter

η

θ

ξ

SP

EP

IP

LP

SS

ES

IS

LS

−0.96922∗

0.96669∗

0.96894∗

0.97071∗

−0.96925∗

0.96219∗

0.95849∗

0.97275∗

0.91645∗

−0.92370∗

−0.93139∗

−0.91427∗

0.91629∗

−0.91515∗

−0.97008∗

−0.90685∗

0.06639

−0.06601

0.06129

−0.04343

0.05812

−0.06501

−0.19833

−0.10692

Minimum

Baseline

Maximum

0.001

1.04

2

0.001

0.75

1

0.001

1.05

2

∗Significant (p-value < 0.01)

are imposed. The estimated number of infected in E + I satisfying system (1) for the long-term period (up to the year 2100) is depicted in Figure 8. The number of infected 12

individuals increases considerably without additional TB and/or smoking intervention efforts. This demonstrates that intervention strategies and control plans are definitely required in South Korea. 5

10

x 10

9 8 7

E+I

6 5 4 3 2 1 0 2014 2020

2030

2040

2050 2060 Time (year)

2070

2080

2090

2100

Figure 8: The prediction of the state variable I + E up to 2100. We assumed same parameter values determined for [2001, 2012].

3.1. Smoking-TB model with controls We denote by u1 (t), u2 (t), u3 (t), u4 (t), and u5 (t) the TB distancing control, TB casefinding control, TB case-holding control, smoking distancing control, and quit smoking control, respectively. The TB distancing control u1 (t) represents efforts intended to reduce the number of susceptible individuals who become infected. Examples include the isolation of infectious people and educational programs or campaigns for healthy individuals. Because TB is transmitted aerially, such campaigns may promote the wearing of face masks and discourage careless coughing or sneezing. The control factor (1 − u1 (t)) reduces the TB transmission rate β T . When u1 = 0, no external control on TB distancing is enforced. Values of u1 approaching 1 signify that full TB distancing efforts are being enforced. The TB case finding control u2 (t) describes efforts made to identify latent individuals at high risk of developing TB, and the treatment of such individuals, who may benefit from preventive intervention. Directly observed therapy (DOTS) has long been recommended by the WHO as a case-finding control for TB. The regulation of periodic checkups, massive X-ray programs at school, easy access to clinical investigations or sputum collection points, and even door-to-door investigations in high-risk areas have been presented as active TB case-finding policies ([27, 28, 50]). The rate modification factor is (1 + u2 (t)) and hence the full effort is quantified as twice the rate of transmission from E = EP + ES to L = L P + LS . 13

The TB case holding control u3 (t) represents efforts to prevent the failure of treatment in TB infectious individuals in order to reduce the reinfection rate. Because TB treatment takes rather a long period and a high cost, control efforts aim to ensure that patients complete their treatment, such as by periodically confirming their drug intake and providing financial support to help pay for the treatment. The rate modifying factor is of the form (1 − u3 (t)) multiplied to p; thus, if the case-holding control u3 (t) is near 1, the treatment failure rate is very low. The smoking distancing control u4 (t) represents the effort of reducing the number of potential smokers that start smoking. One efficient method is to prevent exposure to scenes that glamorize smoking in the mass media, such as in books, columns, magazines, dramas, shows, movies or advertisements, as these can cause susceptible individuals to take up smoking [13]. Another approach may be to emphasize the risks of smoking by highlighting the development of various respiratory and cardiovascular diseases such as lung cancer and heart attacks, as well as cancer in various organs of the human body. The control factor (1 − u4 (t)) is multiplied by the smoking transmission rate β S . The smoking cessation control u5 (t) represents the effort to encourage smokers to quit through governmental policies and enforcements, and via educational programs and/or campaigns aimed at smokers. For instance, non-smoking areas can be extended, with severe fines imposed for violations, smoking cessation clinics could be expanded, cigarette taxes can be raised, and health warnings can be displayed on cigarette packs [47, 54]. We use the control factor (1 + gu5 (t)), where g measures the effectiveness of u5 . When u5 (t) is close to 1, the rate at which smokers quit smoking increases to around (1 + g). We assume that g = 1, and hence the full effort has the potential to double the number of smoking individuals who quit. The system of differential equations for the smoking-TB model with five controls can

14

be written as: ⎧ SP SP ⎪ ⎪ SP = bN − (1 − u1 ) β T ( IP + η IS ) − (1 − u4 ) β S NS − μSP + (1 + u5 )qSS , ⎪ ⎪ N N ⎪ ⎪ ⎪ ⎪ S E P P ⎪  ⎪ ⎪ EP = (1 − u1 ) β T ( IP + η IS ) − (1 − u4 ) β S NS − ((1 + u2 )α + k + μ) EP ⎪ N N ⎪ ⎨ + (1 + u5 )qES + (1 − u3 ) prIP + σL P , ⎪  ⎪ I = kEP − (μ + r + d) IP + ξ (1 + u5 )qIS , ⎪ P ⎪ ⎪ ⎪ ⎪ LP ⎪  ⎪ NS = ( 1 − ( 1 − u ) p ) rI + ( 1 + u ) αE − ( σ + μ ) L − ( 1 − u ) β L ⎪ 3 2 P P P 4 S P ⎪ ⎪ N ⎪ ⎩ + (1 + u5 )qLS , ⎧ SP SS ⎪  ⎪ S N = ( 1 − u ) β − φ ( 1 − u ) β ( IP + η IS ) − ((1 + u5 )q + μ)SS , ⎪ T 4 1 S S S ⎪ N N ⎪ ⎪ ⎪ ⎪ SS EP ⎪  ⎪ E NS − ((1 + u2 )α + ωk + μ) ES = φ ( 1 − u ) β ( I + η I ) + ( 1 − u ) β ⎪ T P 1 4 S S S ⎪ N N ⎪ ⎨ − (1 + u5 )qES + (1 − u3 ) prθ IS + σLS , ⎪  ⎪ IS = ωkES − (ξ (1 + u5 )q + μ + εd + θr ) IS , ⎪ ⎪ ⎪ ⎪ ⎪ LP ⎪ ⎪ LS = (1 − (1 − u3 ) p)rθ IS + (1 + u2 )αES + (1 − u4 ) β S NS ⎪ ⎪ ⎪ N ⎪ ⎩ − ( μ + (1 + u5 ) q + σ ) L S , where NP = SP + EP + IP + L P ,

NS = SS + ES + IS + LS ,

(3)

N = NP + NS .

Optimal intervention strategies are obtained by minimizing the active TB infectious and TB-latent individuals, while keeping the intervention cost low under the proposed scenarios. We minimize the sum of infectious and latent individuals (E + I = EP + ES + IP + IS ) and the linear combination of quadratic terms u2i (t) (i = 1, 2, 3, 4, 5), as the cost of controls within the given time period. Hence, the objective functional to be minimized is given by J ( u1 , u2 , u3 , u4 , u5 )  t  f 1 1 1 1 1 2 2 2 2 2 EP + IP + ES + IS + B1 u1 + B2 u2 + B3 u3 + B4 u4 + B5 u5 dt, = 2 2 2 2 2 t0

(4)

where t0 = 2014 and t f = 2030 are the initial and the final simulation times, respectively. The constants Bi (i = 1, 2, 3, 4, 5) are the relative weight coefficients that balance the cost and importance of intervention strategies. The optimal control problem aims to find a vector of optimal control functions (u1∗ (t), u2∗ (t), u3∗ (t), u4∗ (t), u5∗ (t)) such that J (u1∗ , u2∗ , u3∗ , u4∗ , u5∗ ) = min J (u1 , u2 , u3 , u4 , u5 ), Ω

L1 ( t

where Ω = {(u1 , u2 , u3 , u4 , u5 ) ∈ 0 , t f ) | 0 ≤ ui ≤ 1, i = 1, 2, 3, 4, 5}. Pontryagin’s Maximum Principle is used to solve this optimal control problem, giving the necessary conditions for optimal control. The characteristics of optimal control, which guide the numerical computation to obtain optimal control, are provided in Appendix A. 15

3.2. Optimal intervention strategies We investigate optimal treatment intervention in the following four cases: ⎧Strategy I: TB distancing, TB case finding, TB case holing, smoking distancing, ⎪ ⎪ ⎪ ⎪ and quit-smoking controls (u1 (t), · · · , u5 (t)) ⎨ Strategy II: smoking distancing and quit smoking controls (u4 (t) and u5 (t)) ⎪ ⎪ ⎪ ⎪ ⎩Strategy III: smoking distancing control (u4 (t))

Strategy IV: quit-smoking control (u5 (t)) We assume that the baseline of the balancing coefficients Bi in (4) associated with controls ui (t) takes the order of 104 for i = 1, 2, 3, 4, 5 throughout the strategies. In Strategy I we try other sets of coefficients as well. Figure 9 and Figure 11 illustrate optimal controls and the corresponding state variables (combined as sums of smokers and potential smokers) as functions of time under Strategy I and Strategy II, respectively. In the state variable graphs (the four frames at the bottom of Figures 9 and 11) , the solid black curves represent the number of individuals in S, E, I, and L without controls. These can be compared with the dashed blue state curves with optimal controls. In Figure 10, optimal control problem for Strategy I is implemented with various sets of weight constants (scenarios). Four sets of B j = ( B1 , B2 , B3 , B4 , B5 ), j = 1, 2, 3, 4 are chosen as shown in Table 5. The figure shows optimal controls and the corresponding state values (E + I) for each set of weights (blue dashed curves for B1 , red dotted for B2 , green dash dotted for B3 , magenta dotted with large dots for B4 ), as well as the number of individuals in class E + I in the no-control case (solid black curve). In Strategy I with scenario B1 (baseline: blue dashed curve in Figure 9 and 10), the full TB distancing and TB case-finding controls (u1 (t) and u2 (t)) are applied until around 2021 and 2018, respectively. The curves of these optimal controls then decrease to the lower bound. In contrast, the TB case holding, smoking distance, and smoking cessation controls (u3 (t), u4 (t), and u5 (t)) are never applied to their full extent. This suggests that TB isolation and case-finding efforts play a significant role (particularly TB isolation) in reducing the number of infected and infectious individuals. The TB case-holding control is considered to be a less important factor. Controls on smoking habits (u4 (t) and u5 (t)) Table 5: Description of Strategies (use of controls) and Scenarios (values of weight constants)

Controls used Scenarios

Strategy I

Strategy II

Strategy III

Strategy IV

u1 , u2 , u3 , u4 , u5

u4 , u5

u4

u5

B1

B1

B1

B1 =

(104 , 104 , 104 , 104 , 104 )

B2 =

(5 · 104 , 104 , 104 , 104 , 104 )

B3 = (104 , 104 , 104 , 103 , 103 ) B4 = (104 , 104 , 104 , 102 , 102 )

16

remain near the lower bound throughout the time period, indicating that these controls play a minor role in minimizing the number of TB infected and infectious individuals. In the scenarios B2 , B3 , B4 , the weight constants for the TB controls are bigger than those for the smoking controls, because we are trying to determine cases in which smoking controls are used more. Only scenario B2 , where the weight constant of u1 (t) is especially high, sees the use of u2 (t) outstrip that of u1 (t). Smoking controls u4 (t) and u5 (t) are only fully used in scenario B4 , where their weight constants make them extremely cheap. In addition, the cost of TB case-holding control u3 (t) remains steady in the scenarios considered here. The state value of E + I is stable throughout the scenarios. From the various scenarios in Strategy I, the TB isolation control was found to be most efficient if the weight constants are around even, and moreover, the smoking controls start working if we give very low weight constants to them. Therefore, it is only necessary to explore the smoking-control related strategies to investigate whether, and to what extent, u4 (t) and u5 (t) influence TB transmission (Strategy II, Strategy III, and Strategy IV). We now consider Strategy II, which Figure 11 depicts. The full effort of both controls is applied until 2024, and then the controls decrease to the lower bound. With the application of smoking-control related strategies and optimum controls, it can be seen that E(t) and I (t) gradually but meaningfully decrease, although this is less effective than the additional use of controls targeting TB (Strategy I). In Figure 12, optimal controls u4 (t) (Strategy III) and u5 (t) (Strategy IV) are plotted in the first two graphs. For the corresponding states, we refer to the third graph in Figure 12, which illustrates the ratio of infected individuals ( E + I )/N for Strategy I (blue dashed curve), Strategy II (red dotted curve), Strategy III (green dash-dotted curve), Strategy IV (magenta dotted curve with bigger dots), and the case of no controls (black solid curve). The graph enables us to compare the effectiveness of these five cases. Note that, when only one smoking control is considered (i.e., Strategy III or IV), the ratio ( E + I )/N decreases over time, and the number of infected and infectious individuals also decreases. More specifically, the effect of optimal controls is most effective in Strategy I, when all five TB and smoking controls are used. The ratio of infected people in the population is just 1.13% of those in the no-control case in the year 2030. In Strategy II (using u4 and u5 ), the proportion of infected people in the population is 42.55% of that in the no-control case. In Strategies III and IV, the percentage of the population ratio that are infected are 55.72% and 62.77% of that with no control. The results show that interventions involving only smoking could reduce the ratio of infected individuals. For a vector of weight constants (a scenario) B = ( B1 , B2 , B3 , B4 , B5 ), the relative 1 2 total cost (TC) is defined as TC = ∑5i=1 ∑m j=1 2 Bi ui ( t j ) dt where the given time period [2014, 2030] is divided into m subintervals [t j − 1, t j ], j = 1, 2, . . . , m, with the time step 1 2 dt = 0.01. Also, the relative cost of each control ui (t) is defined as ∑m j=1 2 Bi ui ( t j ) dt. The graph at the bottom of Figure 12 displays the relative cost of controls in each strategy (baseline) in the form of a stacked bar graph, while the values of relative TCs for all scenarios are given in the Table 6. At the graph, we can see that applying only one smoking 17

control (Strategy III or IV) costs less than using all the controls, but using two smoking controls (Strategy II) costs nearly twice as much as Strategy III or Strategy IV, and even more than the case Strategy I, while the effect is only slightly better than Strategy III or Strategy IV (see the third graph in Figure 12). In view of the cost, Strategy III and Strategy IV show more efficiency than Strategy II. To summarize, the numerical results for various strategies and scenarios indicate that the TB isolation strategy is the most effective, since the control u1 is used the most, unless the weight constant for u1 is excessively higher than others. The ratio ( E + I )/N decreases to almost zero in Strategy I. It is remarkable that the number of infected and infectious individuals decreases in the strategies that use only smoking-related controls. These results should be considered by governments when developing TB control policies. Table 6: The relative TC over [2014, 2030] are calculated for the four strategies, and four scenarios in Strategy I.

Strategy I

TC

Strategy II

Strategy III

Strategy IV

B1

B2

B3

B4

B1

B1

B1

8.672 · 104

1.9886 · 105

8.58 · 104

8.262 · 104

1.1258 · 105

6.3508 · 104

6.4287 · 104

4. Conclusions It has been reported that smoking increases the risk of TB. Few studies have set up a mathematical model of TB associated with smoking, and no previous research has investigated optimal control strategies based on mathematical modeling. This paper has described a smoking-TB dynamic model and presented optimal control strategies on smoking and TB with regard to South Korea. A comprehensive and continuous model for the transmission of TB has been proposed that reflects the influence of smoking. The model includes eight compartments of individuals according to their TB transmission stage and smoking status. Parameter values were determined using WHO-provided Korean TB incidence and relapse data from 2001 to 2011, as well as KCDC smoking habit survey results from 2001, 2005, and each year from 2007 to 2012. The main parameters were derived from the real data, and then the confidence level of the parameters was reinforced by bootstrapping. Parameters that could not be obtained from this data were inferred from previous research. We investigated optimal treatment strategies for smoking and TB in South Korea. Five control functions were introduced to minimize the number of exposed and infected individuals and actively infectious individuals. Three controls were related to the transmission of TB, and the other two involved smoking. Four optimal control strategies were investigated, and four scenarios with various combinations of weight constants were explored. The numerical simulations on Strategies 18

I-IV represented by Figure 9-12 shows that TB isolation control is the most efficient, but smoking isolation and cessation controls can also reduce the numbers of TB infected and infectious individuals. The importance of this work is twofold. First, a trailblazing Korea-specific smokingTB model has been developed. Second, the results highlight that optimal strategies to decrease smoking can reduce the prevalence of infected and infectious individuals. The serious nature of TB in South Korea prompted the Korean government to develop a national-level TB eradication plan that is continually adjusted for efficiency. Currently this plan does not focus on smoking-related efforts. Thus the authors hope that this study will be helpful in providing a scientific basis for the determination of future governmental TB policies. Appendix A. Characteristics of Optimal Control The optimal controls are obtained from the necessary conditions, which come from the Pontryagin’s maximum Principle. The conditions are interpreted to minimize the Hamiltonian H with respect to u1 , u2 , u3 , u4 and u5 , where the Hamiltonian is defined as follows: H = EP (t) + IP (t) + ES (t) + IS (t) 1 1 1 1 1 + B1 u21 (t) + B2 u22 (t) + B3 u23 (t) + B4 u24 (t) + B5 u25 (t) 2 2 2 2 2 + λ1 (t)SP (t) + λ2 (t) EP (t) + λ3 (t) IP (t) + λ4 (t) LP (t)

(A.1)

+ λ5 (t)SS (t) + λ6 (t) ES (t) + λ7 IS + λ8 (t) LS (t).

Our objective functional J is defined as J ( u1 , u2 , u3 , u4 , u5 ) =

 t f t0

H (t)dt,

in the domain Ω = {(u1 , u2 , u3 , u4 , u5 ) ∈ L1 (t0 , t f ) | 0 ≤ ui ≤ 1, i = 1, 2, 3, 4, 5}, where t0 is the initial time and t f is the final time. Theorem 1. There exist optimal controls u1∗ , u2∗ , u3∗ , u4∗ , u5∗ , and state variables S∗P , EP∗ , IP∗ , L∗P , SS∗ , ES∗ , IS∗ , and L∗S that minimize the objective functional J (u2 , u2 , u3 , u4 , u5 ) over Ω. Given

19

these optimal controls, there exist adjoint variables λ1 , λ2 , . . . , λ8 satisfying  dλ1 ( IP + η IS ) NS = λ1 (1 − u1 ) β T + (1 − u4 ) β S +μ dt N N ( I + η IS ) N − λ2 (1 − u1 ) β T P − λ5 (1 − u4 ) β S S , N N  dλ2 N = − 1 + λ2 (1 − u4 ) β S S + ((1 + u2 )α + k + μ) dt N N − λ3 k − λ4 (1 + u2 ) α − λ6 (1 − u4 ) β S S , N  dλ3 SP SP = − 1 + λ1 (1 − u1 ) β T − λ2 (1 − u1 ) β T + (1 − u3 ) pr + λ3 (μ + r + d) dt N N S S − λ4 (1 − (1 − u3 ) p )r + λ5 φ (1 − u1 ) β T S − λ6 φ (1 − u1 ) β T S , N N  dλ4 NS NS = − λ2 σ + λ4 ( σ + μ ) + (1 − u4 ) β S − λ8 (1 − u4 ) β S , dt N N  dλ5 S E L =λ1 (1 − u4 ) β S P − (1 + gu5 )q + λ2 (1 − u4 ) β S P + λ4 (1 − u4 ) β S P (A.2) dt N N N  S ( I + η IS ) − λ5 (1 − u4 ) β S P − φ (1 − u1 ) β T P − ((1 + gu5 )q + μ) N N  ( IP + η IS ) EP L − λ6 φ (1 − u1 ) β T + (1 − u4 ) β S − λ8 (1 − u4 ) β S P , N N N  dλ6 S E = − 1 + λ1 (1 − u4 ) β S P + λ2 (1 − u4 ) β S P − (1 + gu5 )q dt N N L S + λ4 (1 − u4 ) β S P − λ5 (1 − u4 ) β S P N N  EP − λ6 (1 − u4 ) β S − (ωk + (1 + u2 )α + μ) − (1 + gu5 )q − λ7 ωk N  LP , − λ8 (1 + u2 ) α + (1 − u4 ) β S N

20

 dλ7 SP SP = − 1 + λ1 (1 − u1 ) β T η + (1 − u4 ) β S dt N N  S E − λ2 (1 − u1 ) β T P η − (1 − u4 ) β S P − λ3 ξ (1 + gu5 )q N N  LP SP SS + λ4 (1 − u4 ) β S − λ5 (1 − u4 ) β S − φ (1 − u1 ) β T η N N N  S E − λ6 φ(1 − u1 ) β T S η + (1 − u4 ) β S P + (1 − u3 ) prθ N N + λ7 (ξ (1 + gu5 )q + μ + εd + θr )  LP − λ8 (1 − (1 − u3 ) p)rθ + (1 − u4 ) β S , N  SP EP LP dλ8 = λ1 (1 − u4 ) β S + λ2 (1 − u4 ) β S + λ4 (1 − u4 ) β S − (1 + gu5 )q , dt N N N  SP EP − λ5 (1 − u4 ) β S − λ6 (1 − u4 ) β S +σ N N  LP − λ8 (1 − u4 ) β S − (μ + (1 + gu5 )q + σ) N with the transversality conditions λi (t f ) = 0, i = 1, . . . , 8.

(A.3)

The optimal controls can be obtained from the following equations: u1∗ u2∗ u3∗ u4∗ u5∗



1 SP SS = (λ2 − λ1 ) β T ( IP + η IS ) + (λ6 − λ5 )φβ T ( IP + η IS ) B1 N N 1 = {(λ2 − λ4 )αEP + (λ6 − λ8 )αES } B2 1 = {(λ2 − λ4 ) IP + (λ6 − λ8 )θ IS } pr B3 1 β = {(λ5 − λ1 )SP + (λ6 − λ2 ) EP + (λ8 − λ4 ) L P } S NS B4 N gq = {(λ5 − λ1 )SS + (λ6 − λ2 ) ES + (λ7 − λ3 )ξ IS + (λ8 − λ4 ) LS } B5

21

(A.4)

Proof. By differentiating the Hamiltonian, we obtain the adjoint system: dλ1 dt dλ2 dt dλ3 dt dλ4 dt dλ5 dt dλ6 dt dλ7 dt dλ8 dt

∂H = − ∂S , λ1 (t f ) = 0, P ∂H = − ∂E , λ2 (t f ) = 0, P ∂H = − ∂I , λ3 (t f ) = 0, P ∂H = − ∂L , λ3 (t f ) = 0, P ∂H = − ∂S , λ5 (t f ) = 0, S ∂H , λ6 (t f ) = 0, = − ∂E S ∂H = − ∂I , λ7 (t f ) = 0, S ∂H = − ∂L , λ8 (t f ) = 0. S

Therefore, the adjoint variables λ1 (t), · · · , λ8 (t) satisfy (A.2). By differentiating the Hamiltonian with respect to the controls, we have the following optimality conditions at the optimal controls u1∗ , · · · , u5∗ : ∂H S S = B1 u1 + λ1 β T P ( IP + η IS ) − λ2 β T P ( IP + η IS ) ∂u1 N N S S + λ5 φβ T S ( IP + η IS ) − λ6 φβ T S ( IP + η IS ) N N ∂H = B2 u2 − λ2 αEP + λ4 αEP − λ6 αES + λ8 αES ∂u2 ∂H = B3 u3 − λ2 prIP + λ4 prIP − λ6 prθ IS + λ8 prθ IS ∂u3 ∂H S E L = B4 u4 + λ1 β S P NS + λ2 β S P NS + λ4 β S P NS ∂u4 N N N S E L − λ5 β S P NS − λ6 β S P NS − λ8 β S P NS N N N ∂H = B5 u5 + λ1 gqSS + λ2 gqES + λ3 ξgqIS + λ4 gqLS ∂u5 − λ5 gqSS − λ6 gqES − λ7 ξgqIS − λ8 gqLS

22

=0 =0 =0

=0

=0

(A.5)

From (A.5), we derive the optimal control condition: 

1 SP SS ∗ u1 = (λ2 − λ1 ) β T ( IP + η IS ) + (λ6 − λ5 )φβ T ( IP + η IS ) B1 N N 1 u2∗ = {(λ2 − λ4 )αEP + (λ6 − λ8 )αES } B2 1 u3∗ = {(λ2 − λ4 ) IP + (λ6 − λ8 )θ IS } pr B3 1 β u4∗ = {(λ5 − λ1 )SP + (λ6 − λ2 ) EP + (λ8 − λ4 ) L P } S NS B4 N gq ∗ u5 = {(λ5 − λ1 )SS + (λ6 − λ2 ) ES + (λ7 − λ3 )ξ IS + (λ8 − λ4 ) LS } B5

(A.6)

We modify the equation (A.4) for the optimal controls considering the bounds. Let a and b be the lower and upper bounds, respectively. In this manuscript, the values a = 0 and b = 1 are adopted. The optimal controls are characterized as follows:    S 1 ∗ u1 = min b, max a, (λ2 − λ1 ) β T P ( IP + η IS ) B1 N

SS +(λ6 − λ5 )φβ T ( IP + η IS ) N   1 ∗ u2 = min b, max a, {(λ2 − λ4 )αEP + (λ6 − λ8 )αES } B2   1 ∗ u3 = min b, max a, (A.7) {(λ2 − λ4 ) IP + (λ6 − λ8 )θ IS } pr B3   βS 1 ∗ u4 = min b, max a, {(λ5 − λ1 )SP + (λ6 − λ2 ) EP + (λ8 − λ4 ) L P } NS B4 N   gq u5∗ = min b, max a, {(λ5 − λ1 )SS + (λ6 − λ2 ) ES + (λ7 − λ3 )ξ IS B5 +(λ8 − λ4 ) LS } Acknowledgements The research work of Jung was supported by the Korea National Research Foundation (NRF) grant funded by the Korea government (MEST) (No. 2012R1A2A2A01011725). The research work of Lee was supported by the Hongik University new faculty research support fund.

23

References [1] Acevedo-Estefania, C. A., Gonzalez, C., Rios-Soto, K. R., Summerville, E. D., Song, B., and Castillo-Chavez, C. A mathematical model for lung cancer: The effects of second-hand smoke and education. Biometrics Unit Technical Report BU-1525-M, Cornell University, 2000. [2] Alcaide, J., Altet, M. N., Plans, P., Parron, I., Folguera, L., Salto, E., Dominguez, A., Pardell, H., and Salleras, L. Cigarette smoking as a risk factor for tuberculosis in young adults: a casecontrol study. Tubercle and Lung Disease, 77(2):112–116, 1996. [3] Alkhudhari, Z., Al-Sheikh, S., and Al-Tuwairqi, S. Global dynamics of a mathematical model on smoking. ISRN Applied Mathematics, 2014, 2014. [4] Alkhudhari, Z., Al-Sheikh, S., and Al-Tuwairqi, S. Stability analysis of a giving up smoking model. International Journal of Applied Mathematical Research, 3(2):168–177, 2014. [5] d’Arc Lyra Batista, J., de Alencar Ximenes, R. A., Rodrigues, L. C., et al. Smoking increases the risk of relapse after successful tuberculosis treatment. International journal of epidemiology, 37(4):841–851, 2008. [6] Basu, S., Stuckler, D., Bitton, A., and Glantz, S. A. Projected effects of tobacco smoking on worldwide tuberculosis control: mathematical modelling analysis. BMJ, 343, 2011. [7] Bates, M. N., Khalakdina, A., Pai, M., Chang, L., Lessa, F., and Smith, K. R. Risk of tuberculosis from exposure to tobacco smoke: a systematic review and metaanalysis. Archives of internal medicine, 167(4):335–342, 2007. [8] Bhunu, C. P. and Mushayabasa, S. A theoretical analysis of smoking and alcoholism. Journal of Mathematical Modelling and Algorithms, 11(4):387–408, 2012. [9] Bhunu, C. P., Mushayabasa, S., and Tchuenche, J. M. A theoretical assessment of the effects of smoking on the transmission dynamics of tuberculosis. Bulletin of mathematical biology, 73(6):1333–1357, 2011. [10] Castillo-Garsow, C., Jordán-Salivia, G., and Rodriguez-Herrera, A. Mathematical models for the dynamics of tobacco use, recovery, and relapse. Public Health, 84(4): 543–547, 1997. [11] Choi, S. and Jung, E. Optimal tuberculosis prevention and control strategy from a mathematical model based on real data. Bulletin of mathematical biology, pages 1–24, 2014.

24

[12] Enarson, D. A., Slama, K., and Chiang, C.-Y. Providing and monitoring quality service for smoking cessation in tuberculosis care [Educational Series: tobacco and tuberculosis. Serialised guide. Tobacco cessation interventions for tuberculosis patients. Number 6 in the series]. The International Journal of Tuberculosis and Lung Disease, 11(8):838–847, 2007. [13] Fielding, J. E. Smoking: health effects and control. (Second of two parts). New England Journal of Medicine, 313(9):555–61, 1985. [14] Fox, J. Bootstrapping regression models. An R and S-PLUS Companion to Applied Regression: A Web Appendix to the Book. Sage, Thousand Oaks, CA., 2002. URL http://cran.r-project.org/doc/contrib/Fox-Companion/ appendix-bootstrapping.pdf . [15] Gajalakshmi, V., Peto, R., Kanaka, T. S., and Jha, P. Smoking and mortality from tuberculosis and other diseases in India: retrospective study of 43 000 adult male deaths and 35 000 controls. The Lancet, 362(9383):507–515, 2003. [16] Hassmiller, K. M. The association between smoking and tuberculosis. salud pública de méxico, 48:s201–s216, 2006. [17] Huo, H.-F. and Zhu, C.-C. Influence of relapse in a giving up smoking model. In Abstract and Applied Analysis, volume 2013. Hindawi Publishing Corporation, 2013. [18] Jee, S. H., Golub, J. E., Jo, J., Park, I. S., Ohrr, H., and Samet, J. M. Smoking and risk of tuberculosis incidence, mortality, and recurrence in South Korean men and women. American journal of epidemiology, 170(12):1478–1485, 2009. [19] Kolappan, C. and Gopi, P. G. Tobacco smoking and pulmonary tuberculosis. Thorax, 57(11):964–966, 2002. [20] Korea Center for Disease Control and Prevention. Korea national health and nutrition examination survey (KNHANES II, KNHANES III, KNHANES IV-1, KNHANES IV-2, KNHANES IV-3, KNHANES V-1, KNHANES V-2, KNHANES V-3), 2001, 2005, 2007, 2008, 2009, 2010, 2011, 2012. URL https://knhanes.cdc.go.kr/ . [21] Leung, C. C., Li, T., Lam, T. H., Yew, W. W., Law, W. S., Tam, C. M., Chan, W. M., Chan, C. K., Ho, K. S., and Chang, K. C. Smoking and tuberculosis among the elderly in Hong Kong. American journal of respiratory and critical care medicine, 170 (9):1027–1033, 2004. [22] Lin, H.-H., Ezzati, M., and Murray, M. Tobacco smoke, indoor air pollution and tuberculosis: a systematic review and meta-analysis. PLoS medicine, 4(1):e20, 2007. [23] Lin, H.-H., Murray, M., Cohen, T., Colijn, C., and Ezzati, M. Effects of smoking and solid-fuel use on COPD, lung cancer, and tuberculosis in China: a time-based, multiple risk factor, modelling study. The Lancet, 372(9648):1473–1483, 2008. 25

[24] Lin, H.-H., Ezzati, M., Chang, H.-Y., and Murray, M. Association between tobacco smoking and active tuberculosis in Taiwan: prospective cohort study. American journal of respiratory and critical care medicine, 180(5):475–480, 2009. [25] Marino, S., Hogue, I. B., Ray, C. J., and Kirschner, D. E. A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of theoretical biology, 254(1):178–196, 2008. [26] Maurya, V., Vijayan, V., Shah, A., et al. Smoking and tuberculosis: an association overlooked. The International journal of tuberculosis and Lung Disease, 6(11):942–951, 2002. [27] Miller, A. C., Golub, J. E., Cavalcante, S. C., Durovni, B., Moulton, L. H., Fonseca, Z., Arduini, D., Chaisson, R. E., and Soares, E. C. C. Controlled trial of active tuberculosis case finding in a Brazilian favela. The International Journal of Tuberculosis and Lung Disease, 14(6):720–726, 2010. [28] Murray, C. J. L. and Salomon, J. A. Expanding the WHO tuberculosis control strategy: rethinking the role of active case-finding [The Pittsfield Lecture]. The International Journal of Tuberculosis and Lung Disease, 2(9s1):S9–S15, 1998. [29] Organization for Economic Cooperation and Development (OECD) Staff. Health at a Glance 2001. 2001. doi: http://dx.doi.org/10.1787/health_glance-2001-en. [30] Organization for Economic Cooperation and Development (OECD) Staff. Health at a Glance 2003: OECD Indicators. 2003. doi: http://dx.doi.org/10.1787/health_ glance-2003-en. [31] Organization for Economic Cooperation and Development (OECD) Staff. Health at a Glance 2005: OECD Indicators. 2005. doi: http://dx.doi.org/10.1787/ 9789264012639-en. [32] Organization for Economic Cooperation and Development (OECD) Staff. Health at a Glance 2007: OECD Indicators. 2007. doi: http://dx.doi.org/10.1787/health_ glance-2007-en. [33] Organization for Economic Cooperation and Development (OECD) Staff. Health at a Glance 2009: OECD Indicators. 2009. doi: http://dx.doi.org/10.1787/health_ glance-2009-en. [34] Organization for Economic Cooperation and Development (OECD) Staff. Health at a Glance 2011: OECD Indicators. 2011. URL http://www.oecd-ilibrary. org/social-issues-migration-health/health-at-a-glance-2011_health_ glance-2011-en.

26

[35] Organization for Economic Cooperation and Development (OECD) Staff. Health at a Glance 2013: OECD Indicators. 2013. doi: http://dx.doi.org/10.1787/health_ glance-2013-en. [36] Reed, G. W., Choi, H., Lee, S. Y., Lee, M., Kim, Y., Park, H., Lee, J., Zhan, X., Kang, H., Hwang, S., et al. Impact of diabetes and smoking on mortality in tuberculosis. PloS one, 8(2):e58044, 2013. [37] Saad, T. and Tirkey, A. S. Association between pulmonary tuberculosis and smoking: a case control study. Indian Journal of Community Health, 25(4):340–347, 2013. [38] Schneider, N. K. and Novotny, T. E. Addressing smoking cessation in tuberculosis control. Bulletin of the World Health Organization, 85(10):820–821, 2007. [39] Sereno, A. B., Soares, E. C. C., Lapa e Silva, J. R., Nápoles, A. M., Bialous, S. A., da Costa e Silva, V. L., and Novotny, T. E. Feasibility study of a smoking cessation intervention in Directly Observed Therapy Short-Course tuberculosis treatment clinics in Rio de Janeiro, Brazil. Revista Panamericana de Salud Pública, 32(6):451–456, 2012. [40] Sharomi, O. and Gumel, A. B. Curtailing smoking dynamics: a mathematical modeling approach. Applied Mathematics and Computation, 195(2):475–499, 2008. [41] Singh, P. N., Yel, D., Kheam, T., Hurd, G., and Job, J. S. Cigarette smoking and tuberculosis in Cambodia: findings from a national sample. Tob Induc Dis, 11(8), 2013. [42] Slama, K., Chiang, C.-Y., and Enarson, D. A. Helping patients to stop smoking [Educational Series: tobacco and tuberculosis. Serialised guide. Tobacco cessation interventions for tuberculosis patients. Number 5 in the series]. The International Journal of Tuberculosis and Lung Disease, 11(7):733–738, 2007. [43] Slama, K., Chiang, C.-Y., and Enarson, D. A. Introducing brief advice in tuberculosis services [Educational Series: tobacco and tuberculosis. Serialised guide. Tobacco cessation interventions for tuberculosis patients. Number 3 in the series]. The International Journal of Tuberculosis and Lung Disease, 11(5):496–499, 2007. [44] Slama, K., Chiang, C.-Y., and Enarson, D. A. Tobacco cessation and brief advice [Educational Series: tobacco and tuberculosis. Serialised guide. Tobacco cessation interventions for tuberculosis patients. Number 4 in the series]. The International Journal of Tuberculosis and Lung Disease, 11(6):612–616, 2007. [45] Slama, K., Chiang, C.-Y., Enarson, D. A., Hassmiller, K., Fanning, A., Gupta, P., and Ray, C. Tobacco and tuberculosis: a qualitative systematic review and metaanalysis [Review Article]. The International Journal of Tuberculosis and Lung Disease, 11(10):1049–1061, 2007. 27

[46] van Voorn, G. A. K. and Kooi, B. W. Smoking epidemic eradication in a ecoepidemiological dynamical model. Ecological Complexity, 14:180–189, 2013. [47] Wakefield, M., Morley, C., Horan, J., and Cummings, K. The cigarette pack as image: new evidence from tobacco industry documents. Tobacco Control, 11(suppl 1):i73– i80, 2002. [48] Wen, C.-P., Chan, T.-C., Chan, H.-T., Tsai, M.-K., Cheng, T.-Y., and Tsai, S.-P. The reduction of tuberculosis risks by smoking cessation. BMC infectious diseases, 10(1): 156, 2010. [49] Whang, S., Choi, S., and Jung, E. A dynamic model for tuberculosis transmission and optimal treatment strategies in South Korea. Journal of theoretical biology, 279(1): 120–131, 2011. [50] World Health Organization and others. Guidelines for intensified tuberculosis casefinding and isoniazid preventive therapy for people living with HIV in resourceconstrained settings. 2011. [51] World Health Organization and others. Tuberculosis prevalence surveys: a handbook. Geneva: World Health Organization, 2011. [52] World Health Organization and others. Tuberculosis. Fact sheet N 104. , 2014. URL http://www.who.int/mediacentre/factsheets/fs104/en/ . [53] World Health Organization and others. Global health observatory data repository, Tuberculosis, 2014. URL http://apps.who.int/gho/data/node.main.1315?lang= en. [54] World Health Organization and others. Tobacco. Fact sheet N 339. , 2014. URL http://www.who.int/mediacentre/factsheets/fs339/en/ . [55] Yen, Y.-F., Yen, M.-Y., Lin, Y.-S., Lin, Y.-P., Shih, H.-C., Li, L.-H., Chou, P., and Deng, C.-Y. Smoking increases risk of recurrence after successful anti-tuberculosis treatment: a population-based study. The International Journal of Tuberculosis and Lung Disease, 18(4):492–498, 2014. [56] Zaman, G. Optimal campaign in the smoking dynamics. Computational and mathematical methods in medicine, 2011, 2011. [57] Zaman, G. Qualitative behavior of giving up smoking model. Malaysian Mathematical Sciences Society, 34(2):403–415, 2011.

28

Bulletin of the

Strategy I u1(t)

1 0.5

u2(t)

0 2014 1

u3(t)

2020

2022

2024

2026

2028

2030

2016

2018

2020

2022

2024

2026

2028

2030

2016

2018

2020

2022

2024

2026

2028

2030

2016

2018

2020

2022

2024

2026

2028

2030

2016

2018

2020

2022

2024

2026

2028

2030

0.5 0 2014 0.4

u4(t)

2018

0.5 0 2014 1

0.2 0 2014 0.4

u5(t)

2016

0.2 0 2014

Time (year)

7

5

x 10

7

4

x 10

6 5

3

E

S

4 2

3 2

1

without control with control

0 2014

2018

2022

2026

1 0 2014

2030

4

2022

2026

2030

2022

2026

2030

7

x 10

x 10 2

5 4

1.5

3

L

I

2018

1

2 0.5

1 0 2014

2018

2022

2026

2030

Time (year)

0 2014

2018

Time (year)

Figure 9: Optimal controls with the weight coefficients B1 = (104 , 104 , 104 , 104 , 104 ) and state variables for Strategy I (with five optimal controls) from 2014 to 2030. The above five graphs show the optimal control variables depending on time. In the below four graphs, the state variables (blue dotted lines) in Strategy I are compared to those without controls (black solid lines).

29

u1(t)

1 0.5

u2(t)

0 2014 1

u3(t)

2020

2022

2024

2026

2028

2030

2016

2018

2020

2022

2024

2026

2028

2030

2016

2018

2020

2022

2024

2026

2028

2030

2016

2018

2020

2022

2024

2026

2028

2030

2016

2018

2020

2022

2024

2026

2028

2030

0.5 0 2014 1

u4(t)

2018

0.5 0 2014 1

0.5 0 2014 1

u5(t)

2016

0.5 0 2014

Time (year)

5

8

x 10

7

6

without controls B1 =(B1=104, B2=104, B3=104, B4=104, B5=104) B2 =(B1=5⋅104, B2=104, B3=104, B4=104, B5=104)

5

E+I

B3 =(B1=104, B2=104, B3=104, B4=103, B5=103) B4 =(B1=104, B2=104, B3=104, B4=102, B5=102)

4

3

2

1

0 2014

2018

2022

2026

2030

Time (year)

Figure 10: Optimal controls and the corresponding E + I for Strategy I as functions of time are displayed in the top frame and the bottom frame, respectively. The blue dashed, red dotted, green dashdotted and magenta dotted (with large dots) curves represent the cases with the weight coefficients B1 = (104 , 104 , 104 , 104 , 104 ), B2 = (5·104 , 104 , 104 , 104 , 104 ), B3 = (104 , 104 , 104 , 103 , 103 ) and B4 = (104 , 104 , 104 , 102 , 102 ), respectively. The black solid curve represents the case without controls.

30

Strategy II 1

u4(t)

0.8 0.6 0.4 0.2 0 2014

2016

2018

2020

2016

2018

2020

2022

2024

2026

2028

2030

2022

2024

2026

2028

2030

1

u5(t)

0.8 0.6 0.4 0.2 0 2014

Time (year)

7

5

x 10

7

4

x 10

6 5

3

E

S

4 2

3 2

without control with control

1 0 2014

2018

2022

2026

1 0 2014

2030

4

2022

2026

2030

2022

2026

2030

7

x 10

x 10 2

5 4

1.5

3

L

I

2018

1

2 0.5

1 0 2014

2018

2022

2026

2030

Time (year)

0 2014

2018

Time (year)

Figure 11: Optimal controls with the weight coefficients B1 = (104 , 104 , 104 , 104 , 104 ) and state variables for Strategy II compared to the case of no controls from 2014 to 2030.

31

Strategy III u4(t)

1

0.5

0 2014

2016

2018

2020

2022

2024

2026

2028

2030

2024

2026

2028

2030

Strategy IV u5(t)

1

0.5

0 2014

2016

2018

2020

2022

Time (year)

0.014

0.012

0.01

37.23%

(E+I)/N

0.008

44.28% 0.006 57.45% without control Strategy I : with five optimal controls Strategy II : with optimal controls u and u

0.004

4

5

Strategy III : with optimal control u

4

0.002

Strategy IV : with optimal control u

5

98.87%

0 2014

2016

2018

2020

2022

2024

2026

2028

2030

Time (year) 4

12

x 10

u1(t) u2(t)

10

u3(t) u4(t)

8

u5(t)

6 4 2 0

Strategy I

Strategy II

Strategy III

Strategy IV

Figure 12: First two graphs plot optimal control u4 in Strategy III and optimal control u5 in Strategy IV. The third graph indicates the ratios of infected individuals ( E + I )/N in all four strategies (using five optimal controls, using two optimal controls u4 and u5 , using one optimal control u4 , and using one optimal control u5 ) comparing with ( E + I )/N in the case without controls, for the period from 2014 to 2030. Reduction rates of the values in 2030 are indicated for each strategy. The last graph depicts the relative cost of controls in each strategy with weight constant Bi = 104 for i = 1, 2, 3, 4, 5.

32