Optimal Control of a Model of Gambiense Sleeping Sickness in Humans and Cattle
Ndondo Mboma Apollinaire^{1}, Walo Omana Rebecca^{1}, Maurice Yengo Vala-ki-sisa^{2}
^{1}Faculty of Sciences, Regional Center for Doctoral Education in Mathematics and Computer Science, University of Kinshasa, Kinshasa, D. R. Congo
^{2}Faculty of Sciences, Department of Mathematics and Computer Science, University of Kinshasa, Kinshasa, D. R. Congo
Email address:
To cite this article:
Ndondo Mboma Apollinaire, Walo Omana Rebecca, Maurice Yengo Vala-ki-sisa. Optimal Control of a Model of Gambiense Sleeping Sickness in Humans and Cattle. American Journal of Applied Mathematics. Vol. 4, No. 5, 2016, pp. 204-216. doi: 10.11648/j.ajam.20160405.12
Received: July 5, 2016; Accepted: July 18, 2016; Published: August 31, 2016
Abstract: Human African trypanosomiasis (HAT) generally known as sleeping sickness is a fatal parasitic disease which appears mostly in sub-Saharan Africa, threatening millions of people and animals. Sleep disorders are a major feature of the (most) advanced stage of the disease, when the central nervous system is affected. In the absence of treatment, the outcome is always fatal. The parasite is transmitted to humans or animals through the bite of a tsetse fly previously infected by humans or animals carrying the parasite. We look for different scenarios to control the epidemic by integrating in our model terms that model the different control techniques.
Keywords: Trypanosoma Brucei Gambiense, Sleeping Sickness, Glossina, Optimization, Control, Modeling, Optimal Control
1. Introduction
Recently, a number of studies were carried out to explore the optimal control theory in some mathematical models for infectious diseases including the HIV virus [45,46], tuberculosis [48] and vector-borne diseases [47]. Authors in [47] derive the optimal control efforts for treatment and prevention in order to prevent the spread of a vector-borne disease using a system of ordinary differential equations (ODEs) for the host and vector populations. Authors in [9] investigated such optimal strategies for prevention, treatment and vector control using two systems of ODEs which consist of a stage-structured model for the vector and a SI/SIR-type model for the vector/host population. In this paper, using models described in [39] for the tsetse fly population dynamics and the transmission virus, we formulate the associated control model in order to derive an optimal prevention and treatment strategies with minimal implementation costs. Controls used here are based on five main actions applied in the epidemics. There are optimal strategies for prevention, treatment and vector control for humans and cattle using two systems of ODEs which consist of a stage-structured model for the vector and a SEI/SEIRS-type model for the vector/host populations (human and cattle).
The paper is organized as follows. In section 2, we present the compartmental models used in [39] to describe the tsetse fly population dynamics and HAT transmission between tsetse flies, humans and cattle. In section 3, we formulate an optimal control problem; first, we investigate the existence of an optimal control, then we derive the optimality system which characterizes the optimal control using Pontryagin’s Maximum Principle [49]. In section 4 numerical results illustrate our theoretical results.
2. The Basic Model
2.1. Modeling the Growth Dynamics of the Tsetse Fly
We first model the growth dynamics of the tsetse fly. From the life cycle of the tsetse fly, it is sufficient to consider two life stages, namely the pupae and the adult flies.
Let L (t) be the number of pupae at time t and A (t) be the number of (male and female) adult flies at time t. The dynamics of L and A is modeled by the following system
(1)
(2)
Here, is the rate at which female flies give birth to larvae; W is the proportion of female flies in the population of adult flies; is the pupae carrying the capacity of the nesting site; is the transfer rate from pupae to adult tsetse flies, so is the average time as a pupa; d_{L} and d_{F} are the mortality rate of pupae and adult flies, respectively; with all parameters assumed positive.
2.1.1. Equilibrium Points
The threshold defined by
(3)
is important when calculating the equilibrium points of system (1)–(2), as shown in the following result. Parameter r can be interpreted as the probability of surviving the pupal stage multiplied by the birthrate divided by the death rate. If r <1, the population of tsetse flies will be extinguished, otherwise they evolve toward an equilibrium given by
(4)
Details as regards the stability of this model are given in [39].
2.1.2. Formulation of the Full Model
For our full model, we assume that r >1 and that the flies are in equilibrium Trypanosomiasis in the fly population is modeled by an SEI compartmental model. It is assumed that a fly once infected will never recover or be removed. So we subdivide the adult fly population into three compartments, S_{F} susceptible tsetse flies, E_{F} exposed tsetse flies infected but not yet infectious, and I_{F} infectious tsetse flies that are able to transmit the disease once they bite a susceptible host. Thus the total adult fly population is
S_{F }+ E_{F }+ I_{F} (5)
The human and cattle host populations are described by a Malthus model. We denote by N_{H} and N_{C} the total size of the human and cattle host populations, respectively, at time t and b_{H}, b_{C}, d_{H}, d_{C} are the rates of birth and mortality of the human and cattle host populations, respectively. The dynamics of N_{H }and N_{C}_{ }is governed by
(6)
(7)
Where α_{H} = b_{H} − d_{H} and α_{C }= b_{C }−_{ }d_{C} are the growth rates of the human and cattle population respectively. If α_{H}<0 (α_{C}<0), the human (cattle) population will be extinguished, it will remain constant if α_{H}=0 (α_{C}=0), and will grow exponentially if α_{H}>0 (α_{C}>0). We assume that α_{H}=0 (α_{C}=0), i.e. b_{H}=d_{H} (b_{C}=d_{C}), so that the human (cattle) population is constant over the period of the study and that there is no human and cattle death due to HAT. Trypanosomiasis in the human and cattle host populations is modeled by SEIRS compartmental models, each with four compartments:
• susceptible hosts S_{H} (S_{C}): humans (cattle) at risk and disease free;
• Exposed hosts E_{H} (E_{C}): humans (cattle) in the latent stage of the disease. they are infected but unable to transmit the disease;
• Infectious hosts I_{H} (I_{C}), humans (cattle) able to transmit the disease to tsetse flies if they are bitten. These compartments contain hosts in the first stage of the disease with only minor symptoms or not aware if they are infected;
• Removed hosts R_{H} (R_{C}) consist of humans (cattle) in the second stage of the disease, very sick and not exposed to flies, so that they do not pass on infection, as well as humans (cattle) undergoing treatment and not exposed to flies. We assume that treatment starts at the beginning of stage 2, since this is usually when hosts become symptomatic. These compartments also contain removed humans (cattle) that have developed temporary immunity after recovery from stage 2 or treatment and they can neither transmit nor acquire HAT, but they will become susceptible again after the period of temporary immunity has passed.
The constant total human and cattle populations are defined by:
(8)
(9)
The dynamics of T. b. Gambiense in the tsetse fly population, assuming that transmission to flies occurs from humans and cattle in only the first stage of HAT, is given by the system
(10)
(11)
(12)
Here a is the vector blood feeding rate,is the probability that a fly becomes infected after biting an infectious human, v is the probability that a fly becomes infected after biting infectious cattle, is the incubation period in the fly, is the natural mortality rate of adult flies and is the proportion of tsetse fly bites on cattle (thus (1−p) is the proportion of bites on humans). This proportion is assumed to be constant as in Funk et al. [40]. For a discussion of this assumption see Rock et al. [10, Section 3.3].
The dynamics of T. b. gambiense in the human host population is governed by the system
Where is the probability that an infectious fly infects a human host, is the birth rate of the human population, is the human mortality rate, is the average incubation period for a human host, is the average length of stage 1 for humans corresponding to the infectious period. For untreated humans, is the sum of the average length of stage 2 and the average temporary immunity period. For treated humans, is the sum of the average length of treatment and the average temporary immunity period. Note that we assume that the average length of treatment is equal to the average length of stage 2. Similarly, the dynamics of T. b. Gambiense in the cattle host population is governed by the system
where u is the probability that an infectious fly infects a cattle host, is the birth rate of the cattle population, = is the cattle mortality rate, is the average incubation period for cattle, is the average length of stage 1 for cattle corresponding to the infectious period. For the untreated cattle, is the sum of the average length of stage 2 and the average temporary immunity period. For the cattle treated, is the sum of the average length of treatment and the average temporary immunity period. As for humans, we assume that the average length of treatment is equal to the average length of stage 2 for cattle.
Thus, the dynamics of the transmission of sleeping sickness is then described by the system of equations below, where we have assumed that there is no death due to the disease, no vertical transmission, and all parameters are positive, except and which are non negative.
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
(23)
Where for i ∈ {H, C} and
Figure 1 shows a flow diagram for this system and Table 1 describes the model parameters. Note that all cross transmission terms are normalized with respect to the host population as it is common in vector-borne disease models [19, Section 14.3]. Non negative initial conditions with positive and small complete the formulation of our HAT model in the invariant region
(24)
2.2. Model Equilibrium and Stability
In [39] A. M. Ndondo et al. showed that system (13)–(23) always reaches the disease free equilibrium (DFE), . Considering the local stability of , they followed the notation of [7] and only considered the infected compartments given by (13)−(18) and the calculation of the basic reproduction number gave
3. Optimal Control
In this section, we discuss the control measures to limit the impact of the tsetse fly in the transmission of HAT. We consider model (13)–(23) with different control measures:
• The first control measure noted models efforts to prevent the spread of infection. It brings together the various strategies used in order to reduce the number of vector-human contacts. We can cite in particular the use of repellents, mosquito nets (since the tsetse flies are mostly active in the late afternoon and early morning) or suggest that humans wear long-sleeved clothing to cover the exposed parts of the body, and prevent cattle from being bitten by dipping them in fly repellent solutions. Moreover, control µ1 also takes into account the municipality efforts to raise awareness of the danger of the tsetse fly and HAT in the population. Our action is limited in the interval [0, T].
• The second control measure is the treatment of human patients on the interval [0, T]. It also models the isolation of patients in treatment places (hospitals and isolation areas) to avoid all possible new contamination cases.
• The third control measurespecifically models the vector control on the interval [0, T]. It includes actions under the responsibility of vector control services and is intended to kill the pupae through the use of larvicides whose action is specifically directed against the pupae. It also includes the community effort to destroy potential breeding sites around homes and deprive the tsetse flies of the breeding sites in which their larvae can develop (in hot weather, at midday, tsetse hide in shelters to avoid high temperatures). Houses and huts also provide a fresh environment and many tsetse flies find refuge there in the dry season. Tsetse flies normally leave their refuge in the late afternoon to find blood meal [15], thus to avoid being bitten or to prevent tsetse flies from entering or facilitating their exit from habitations. The control measure can also model the efforts to trap tsetse flies, reducing the larvae capacity . There is also the idea of introducing a sterile population of flies to reduce the number of clutches (since the females are in general fertilized only once).
• The fourth control noted models efforts to prevent the spread of infection. It brings together the various strategies used in order to reduce the number of cattle-vector contacts. It brings together the different strategies used in order to reduce the number of contacts between vectors and cattle on the interval [0, T].
• The fifth control noted measures the treatment of infected cattle on the interval [0, T]. It also models the isolation of infected cattle in treatment places (isolated areas) to avoid all possible new contamination cases.
We will consider a situation where:
is greater than 1, and we are not at a stable disease free equilibrium??, and there is need to apply a treatment to control the epidemic. Under these conditions, the optimal control problem is formulated as follows:
(25)
(26)
(27)
(28)
(29)
(30)
(31)
(32)
(33)
(34)
(35)
(36)
(37)
Where
∈ [0, 1] corresponds to prevention. So if = 1, then the number of human-vector contacts is zero and if = 0, the infection rate is maximum and equal to or;
∈ [0, 1] corresponds to the treatment
is the proportion of actual treatment (so gives the proportion of humans cured with treatment)
∈ [0, 1] corresponds to the vector control. is the mortality rate of tsetse flies due to the use of chemicals or destruction of breeding sites.
∈ [0, 1] corresponds to prevention. So if = 1, then the number of cattle-vector contacts is zero and if = 0, the infection rate is maximum and equal to or;
∈ [0, 1] corresponds to the cattle treatment and is the proportion of actual treatment (so gives the proportion of cattle cured with treatment);
Proposition 3 D × Γ is positively invariant for system (17-29)
Proof
First we have:
Using Grönwall’s inequality, we deduce easily that while all variables in this system are positive, we have moreover
The right side of these inequalities is exactly the HAT transmission model (13-23) without control. Since the solutions of system (13-23) were defined in D × Γ. So, as done before, using Grönwall’s inequality, we deduce that the solutions of the last system are bounded.
We associate to system (27-37), the optimal control problem coupled with the following cost function:
(38)
The first terms represent the costs for populations, and that we want to reduce. The constants , , and are positive and correspond to the ratio used to regulate the control for the prevention, treatment and vector control in human and cattle population respectively. As given in the literature, the cost functions are assumed to be quadratic functions. In fact, a quadratic cost function is a natural way that allows the analogy with the expended energy for all those dimensions.
The objective is to limit the spread of the disease by reducing the number of flies and infected humans. We then look for the control terms (, , , , ) that minimize the cost:
where
with a piecewise continuous function on [0, T]
and i = 1, 2, 3, 4, 5
The aim is not only to reduce the populations of infected humans and infected vectors , and after a time T but also over [0, T] to act simultaneously on prevention (contact vectors-infected humans/cattle), treatment (infected humans/cattle), , and the vector control. The First, second and third terms of the functional J model the human population, the cattle population and the accumulated infected vectors from the initial time T_{0} = 0 to the final time = T.
The fourth, fifth, sixth, seventh and eighth terms model the cumulative costs of prevention, treatment and vector control. All of them vary according to time.
The choice of positive parameters A_{1}, A_{2}, A_{3}, B_{1}, B_{2}, B_{3}, B_{4 }and B_{5} depends on a subjective assessment of the relative importance given by the technical staff to reduce the population of infected human, infected cattle, infected flies by prevention, treatment of infected humans and cattle, and vector control (the weight or weighting for regulating the control). Λ is the set of controls and are constants belonging to [0, 1], i = 1, 2, 3, 4, 5
The optimal control problem is solved when we determine (Λ that minimizes the function (38). Hence our work will consist of the following steps:
• to show the existence of an optimal control;
• to give a characterization of the optimal control;
• to obtain numerical representations
3.1. Existence of the Optimal Control
Proposition 3.1 Consider the control problem associated with problem (25-37). There is a control (and a corresponding solution
which minimizes J (µ_{1}, µ_{2}, µ_{3}, µ_{4}, µ_{5}) on Λ such that
Proof
According to W. H. Fleming [1], Theorem 3.4.1, we need to check the following conditions:
1. The set of controls and corresponding solutions is not empty;
2. The set of controls of Λ is convex and closed in L^{2} (0, T);
3. The vector set of the state system is bounded by a linear function of control;
4. the integrand of the objective function is convex;
5. there exist constants and β > 0 such that the objective function integrand is bounded by
We check that µ_{1 }= µ_{2 }= µ_{3 }= µ_{4 }= µ_{5 }= 0 is a control in Λ and is a solution corresponding to control µ_{1 }= µ_{2 }= µ_{3 }= µ_{4 }= µ_{5 }= 0, so the set of all controls and corresponding solutions is not empty, which satisfies condition 1.
• The set Λ is bounded by definition, so condition 2 is satisfied.
• The system vector field (25-37) satisfies condition 3 because it is bounded.
• there exist c_{1}, c_{2 }>0 and β>1 satisfying
Since variables are bounded, we deduce then the existence of an optimal control which minimizes the cost function J (µ_{1}, µ_{2}, µ_{3}, µ_{4}, µ_{5}).
In summary, for this minimization problem, the necessary convexity condition required for the functional J with parameters I_{H}, I_{F}, I_{C}, µ_{1}, µ_{2}, µ_{3}, µ_{4} and µ_{5} holds. The right set of the system of equations (25-37) is linearly bounded because of the fact that a priori variable T is bounded, which implies that other state variables are also bounded. Boundedness and the fact that the bounds are finite guarantee the required compactness for the existence of optimal control. The initial conditions are .
Since now we have assured the existence of an optimal control, we can use Pontryagin’s Maximum principle to solve the optimal control problem.
3.2. Characterization of Optimal Control
Theorem 4 The optimal control which minimizes the functional J given in (38) under the constraints given by the system of differential equations (25-37) is given by:
Proof
Let and Π = (λ_{1}, λ_{2}, λ_{3}, λ_{4}, λ_{5}, λ_{6}, λ_{7}, λ_{8}, λ_{9}, λ_{10}, λ_{11}, λ_{12}, λ_{13}) be the adjoint variables. We define the Lagrangian (Hamiltonian increased penalties) associated with the problem defined above:
Where ω_{1}, ω_{2}, ω_{3}, ω_{4}, ω_{5}, ω_{6}, ω_{7}, ω_{8}, ω_{9 }and ω_{10 }are penalty variables (multipliers) attached to control µ_{1}, µ_{2}, µ_{3}, µ_{4 }and µ_{5}. These penalty multipliers must meet the following conditions:
ω_{1}(t)(µ_{1}(t) −a_{1}) = 0; ω_{2}(t)(b_{1 }−µ_{1}(t)) = 0
ω_{3}(t)(µ_{2}(t) −a_{2}) = 0; ω_{4}(t)(b_{2 }−µ_{2}(t)) = 0
ω_{5}(t)(µ_{3}(t) −a_{3}) = 0; ω_{6}(t)(b_{3 }−µ_{3}(t)) = 0
ω_{7}(t)(µ_{4}(t) −a_{4}) = 0; ω_{8}(t)(b_{4 }−µ_{4}(t)) = 0
ω_{9}(t)(µ_{5}(t) −a_{5}) = 0; ω_{10}(t)(b_{5 }−µ_{5}(t)) = 0
In addition, the differential equations which govern the adjoint variables are obtained by differentiating the Lagrangian (as per Maximum Principle):
For these adjoint variables, we must have λi (T) = 0, i = 1, 2,···, 9, these are the transversality conditions called also marriage conditions.
The value of the optimal control can be characterized at each time t ∈ [0, T] noting that it minimizes the Lagrangian (Pontryagin’s Maximum Principle) and that’s at this optimal control, variables must satisfy the necessary condition:
.
Given that
Differentiating with respect to µ_{1}, µ_{2}, µ_{3}, µ_{4} and µ_{5 }gives respectively:
At , we have
These equalities give:
Hence, we get the optimum
For a more explicit formula for optimal controls without ω_{1}, ω_{2}, ω_{3}, ω_{4}, ω_{5 }and ω_{6}, we use standard techniques. For this purpose, we consider 15 cases, with 3 cases for each optimal control:
1) Case 1
In the set
So the optimal control is
2) Case 2
In the set
So the optimal control is
With
3) Case 3
In the set from where
and therefore
With
4) Case 4
In the set
So the optimal control is
5) Case 5
In the set , hence
.
.
6) Case 6
In the set , hence
.
.
With these three cases, we rewrite the expression of the second control:
7) Case 7
In the set
So the optimal control is
8) Case 8
In the set , hence
.
. .
9) Case 9
In the set , hence
.
. .
With these nine cases, we rewrite the expression of the third control:
10) Case 10
In the set
So the optimal control is
11) Case 11
In the set, hence
s
So the optimal control is
12) Case 12
In the set , hence
s
So the optimal control is
With these three cases, we rewrite the expression of the third control:
13) Case 13
In the set
So the optimal control is
14) Case 14
In the set, hence
.
15) Case 15
In the set , hence
.
.
With these three cases, we rewrite the expression of the third control:
In summary, we obtain the control (, , , , ) that optimizes the functional J under the constraints given by the system of differential equations (25-37), its expression is given by the optimal quintuplet defined by:
This completes the proof of the Theorem.
4. Numerical Results
Figures 2, 3, 4 and 5 are plots with different values of a and p. The value a = 0.25 means that flies take one blood meal in four days and the value a = 1 means flies take one blood meal very day. In the other hand, the value p = 0.7 shows the preference of blood meal taken by flies on cattle as given below. We assume the carrying capacity of the tsetse pupae K_{L }= 300000, the human population N_{H} = 300, and the cattle population N_{C} = 50. Baseline parameter values given in Table 1 were collected from the literature on HAT in West Africa as cited, and values that were not found in the literature were estimated. Values from Table 1 gives r = 1.0154, the number of larvae 4545, and the number of adult flies = 5000. Note that since by our assumptions compartments R_{H }and R_{C }contain hosts in stage 2 (or in treatment) and recovered hosts, = 30 + 90 days and = 25 + 75 days have the same values as given by Rogers [35] for the sums of the duration of infection and immunity in species 1 and 2, although the definitions of our parameters are different. Stage 1 of T. b. gambiense HAT in humans in Africa may last for several months [2] (i.e. may be much smaller than the above value). Thus our baseline values apply more to our model with treatment giving control reproduction numbers.
5. Concluding
We have built a model that includes the most important mechanisms of transmission of sleeping sickness between a host population and a vector population. First fly growth dynamics is modeled. We use it because unlike most flies, once fertilized the tsetse fly generally keeps its eggs and each birth, one larva is hatched. In addition, the tsetse flies of both sexes are blood-sucking insects. A SEI model was used to express the different states of infection of the tsetse fly. For the host population, we modeled the dynamics of a SEIRS model. A coupled model was obtained and analyzed. A mathematical analysis of the model was made. The equilibrium points of the model were calculated and presented. The study of the system’s stability of these equilibrium points was presented. A numerical simulation of the system illustrates the theoretical results obtained and discussed throughout this study.
We built optimal control problem to assess the efficiency of three control measures to reduce the number of infections in the human population. The study of the proposed model allows us to know when the controls efforts should be applied, and the importance to be given to each of these vector control efforts, the treatment to administer to patients infected by the sleeping disease and the various means of prevention to avoid or reduce the host-vector contact. These efforts can obviously be applied to a finite time interval [0, T]. The model does not allow us to make long-term predictions of the disease dynamics. We have shown the existence of optimal controls for the cost minimization problem (quadratic) to reduce the densities of infected individuals in the population, namely the individuals who constitute the germ reservoir and the essential agents in the spread of infection.
References