Applied Mathematical Sciences, Vol. 3, 2009, no. 5, 231 - 240
Optimal Control of Tuberculosis with
Exogenous Reinfection
K. Hattaf∗ , M. Rachik, S. Saadi, Y. Tabit and N. Yousfi†
Laboratoire d’Analyse Modelisation et Simulation,
Departement de Mathematiques et d’Informatique,
Faculte des Sciences Ben M’Sik, Universit Hassan II
Mohammedia, B.P 7955, Sidi Othman, Casablanca, Maroc
Abstract
The aim of this work is the application of optimal control for the
system of ordinary differential equations modeling a tuberculosis disease
with exogenous reinfection. Seeking to reduce the infectious group by
the reduction of the contact between infectious and exposed individuals,
we use control representing the prevention of exogenous reinfection. The
Pontryagin’s maximum principle is used to characterize the optimal
control. The optimality system is derived and solved numerically.
Keywords: Optimal control, exogenous reinfection, Tuberculosis, TB model
with exogenous reinfection
1
Introduction
Tuberculosis (TB) is a bacterial disease with about one third of the world human population as its reservoir (see [1]and [9]) but only a small proportion (approximately 10%) of individuals develop the progressive disease (active TB).
Most people are assumed to mount an effective immune response to the initial
infection that limits proliferation of the bacilli and leads to long-lasting partial immunity both to further infection and to the reactivation of latent bacilli
remaining from the original infection (see [10]). Individuals who have a latent
infection are not clinically ill or capable of transmitting TB (see [9]). Exposed
individuals may remain in this latent stage for long and variable periods of
time (in fact, many die without ever developing active TB). At greater ages,
∗
†
Corresponding author. Email: k.hattaf@yahoo.fr
Corresponding author. Email: nourayousfi@gmail.com
232
N. Yousfi et al
the immunity of persons who have been previously infected may wane, and
they may be then at risk of developing active TB as a consequence of either
exogenous reinfection (i.e., acquiring a new infection from another infectious
individual) or endogenous reactivation of latent bacilli (i.e., re-activation of a
pre-existing dormant infection) (see [11] and [10]).
The exogenous reinfection plays a key role in tuberculosis transmission in
high-incidence regions, particularly in Africa (where HIV is high) and in inner
cities of developed countries.
In literature, Zhilan Feng, Carlos Castillo-Chavez and Angel F. Capurro
(see [3]) have incorporated exogenous reinfection into a epidemiological model
for the transmission dynamics of tuberculosis and have discussed control of
the disease by looking at the role of disease transmission parameters in the
reduction of R0 and the prevalence of the disease. However, this model did
not account for time dependent control strategies since their discussions are
based on prevalence of the disease at equilibria.
In this article we consider optimal control strategies associated with preventing exogenous reinfection based on a exogenous reinfection tuberculosis
model developed in [3]. This model is governed by the following system of
ordinary differential equations
dS
dt
dE
dt
dI
dt
dT
dt
I
,
N
I
I
I
= βcS − pβcE − (µ + k)E + σβcT ,
N
N
N
I
= pβcE + kE − (µ + d + r)I,
N
I
= rI − µT − σβcT ,
N
= Λ − µS − βcS
(1)
(2)
(3)
(4)
where the four epidemiological classes and parameters model are defined in
section 2.
The term pβcE NI models exogenous reinfection, that is, the potential reactivation of TB by continuous exposure of latently-infected individuals to those
who have active infections, and p represents the level of reinfection. When
p = 0, system (1)-(4) reduces to our earlier TB model (Castillo-Chavez and
Feng, 1997). A value of p ∈ [0, 1] implies that reinfection is less likely than a
new infection. In fact, a value of p ∈ [0, 1] implies that a primary infection
provides some degree of cross immunity to exogenous reinfections. A value of
p ∈ [1, +∞[ implies that TB infection increases the likelihood of active TB.
We introduce into this model a control to the above mentioned model simulating effect of exogenous reinfection. This exogenous reinfection effect is
incorporated by adding a term that may lower the exogenous reinfection rate
reducing the contact between infectious and exposed individuals so that the
233
Control of tuberculosis
number of infectious individuals with tuberculosis will be reduced. Our objective functional balances the effect of minimizing the cases of infectious TB
and minimizing the cost of implementing the control.
The paper is organized as follows. Section 2 describes a model with exogenous reinfection with control variable. The analysis of optimization problems
is presented in section 3. In section 4, we present a numerical appropriate
method and the simulation corresponding results. Finally, the conclusion are
summarized in Section 5.
2
A TB model with exogenous reinfection
In this section, we present a TB model with exogenous reinfection [3]. The host
population is divided into the four epidemiological classes: namely susceptible
S(t), exposed (infected but not infectious) E(t), infectious I(t) and treated
(removed) T (t). We assume that an individual can be infected only through
contacts with infectious individuals. Our TB model with exogenous reinfection
is given by the following nonlinear system of differential equations
dS
dt
dE
dt
dI
dt
dT
dt
= Λ − µS − βcS
I
,
N
I
I
I
− pβc(1 − u)E − (µ + k)E + σβcT ,
N
N
N
I
= pβc(1 − u)E + kE − (µ + d + r)I,
N
I
= rI − µT − σβcT ,
N
= βcS
(5)
(6)
(7)
(8)
where S(0) = S0 , E(0) = E0 , I(0) = I0 , T (0) = T0 are given and the definitions
of above model parameters are listed in Tab. 1.
The control variables, u(t), is bounded, Lebesgue integrable function [7].
The coefficient, 1 − u(t), represents the effort that prevents the exogenous
reinfection in order to reduce the contact between the infectious and exposed
individuals, then we decease the number of infectious individuals.
If u = 1, the prevention of exogenous reinfection is 100% effective, whereas if
u = 0, we find the model for tuberculosis with exogenous reinfection (see, [3]).
234
N. Yousfi et al
Parameter
Λ
µ
d
β
σβ,
0≤σ≤1
c
k
r
N
p
Definition
Recruitment rate
Natural mortality rate
TB induced mortality rate
Average number of susceptible individual infected by
one infectious individual per contact per unit of time
Average number of treated individual infected by
one infectious individual per contact per unit of time
Per-capita contact rate
Rate of progression to active TB
Per-capita treatment rate
Total population N = S + E + I + T
level of exogenous reinfection
Table 1: Parameter definitions
3
The optimal control problems
The problem is to minimize the objective functional
tf
[I(t) + Au2 (t)]dt
J(u) =
(9)
t0
where the parameter A represents the weight on the benefit and cost ( A balance the size of the terms). Our target is to minimize the objective functional
defined in equation (9) by minimizing the number of the infectious classes. In
other words, we are seeking optimal control u∗ such that
J(u∗ ) = min{J(u) : u ∈ U},
(10)
where U is the control set defined by
U = {u ∈ L1 (0, tf ) : 0 ≤ u ≤ 1}.
Pontryagin’s Maximum Principle[12] provides necessary conditions for an optimal control problem. This principle converts (5) - (8), (9) and (10) into a
problem of minimizing an Hamiltonian, H, pointwisely with respect to u :
2
H = I(t) + Au (t) +
4
λi fi ,
(11)
i=1
where fi is the right hand side of the differential equation of i-th state variable.
By applying Pontryagin’s Maximum Principle [12] and the existence result for
the optimal control from [4], we obtain the following theorem.
235
Control of tuberculosis
Theorem 3.1 There exists an optimal control u∗ (t) and corresponding solution S ∗ , E ∗ , I ∗ , T ∗ and J ∗ , that minimizes J(u) over U. Furthermore, there
exists adjoint functions, λ1 , λ2 , λ3 , λ4 satisfying the equations
′
λ1 = λ1 µ + (λ1 − λ2 )βc
I∗
,
N
′
λ2 = λ2 µ + (λ2 − λ3 )(k + p(1 − u)βc
′
λ3 = −1 + λ3 (µ + d) + (λ1 − λ2 )βc
I∗
),
N
S∗
E∗
+ (λ2 − λ3 )p(1 − u)βc
N
N
T∗
+(λ4 − λ2 )σβc
+ (λ3 − λ4 )r,
N
I∗
′
λ4 = λ4 µ + (λ4 − λ2 )σβc ,
N
with transversality conditions
λi (tf ) = 0, i = 1, ..., 4.
Moreover, the optimal control is given by
1
(λ3 − λ2 )pβcE ∗ I ∗ ))
(12)
u∗ = min(1, max(0,
2AN
Proof.
Due to the convexity of integrand of J with respect to u, a priori boundedness of
the state solutions, and the Lipschitz property of the state system with respect
to the state variables. The existence of an optimal control has been given by
[4] (see Corollary 4.1). The adjoint equations and transversality conditions
can be obtained by using Pontryagin’s Maximum Principle such that
∂H
, λ1 (tf ) = 0,
∂S
∂H
′
, λ2 (tf ) = 0,
λ2 = −
∂E
∂H
′
λ3 = −
, λ3 (tf ) = 0,
∂I
∂H
′
, λ4 (tf ) = 0.
λ4 = −
∂T
The optimal control u can be solve from the optimality condition,
′
λ1 = −
∂H
= 0.
∂u
That is
∂H
= 2Au + (λ2 − λ3 )pβcE ∗ I ∗ /N = 0.
∂u
By the bounds in U of the controls, it is easy to obtain u∗ in the form of
(12).
236
4
N. Yousfi et al
Numerical simulations
The numerical algorithm presented below is a semi-implicit finite difference
method.
We discretize the interval [t0 , tf ] at the points ti = t0 +ih (i = 0, 1, ..., n), where
h is the time step such that tn = tf , [5]. Next, we define the state and adjoint
variables S(t), E(t), I(t), T(t), λ1 (t), λ2 (t), λ3 (t), λ4 (t), and the control u in
terms of nodal points Si , Ei , Ii , Ti , λi1 , λi2 , λi3 , λi4 and ui. Now a combination
of forward and backward difference approximation is used as follows :
The Method, developed by [6] and presented in [8], is then read as :
Si+1 − Si
Ii
= Λ − µSi+1 − βcSi+1 ,
h
N
Ii
Ii
Ii
Ei+1 − Ei
= βcSi+1 − p(1 − ui)βcEi+1 − (µ + k)Ei+1 + σβcTi ,
h
N
N
N
Ii+1 − Ii
Ii+1
= p(1 − ui )βcEi+1
+ kEi+1 − (µ + d + r)Ii+1 ,
h
N
Ii+1
Ti+1 − Ti
= rIi+1 − µTi+1 − σβcTi+1
.
h
N
By using a similar technique, we approximate the time derivative of the adjoint
variables by their first- order backward-difference and we use the appropriated
scheme as follows
λ1n−i − λ1n−i−1
Ii+1
= λn−i−1
µ + (λn−i−1
− λ2n−i )βc
,
1
1
h
N
λ2n−i − λ2n−i−1
Ii+1
= λn−i−1
),
µ + (λn−i−1
− λ3n−i )(k + p(1 − ui )βc
2
2
h
N
Si+1
λ3n−i − λ3n−i−1
)βc
− λn−i−1
= −1 + λ3n−i−1 (µ + d) + (λn−i−1
2
1
h
N
Ei+1
n−i−1
n−i−1
+(λ2
− λ3
)p(1 − ui )βc
N
Ti+1
n−i
n−i−1
n−i−1
+(λ4 − λ2
)σβc
− λ4n−i )r,
+ (λ3
N
Ii+1
λ4n−i − λ4n−i−1
= λn−i−1
.
µ + (λn−i−1
− λn−i−1
)σβc
4
4
2
h
N
The algorithm describing the approximation method for obtaining the optimal
control is the following
Algorithm
step 1 :
S(0) = S0 , E(0) = E0 , I(0) = I0 , T (0) = T0 , λi (tf ) = 0 (i=1, ..., 4),
u(0) = 0.
237
Control of tuberculosis
step 2 :
for i=1, ..., n-1, do :
Si+1 =
Ei+1 =
Ii+1 =
Si + hΛ
Ii
)
1 + h(µ + βc N
Ii
Ii
+ σβcTi N
)
Ei + h(βcSi+1 N
Ii
1 + h(µ + k + p(1 − ui)βc N )
Ii + hkEi+1
i+1
)
1 + h(µ + d + r − p(1 − ui)βc EN
Ti+1 =
λ1n−i−1
λ2n−i−1 =
Ti + hrIi+1
)
1 + h(µ + σβc Ii+1
N
=
λ1n−i + hλ2n−i βc Ii+1
N
)
1 + h(µ + βc Ii+1
N
λ2n−i + hλ3n−i (k + p(1 − ui)βc Ii+1
)
N
1 + h(µ + k + p(1 − ui)βc Ii+1
)
N
λn−i−1
= [λ3n−i + h(1 + (λn−i−1
− λn−i−1
)βc
3
2
1
−λn−i−1
p(1 − ui )βc
2
Ei+1
Ti+1
+ (λn−i−1
− λ4n−i )σβc
2
N
N
+λ4n−i r)]/[1 + h(µ + d + r − p(1 − ui )βc
λ4n−i−1 =
Ri+1 =
Si+1
N
σβc Ii+1
λ4n−i + hλn−i−1
2
N
1 + h(µ + σβc Ii+1
)
N
Ei+1
)]
N
,
(λn−i−1
− λn−i−1
)r1 Ei+1 Ii+1
3
2
,
2AN
ui+1 = min(1, max(Ri+1 , 0))
end for
step 3 :
for i=1, ..., n-1, write
S ∗ (ti ) = Si , E ∗ (ti ) = Ei , I ∗ (ti ) = Ii , T ∗ (ti ) = Ti , u∗ (ti ) = ui .
end for
238
N. Yousfi et al
The following parameters and initial values are used for the simulation which
we have taken from [2], [3] and [7]:
µ = 0.016, Λ = µ ∗ N = 192, d = 0.1, k = 0.005, c = 1, β = 13, p = 0.4,
σ = 0.9, r = 2, N = 12000, S0 = 7600, E0 = 3800, I0 = 500 and T0 = 100.
the period of the prevention effort is 12 months, and we take A = 400. In
figure 1, we remark that in absence of prevention the number I (solid curve) of
individuals infectious with Tuberculosis increase in the first five months start
to grow after.
Whereas, in presence of prevention, the number I (dashed curve) of individuals
infectious decreases, in addition the number of individuals I infectious with
TB at the final time tf = 1 (years) is 256 in the case with control and 584
without control, and the total cases of TB prevented at the end of the control
program is 319 (= 584 − 256). Finally, the figure 1 represents the control
optimal u∗ for the effort that prevents the exogenous reinfection in order to
reduce the number of infectious individuals.
600
I without control
I with control
550
500
I(t)
450
400
350
300
250
0
0.1
0.2
0.3
0.4
0.5
0.6
Time (years)
0.7
0.8
0.9
1
Figure 1: Function I with and without control
5
Conclusion
Our numerical results show the effectiveness to introduce the control that prevents the exogenous reinfection which reactivates the bacterium tuberculosis
at the latent individuals. This objective is realizable by sensitizing the latent
individuals not to contact the infectious individuals with TB active , particularly in the closed places. As says the proverb ”to prevent better than to
cure”.
239
Control of tuberculosis
0.7
0.6
Control u
0.5
0.4
0.3
0.2
0.1
0
0
0.1
0.2
0.3
0.4
0.5
0.6
Time (years)
0.7
0.8
0.9
1
Figure 2: The control u
References
[1] B. R. Bloom, Tuberculosis: Pathogenesis, Protection, and Control, ASM
Press, Washington, D.C., 1994.
[2] C. Castillo-Chavez and Z. Feng, To treat or not to treat: the case of
tuberculosis, J. Mathematical Biology, 35 (1997), 629659.
[3] Z. Feng, C. Castillo-Chavez, A.F. Capurro, A model for tuberculosis with
exogenous reinfection, Theor. Pop. Biol.57, (2000), 235.
[4] W. H. Fleming, R. W. Rishel, Deterministic and Stochastic Optimal Control, Springer Verlag, New York, 1975.
[5] A. B. Gumel, K. C. Patidar, and R. J. Spiteri, editors, Asymptotically
Consistent Non-standard Finite-Difference Methods for Solving Mathematical Models Arising in Population Biology, R. E. Mickens and Worl
Scientific, Singapore, 2005.
[6] A. B. Gumel, P. N. Shivakumar, and B. M. Sahai, A mathematical model
for the dynamics of HIV-1 during the typical course of infection, Third
world congress of nonlinear analysts, (2001), 47:20732083.
[7] E. Jung, S. Lenhart, Z. Feng, Optimal control of treatments in a two-strain
tuberculosis model, Discrete and Continuous Dynamical Systems-Series
B, (2002), 2(4): 473482.
240
N. Yousfi et al
[8] J. Karrakchou, M. Rachik, and S. Gourari, Optimal control and Infectiology: Application to an HIV/AIDS Model, Applied Mathematics and
Computation, (2006), 177:807818.
[9] B. Miller, Preventive therapy for tuberculosis, Med. Clin. North Am. 77,
(1993), 1263-1275.
[10] P. G. Smith, A. R. Moss, Epidemiology of tuberculosis, in ”Tuberculosis:
Pathogenesis, Protection, and Control”, B. R. Bloom, Eds., ASM Press,
Washington, 1994.
[11] K. Styblo, Epidemiology of Tuberculosis, VEB Gustav Fischer Verlag
Jena, The Hague, 1991.
[12] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, E. F. Mishchenko,
The Mathematical Theory of Optimal Processes, Wiley, New York, 1962.
Received: July, 2008