Analysis and Optimal control of deterministic Vector-Borne diseases model

by

Hyun J. Son

A dissertation submitted to the Graduate Faculty of

Auburn University
in partial fulfillment of the
requirements for the Degree of
Doctor of Philosophy

Auburn, Alabama
August 4, 2018

Keywords: vertically transmitted, sensitivity, Vector-borne, optimal control, basic

reproduction number,

Copyright 2018 by Hyun J. Son

Approved by

Dr. Yanzhao Cao, Chair, Professor of Mathematics and Statistics

Dr. Xiaoying (Maggie) Han, Professor of Mathematics and Statistics
Dr. Junshan Lin, Assistant Professor of Mathematics and Statistics
Dr. Wenxian Shen, Professor of Mathematics and Statistics
 Abstract

In this dissertation, two systems of deterministic differential equations are introduced to

study the transmission of vector-borne diseases between a host and a vector. The total popu-
lation of the host and the vector are divided into different compartments. The total population
of host is divided into susceptible , exposed, infected, and treated groups. The total population
of vector is divided into susceptible, exposed, and infected. In chapter 2, we introduce a mo-
del to study vertically transmitted vector-borne diseases with nonlinear system of differential
equations. We analyze the model by finding the disease free equilibrium point E0 and deriving
the basic reproductive number R0 by using the next generation matrix method. We study the
local and global stability of E0 and how the stability is related to R0 . We study the sensitivity
of R0 using the normalized forward sensitivity index and find the relation to the parameters in
the model. We have numerical simulations to show the result we get from the analysis based on
the dengue virus. In chapter 3, we introduce a model to study optimal control to find the best
way to control viruses. We introduce two optimal controls, the prevention of contact between
host and vector u1 and the treatment of host u2 in the model given in the chapter 2. We consider
a cost functional related to the cost of the prevention and the treatment. We try to minimize
the number of exposed and infected host groups and maximize the number of susceptible and
treated host groups. We show the existence of u1 and u2 by using Carathodory’s existence theo-
rem. We find the explicit formula of u1 and u2 with the status variables and the adjoint variables
from Hamiltonian by using Pontryagin’s maximum principle. We find the numerical values of
u1 and u2 by solving the given status system and the adjoint system derived from Hamiltonian.
We use the forward-backward sweep method and the Runge-Kutta method in 3-dimension to
solve the status system and the adjoint system. In the numerical simulation, we compare the
result between the controlled case and uncontrolled case for each host and vector groups. Also
we see which control among u1 and u2 is more effective to control the virus. In appendices, we
show the Matlab code used in the numerical simulation.

ii
 Acknowledgments

I want to express my deepest appreciation and gratitude to my advisor Dr. Yanzhao Cao
who guided me with patience and supported during the Ph.D course at Auburn University. I
really thank for his effort, enthusiasm, and support through all the years. He continuously
encouraged and convinced me that I could finish the research to get the Ph.D and I can do
better on teaching. I could not write the dissertation without his guidance and persistent advise.
I also thank the committee Dr. Xiayoing Han, Dr. Junshan Lin, Dr. Wenxian Shen, and the
university reader Dr. Sang-Jin Suh who helped me to revise and improve the dissertation.
I also want to express my deepest appreciation to all members of the Department of Mathe-
matics and Statistics at Auburn University. The faculty at the department have been really kind
and available whenever I needed help and so I take the further steps to improve the research
and the teaching. I am really glad that I am a part of the department.
I want to say thanks to my My wife, Hwanhee Lee, and two kids, Jason and Justin. I
couldn’t have accomplished my study without their love, support, and immense understanding.
I also want to say thanks to my father, mother, and brother who longed for my graduation.
They always encouraged me to keep up the study and supported my life with everything they
have.
I thank my friends who shared the knowledge and the friendship and encouraged me in
the department and out of the department.
Above all, I thank the Almighty God who allow me to have the knowledge, the strength
and the guidance to complete the Ph.D. I can not imagine the completion and the new job
without his help.

iii
 Table of Contents

Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii

1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2 Analysis of Vertically transmitted vector-borne disease model . . . . . . . . . . . . . 11

2.1 Vertically transmitted vector-borne disease Model . . . . . . . . . . . . . . . . 11

2.2 Disease free equilibrium point E0 and R0 . . . . . . . . . . . . . . . . . . . . 15

2.3 Stability of disease free equilibrium point . . . . . . . . . . . . . . . . . . . . 18

2.4 Sensitivity analysis of R0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.5 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3 Optimal control for Vertically transmitted vector-borne disease:treament and preven-

tion of human-mosquito interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.1 A Model for Optimal Control of vertically transmitted vector-borne disease . . 33

3.2 Optimal Control Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.3 Existence of an Optimal Control . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.4 Optimality System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.5 The forward-backward sweep method Algorithm . . . . . . . . . . . . . . . . 44

3.6 The Runge-Kutta method in 3-dimension . . . . . . . . . . . . . . . . . . . . . 45

3.7 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

iv
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

.1 Matlab Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

v
 List of Figures

1.1 The general transfer diagram for the MSEIR model . . . . . . . . . . . . . . . 6

2.1 Compartment between human and vector where Sh = x1 , Eh = x2 , Ih = x3 ,

Th = x4 , Sv = y1 , Ev = y2 , and Iv = y3 . . . . . . . . . . . . . . . . . . . . . . 13

2.2 Solution for hosts, x1 , x2 , x3 , and x4 , with the values of the parameters in Table
2.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.3 Solution for vectors y1 , y2 , and y3 , with the values of the parameters in Table 2.4. 27
1 1 1
2.4 Solution for Exposed, x2 , and Infectious,x3 , host for γ = 10 , γ = 13 , γ = 15 ,
1
γ = 17 with the values of the parameters in Table 2.4. In each cases, R0 ≈
0.0328, R0 = 0.0502, R0 = 0.065, R0 = 0.0838, respectively. . . . . . . . . . . 28

2.5 Solution for Exposed, x2 , and Infectious,x3 , host for δ1 = 0.04, δ1 = 0.06,
1
δ1 = 0.08 with γ = 10 and the values of the parameters in Table 2.4. In each
cases, R0 ≈ 0.0328, R0 ≈ 0.0405, R0 ≈ 0.0571, respectively. . . . . . . . . . . 29

2.6 Solution for Exposed, x2 , and Infectious,x3 , host for φ = 2, φ = 3, φ = 4 with

1
δ1 = 0.05, γ = 17 , and the values of the parameters in Table 2.4. In each cases,
R0 ≈ 0.1223, R0 ≈ 0.1822, R0 ≈ 0.2421, respectively. . . . . . . . . . . . . . 30

2.7 Solution for Exposed, x2 , and Infectious,x3 , host for φ = 17, φ = 18, φ = 19
1
with δ1 = 0.05, γ = 17 , and the values of the parameters in Table 2.4. In each
cases, R0 ≈ 1.0209, R0 ≈ 1.0808, R0 ≈ 1.1407, respectively. . . . . . . . . . . 31

2.8 A phase plane portrait for the Exposed and the Infectious host,x2 + x3 , against
the Exposed and the Infectious vector, y2 + y3 , where φ = 17. . . . . . . . . . 32

3.1 Comparison between Controlled(left) and Uncontrolled(right) for susceptible

host . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.2 Comparison between Controlled(left) and Uncontrolled(right) for exposed host 49

3.3 Comparison between Controlled(left) and Uncontrolled(right) for infected host 49

3.4 Comparison between Controlled(left) and Uncontrolled(right) for treated host . 50

3.5 Comparison between Controlled(left) and Uncontrolled(right) for susceptible

vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

vi
3.6 Comparison between Controlled(left) and Uncontrolled(right) for exposed vector 51

3.7 Comparison between Controlled(left) and Uncontrolled(right) for infected vector 51

3.8 Optimal Control u1 (left) and u2 (right) . . . . . . . . . . . . . . . . . . . . . . 52

vii
 List of Tables

1.1 Average Laboratory vertical infection (v), vertical transmission (vt ), and filial
infection (ft ) rates for Aedes mosquitoes. [16] . . . . . . . . . . . . . . . . . . 3

2.1 Parameters and description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.2 Parameter estimation for dengue fever in [12, 13]. Rate is per day. . . . . . . . 23

2.3 Sensitivity index for parameters in R0 . The parameters are orders from most
1
sensitivity to [Link] are : θ1 = 0.0083, θ2 = 0.00513, µ = 70(365) ,
0.375 1 1
Λ = 70(365) , ρ = 205, δ0 = 1050, β = 0.2, r = 7 , δ1 = 0.0399, d = 10 ,
ψ = 0.00233, ζ2 = 0.0087, α = 0.0001, φ = 2,  = 18 , γ = 17
1
. . . . . . . . . . 24

2.4 Parameter values for a dengue virus. In this case, R0 ≈ 0.0838. . . . . . . . . . 25

3.1 Parameter values are estimated based on a dengue virus. . . . . . . . . . . . . . 48

viii
 Chapter 1

Introduction

In this thesis, we study mathematical models for the dynamics of vector-borne diseases, especi-
ally dengue virus. We will provide a brief explanation for vector-borne diseases, mathematical
modeling of vector-borne diseases, dengue virus, vertical transmission of dengue virus, analy-
sis of the model, and optimal control approach to find the best way to control the virus.
Vectors are living organisms that can transmit infectious diseases between humans or from
animals to humans. Many of these vectors are bloodsucking which ingest disease-producing
microorganisms during a blood meal from an infected host (human or animal) and later inject
it into a new host during their subsequent blood meal.
Vector-borne diseases are human illnesses caused by parasites, virus and bacteria that are
transmitted by mosquitoes, sandflies, triatomic bugs, blackflies, ticks, tsetse flies, mites, snails
and lice. Every year there are more than 700,000 deaths from diseases such as malaria, dengue,
schistosomiasis, human African trypanosomiasis, leishmaniasis, Chagas disease, yellow fever,
Japanese encephalitis and onchocerciasis globally. Since 2014, the major vector-born diseases
are dengue, malaria, chikungunya, yellow fever, and Zika.
Changes in agricultural practices due to variation in temperature and rainfall can affect
the transmission of vector-borne diseases. The growth of urban slums, lacking reliable piped
water or adequate solid waste management, can render large populations in towns and cities at
risk of viral diseases spread by mosquitoes. Together, such factors influence the reach of vector
populations and the transmission patterns of disease-causing pathogens. [72]
Dengue fever is a mosquito-borne viral infection and is a severe, flu-like illness that affects
infants, young children and adults [32, 71]. It occasionally develops into Severe dengue which

1
is the leading cause of serious illness and death among children. Classical dengue fever is
generally observed in older children and adults and is characterized by sudden onset of fever,
frontal headache, nausea, vomiting, and other symptoms. The actual illness last for 3 to 7 days
is usually benign.
It is transmitted by the Aedes aegypti mosquitoes and the Aedes albopictus mosquitoes.
The Aedes aegypti mosquito lives in urban habitats and breeds mostly in man-made containers.
Unlike other mosquitoes the Aedes aegypti mosquito is a day-time feeder. Therefore the peak
biting periods are early in the morning and in the evening before dark. The Aedes albopictus
mosquito lives mostly in Asia but due to the international trade, it has spread to North America
and more than 25 countries in Europe.
The virus is transmitted to humans through the bites of infected female mosquitoes. In-
fected symptomatic or asymptomatic humans are the main carriers and multipliers of the virus,
serving as a source of the virus for uninfected mosquitoes. An infected mosquito is capable
of transmitting the virus for the rest of its life. Mosquitoes’ spread is due to its tolerance to
temperature below freezing, hibernation, and ability to shelter in microhabitats.
At the present, there is no effective vaccine for the dengue fever. In late 2015 and early
2016, the first dengue vaccine was developed, but its efficacy depended on geographical set-
tings. The main method to control or prevent the transmission of dengue virus is through con-
trolling the environment including the mosquito’s habitats, applying appropriate insecticides,
using personal household protection, and monitoring and surveillance of vectors.
One of the difficulties to understand the dengue virus is how the virus can remain in human
population even through long periods of extremely low incidence. One hypothesis is that the
vertical transmission within the mosquito population allows the virus to persist during these
times. Vertical transmission of dengue virus by mosquitoes was discovered at the end of the
late 1970s. However, it is unclear how widespread it is in nature, and its importance in the
epidemiology of the disease. Vertical transmission of dengue virus has been demonstrated in
the lab for the several different mosquitoes. Numerous studies have provided clear evidence
of vertical transmission of dengue in wild Aedes aegypti and Aedes albopictus mosquitoes.
The laboratory experiments reviewed three ways of measuring vertical transmission. The first

2
is the vertical transmission rate (VTR) that is defined as the proportion of infected parents
that produce at least one infected offspring. The second is the filial infection rate (FIR) that
is defined as the proportion of infected progeny produced from infected parents, given that
vertical transmission has occurred. The third is the vertical infection rate (VIR), which is the
VTR multiplied by the FIR [1, 34].

Species vt (%) ft (%) v(%)

Ae aegypti 3.0 0.13 0.039

Ae albopictus 41.2 2.9 .2

Table 1.1: Average Laboratory vertical infection (v), vertical transmis-

sion (vt ), and filial infection (ft ) rates for Aedes mosquitoes. [16]

Ross first described malaria transmission mathematically and Macdonald updated and
extended Ross’s theory and applied it to the Global Malaria Eradication Programme (GMEP).
During the era of Macdonald, a quantitative theory, consisting of a set of linked concepts,
notation and metrics for understanding and measuring mosquito-borne pathogen transmission
and control were fully developed [61]. A number of factors that contribute to the rising of
vector-borne diseases include (1) the ability of the anthropoids to adapt to new habitats, (2)
development of insecticide and drug-resistant vectors, (3) global and rapid human movement
(by jet airplanes), (4) building widespread irrigation and water-impoundment, (5) civil unrest
and wars which lead to displacement of large masses of people who live for long periods of time
under poor conditions, (6) rapid urbanization which concentrates many host on small area, (7)
change in policies that took away resources for vector-control measures. In addition,the impact
of climate change and global warming is a topic of significant debate. The emergence and
reemergence of vector-borne diseases have promoted interest in their mathematical modeling.
[70]
Mathematical models have become important tools for analyzing the spread and control
of infectious diseases. The model formulation process clarifies assumptions, variables, and
parameters. The models provide conceptual results such as thresholds, basic reproduction

3
numbers, contact numbers, and recovered numbers. Also computer simulations are useful ex-
perimental tools for building and testing the theories, assessing quantitative conjectures, deter-
mining sensitivities to changes in parameter values, and estimating key parameters from data.
Mathematical models have been formulated for diseases such as measles, rubella, chickenpox,
whooping cough, diphtheria, smallpox, malaria, onchocerciasis, filariasis, rabies, gonorrhea,
herpes, syphilis, and HIV/AIDS [39].
The basic reproduction rate, R0 , is the number of secondary infections produced by one
primary infection in a totally susceptible population. The traditional threshold condition is
expressed in terms of relationship between S0 and ρ. S0 is the initial population of susceptible
individuals and ρ = βγ , where β is the transmission rate and γ is the recovery rate. If S0 > ρ
S0
the disease persists and if S0 < ρ the disease dies out. Because R0 = ρ
the condition S0 > ρ
is equivalent to the condition R0 > 1. Similarly, the condition S0 < ρ is equivalent to the
condition R0 < 1. If R0 > 1, then each infectious individual will pass the infection to more
than one susceptible individual. Therefore the disease can be maintained in the population. If
R0 < 1, then the disease will die out in the population because it is not able to reproduce itself
at a sufficient rate. This kind of information has proven that R0 is a useful concept to determine
effective control measures.
     
inf ectious contact time
R0 ∝ · · =τ ·c·d
contact time inf ectious
where τ is the transmissibility, i.e., probability of infection given contact between a suscepti-
ble and infected individual, c is the average rate of contact between susceptible and infected
individuals, and d is the duration of infectiousness.
A number of approaches have been used in the development of models analyzing the disea-
ses. These approaches include 1) compartment models, 2) statistical approaches, 3) geographic
approaches, and 4) economic models. Two classic epidemiology models are Epidemic and En-
demic models. Epidemic models are used to describe rapid outbreaks that occur in less than
one year. Endemic model are used for studying diseases over longer periods, during which
there is a renewal of susceptible individuals by births or recovery from temporary immunity.
The most common method in use is the compartment model. For example, there are SI, SIS,

4
SEI, SIR, SEIS, SIRS, SEIR, SEIRS, MSEIR, MSEIR, and MSEIRS models. The models consist
of a number of compartment based on the disease status of an individual.

• Passive immune (M): is composed by newborns that are temporarily passively immune
due to antibodies transferred by their mothers.

• Susceptible (S): is the class of individuals who are susceptible to infection. This can in-
clude the passively immune individuals once they lose their immunity or, more commonly,
any newborn infant whose mother has never been infected and therefore has not passed
on any immunity.

• Exposed or Latent (E): compartment refers to the individuals that despite being infected,
do not exhibit obvious signs of infection.

• Infected (I): in this class, the level of pathogen is sufficiently large within the host and
there is potential for transmitting the infection to other susceptible individuals.

• Recovered or Resistant (R): includes all individuals who have been infected but have
recovered.

The process of building a mathematical model begins with a series of assumptions about how
the disease process works and developing a simplified model to describe the process. The
choice of which compartments to include in a model depends on the characteristics of the
particular disease being studied and the purpose of the model. The exposed compartment is
sometimes neglected when the latent period is very short. Additionally, the compartment of
the recovered individuals cannot always be considered since there are diseases where the host
does not become resistant [62]. The general transfer diagram for the MSEIR is given as the
following.

5
 Figure 1.1: The general transfer diagram for the MSEIR model

The two classical epidemic and endemic SIR model provide an intuitive basis for under-
standing more complex epidemiology modeling results. For the classical SIR Epidemic model
we assume that 1) Constant (closed) population size, N , 2) Constant rates (e.g., transmission,
recovery rates), 3) No demography (i.e., births and deaths), and 4) Well-mixed population,
where any infected individual has a probability of contacting any susceptible individual that is
reasonably well approximated by the average. It is given by the initial value problem
dS −βIS
= , S(0) = S0 ≥ 0
dt N
dI βIS
= − γI, I(0) = I0 ≥ 0 (1.1)
dt N
dR
= γI, R(0) = R0 ≥ 0.
dt
where S(t), I(t), and R(t) are the number of susceptible, infectious, recovered, respectively
and N (t) = S(t) + I(t) + R(t). β is the effective contact rate, γ is the recovery rate. An
dI
epidemic occurs if the number of infected individuals increases, i.e., dt
> 0. Then we have
βIS S
N
− γI > 0 from the model (1.1). We assume that N
≈ 1 because at the outset of an
S
epidemic, nearly everyone is susceptible. Substituting N
= 1, we have the basic reproduction
β
number R0 = γ
> 1.
The classic Endemic SIR model is almost the same as the SIR epidemic model, except that
it has an inflow of newborns into the susceptible class at rate µN and deaths in the classes at

6
rates µS, µI, and µR. It is given by the initial value problem
dS −βIS
= µN − µS , S(0) = S0 ≥ 0
dt N
dI βIS
= − γI − µI, I(0) = I0 ≥ 0 (1.2)
dt N
dR
= γI − µR, R(0) = R0 ≥ 0.
dt
β
We find R0 = γ+µ in the similar way from the endemic SIR model.
dI
Since the two models are simple, we find the basic reproduction number, R0 from dt
. For a
complex model, we use the next generation method. The next generation method introduced by
Diekmann (1990) is a general method of deriving R0 in cases encompassing any situation in
which the populations is divided into discrete, disjoint classes. In the next generation method,
R0 is defined as the spectral radius (dominant eigenvalue) of the next generation operator
(matrix). Let us assume that there are n compartments of which m are infected. We define
the vector x̄ = {xi }ni=1 where xi denotes the number of individuals in the ith compartment.
Let Fi (x̄) be the rate of appearance of new infections into compartment i, and let Vi (x̄) =
Vi− (x̄) − Vi+ (x̄), where Vi+ (x̄) is the rate of transfer of individuals into compartment i by all
other means and V − (x̄) is the rate of transfer of individuals out of the ith compartment. The
difference Fi (x̄) − Vi (x̄) gives the rate of change of xi . We assume that Fi and Vi satisfy the
conditions outlined in Van den Driessche [66]. We can form the next generation matrix F V −1
where
 
∂Fi (x0 )
F = (1.3)
∂xj
and
 
∂Vi (x0 )
V = (1.4)
∂xj
where i, j = 1, . . . , m and x0 is the disease-free equilibrium. The entries of F V −1 give the
rate at which infected individuals in xj produce new infections in xi , times the average length
of time an individual spends in a single visit to compartment j. [38, 39, 41, 64].
We perform sensitivity analyses on a mathematical model to determine the relative im-
portance of model parameters to disease transmission and prevalence. With the sensitivity, we
can reduce human morbidity and control the disease. There are many methods available for

7
conducting sensitivity analysis such as differential analysis, response surface methodology, the
Fourier amplitude sensibility test (FAST) and other variance decomposition, fast probability
integration and sampling-based procedures. Nakul Chitnis [26] have evaluated the sensitivity
indices of the basic reproduction number and the point of endemic equilibrium to the parame-
ters in the model. We defines the normalized forward sensitivity index of a variable, u, that
depends differentiable on a parameter, p, as
∂u p
γpu = × . (1.5)
∂p u
These indices allow us to measure the relative change in a state variable when a parameter
changes [63].
Optimal control theory is a powerful mathematical tool to make decision involving a com-
plex system. For example, what percentage of the population should be vaccinated as times
evolves in a given epidemic model to minimize both the number of infected people and the cost
of implementing the vaccination strategy. Optimal control methods have been used to study the
dynamics of diseases including malaria, yellow fever, and dengue.
A typical optimal control problem requires a performance index or cost functional, J(x(·), u(·));
a set of state variable, x(·) ∈ X; and a set of control variable u(·) ∈ U . The main goal con-
sists in finding a piecewise continuous control u(t), t0 ≤ t ≤ tf , and the associated state
variable x(t), to minimize ( or maximize) the given objective functional. There are three
well known equivalent formulations to describe an optimal control problem, which are the
Lagrange, Mayer, and Bolza forms [24, 73].
The principal technique for an optimal control problem is to solve a set of necessary con-
ditions that an optimal control and corresponding state must satisfy. The necessary conditions
were developed by Pontryagin and his co-workers. Pontryagin introduced the idea of adjoint
functions to append the differential equation to the objective functional. Adjoint functions have
a similar purpose as Lagrange multipliers in multivariate calculus which append constraints
to the functions of several variables to be maximized or minimized. We need to find the appro-
priate conditions that the adjoint function should satisfy and derive a characterization of the
optimal control in terms of the optimal state and corresponding adjoint. We find the necessary

8
conditions from the Hamiltonian H, which is defined as follows.

H(t, x, u, λ) = f (t, x, u) + λg(t, x, u)

(1.6)
= integrand + adjoint ∗ RHS of the state equations.
We minimize (or maximize) H with respect to u at u∗ , and the conditions can be written in
terms of the Hamiltonian
∂H
= 0 at u∗ ⇒ fu + λgu = 0 (optimality conditin)
∂u
∂H
λ0 = − ⇒ λ0 = −(fx + λgx ) (adjoint equation) (1.7)
∂λ
λ(t1 ) = 0 (transversality condition)
There are two different ways to solve the optimal control problems in numerically 1) in-
direct method and 2) direct method. To apply an indirect method, it is necessary to expli-
citly get the adjoint equations, the control equations, and all of the transversality conditions if
they exist. To solve the state and adjoint differential equations, we use the backward-forward
sweep method and Rung-Kutta fourth order method . A direct method has been driven by
the industrial need to solve large-scale optimization problems. This method constructs a se-
quence of points x1 , x2 , . . . , x∗ such that the objective function F to be minimized satisfies
F (x1 ) > F (x2 ) > · · · > F (x∗ ) and reformulates the problem as a standard nonlinear op-
timization problem (NLP). There are many well-known software programs that can handle
it [47, 58, 62].
For the future research we consider periodicity and stochastic model. Periodicity and
other oscillatory have been observed in the incidence of many infections diseases, including
measles, mums, rubella, chickenpox, poliomyelitis, diphtheria, pertussis, and influenza. Many
researches show that models with periodic coefficients can explain the periodicity and other
oscillatory. For example, Contact rate vary seasonally for childhood diseases because of ope-
ning and closing of schools. Periodic changes in birth rate of populations are evidence in many
biological works. Vaccinations program is also a source of periodicity. A natural and im-
portant problem associated with periodic epidemic model is to define and compute their basic
reproduction numbers. Bacaër and Guernaoui [6] presented a general definition of the ba-
sic reproductions number in a periodic environment and show that it is a threshold parameter

9
for the local stability of the disease-free periodic solution and for the global dynamics under
certain circumstances [55, 68].
Stochastic models are characterized by randomness, and variable states are described by
probability distribution. If the environment is randomly varying and the population systems are
often subject to environment noise, then parameters involved in epidemic models are not ab-
solute constants, and they may fluctuate around some average values. If the initial population
size is small then a stochastic model is more appropriate, since the likelihood that the popula-
tion becomes extinct due to chance must be considered. Based on these factors, people began
to be concerned about stochastic epidemic models. There are different possible approaches
to including random effects in the model. In the future research, we will study the stochastic
model based on our deterministic model [20, 33, 51].
This dissertation includes the analysis and the optimal control of the deterministic verti-
cally transmitted vector-borne disease model. In chapter 2, we introduce a deterministic verti-
cally transmitted vector-borne disease model that uses the SEIR model for the host and the SEI
model for the vector. The basic reproduction number is derived using the next generation met-
hod and the local and global stability of the disease-free equilibrium point is discussed. Also,
the sensitivity for R0 is discussed. In chapter 3, we present the vertically transmitted vector-
borne epidemic model with two controls to derive an optimal prevention of the contact between
vector and host and an optimal treatment for host with the minimal implementation cost. We
introduce an optimal control problem under the given epidemic model and discuss the exis-
tence of the optimal controls and the optimality system to find the optimal controls. In chapter
4, we consider controlling the number of mosquitoes and prevention of human-mosquito itera-
tion. Similar to chapter 3, We introduce an optimal control problem under the given epidemic
model and discuss the existence of the optimal controls and the optimality system to find the
optimal controls. In each chapter, we provide the numerical results supporting the analytical
conclusions.

10
 Chapter 2

Analysis of Vertically transmitted vector-borne disease model

2.1 Vertically transmitted vector-borne disease Model

In this section, we formulate and analyze the vertically transmitted vector-borne disease model.
We use SEIR type of structure for host and SIR type of structure for vector. The host population
is grouped into four compartments: susceptible host (x1 ), exposed host (no symptom, x2 ),
infectious (x3 ), and treated (x4 ). The total population of host is N = x1 + x2 + x3 + x4 .
The vector population is grouped into three compartments: susceptible vector (y1 ), infectious
vector (y2 ). The total population of vector is P = y1 + y2 + y3 .
To formulate a model of a disease, we introduce parameters and assume the followings. We
consider that the vertical transmission. That is, a disease is transmitted vertically from mother
to child by blood transfusion, breast feeding, or complications during pregnancy and from
mosquito to mosquito’s eggs when they are infected. In the susceptible host, x1 , it is increased
by a result of new recruits and birth from susceptible, exposed, and treated hosts, and treated
hosts. It is decreased as a result of biting from infectious vectors and natural death. In the
φβy3 x1
exposed host, x2 , it is increased by a result of biting from infectious vectors (at a rate of N

) and birth from infected parents (at a rate of Λζ1 x3 ). it is decreased by a result of natural
death and becoming infectious host after the incubation period. In the infectious host, x3 , it is
increased by a result of infectious host from an exposed host after the incubation period. It is
decreased by a result of recovering, death induced by the disease, and natural death. In the
treated host, x4 , it is increased by a result of recovering from infectious host. It is decreased
by a result of becoming susceptible host after treated and natural death. In the susceptible

11
vector, y1 , it is increased by a result of new adult female and maturation from susceptible and
exposed vector. It is decreased by a result of biting the exposed and infectious host and natural
death. In the exposed vector, y2 , it is increased by a result of becoming a exposed vector from
φθ2 x2 φθ1 x3
biting the exposed (at a rate of N
) and infectious host (at a rate of N
) and birth from
infected parents. It is decreased by a result of becoming infectious after the incubation period
and natural death. In the infectious vector, y3 , it is increased by a result of becoming infectious
from exposed vector. It is decreased by natural death.

Parameter Description (rates are per day)

µ host death rate (when density is ignored)

ζ1 vertical disease transmission rate (in the host)

ζ2 vertical disease transmission rate (in the vector)

δ0 average number of new adult female mosquitoes

δ1 a factor for density dependent maturation of mosquitoes to adulthood

Λ per capita(person) birth rate of host

r host recovery rate

β the probability that the disease is transmitted from an infected vector to a host
per contact

ρ host recruitment rate (by birth, assumed susceptible)

φ the number of contacts between a host and a vector

α disease-induced host death rate

ψ fading rate of treatment to make hosts susceptible to the disease

γ density independent death rate of vectors

 incubation rate of the disease in a vector

d incubation rate of the disease in a host

θ1 transmission efficacy of the disease from infectious host to vector

θ2 transmission efficacy of the disease from exposed host to vector

Table 2.1: Parameters and description

12
 We consider the following transmission between human and vector [11], [12].

Figure 2.1: Compartment between human and vector where Sh = x1 ,

Eh = x2 , Ih = x3 , Th = x4 , Sv = y1 , Ev = y2 , and Iv = y3 .

Based on the above figure, we can establish a model with the parameters in Table 2.1 as
the following.
dx1 βφy3 x1
= ρ + Λ(x1 + x2 + x4 ) + Λ(1 − ζ1 )x3 + ψx4 − − µx1
dt N
dx2 βφy3 x1
= + Λζ1 x3 − dx2 − µx2
dt N
dx3
= dx2 − (r + α + µ)x3
dt
dx4
= rx3 − (ψ + µ)x4 (2.1)
dt
dy1 φθ1 x3 y1 φθ2 x2 y1
= δ0 + δ1 (y1 + y2 ) + δ1 (1 − ζ2 )y3 − − − γy1
dt N N
dy2 φθ1 x3 y1 φθ2 x2 y1
= + + δ1 ζ2 y3 − y2 − γy2
dt N N
dy3
= y2 − γy3
dt
with initial condition xi (0) ≥ 0, i = 1, 2, 3, 4 and yj (0) ≥ 0, j = 1, 2, 3 and t ∈ [0, T ].
To show that {x1 , x2 , x3 , x4 , y1 , y2 } are all bounded in a set, we find a positively invariant
set with the state system (2.1). We use the theorem from J.K. Hale. [37].

13
Theorem 2.1 (Differential Inequality). Let ω(t, u) be continuous scalar function on an open
connected set Ω ∈ R2 and such that the initial value problem for the scalar equation
du
= ω(t, u)
dt
has a uniques solution. If u(t) is a solution of the above equation on a ≤ t ≤ b and v(t) is a
solution of
dv
≤ ω(t, v(t))
dt
on a ≤ t ≤ b with v(a) ≤ u(a), then v(t) ≤ u(t) for a ≤ t ≤ b.

Adding the first four equations in (2.1), we have the differential equations of the total
population of host.
dN
= ρ + ΛN − αx3 − µN ≤ ρ + ΛN − µN. (2.2)
dt
dz
Let dt
= ρ + Λz − µz. Then it is one-dimensional differential equation with attracting set
ρ
[0, N ∗ ], where z ∗ = µ−Λ
is a positive equilibrium point. By the theorem 2.1, the equation
(2.2) implies that N (t) ≤ z(t) for N (0) ≤ z(0). Since z(t) is an autonomous equation and
N ∗ is the equilibrium solution, N (t) remains in [0, N ∗ ] for 0 ≤ N (0) ≤ N ∗ . Furthermore, if
N (0) > N ∗ , then N (t) approaches N ∗ .
Similarly, adding the last three equation, yields
dP
= δ0 + δ1 P − γP (2.3)
dt
δ0
with a positive equilibrium P ∗ = γ−δ1
and for 0 ≤ P (0) ≤ P ∗ , P (t) remains in [0, P ∗ ].
Let X = (x1 , x2 , x3 , x4 ) , Y = (y1 , y2 , y3 ), and define
4 3
X
∗
X ρ δ0
Ω = {(X, Y ) ∈ 4
R+ 3
×R+ , xi ∈ [0, N ], y1 ∈ [0, P ∗ ]} = [0, ]×[0, ] (2.4)
i=1 i=1
µ−Λ γ − δ1
ρ ρ
The set Ω is forward invariant and attractor. Also, for N (0) ≥ N ∗ = µ−Λ
, N (t) → N ∗ = µ−Λ
δ0 δ0
and for P (0) ≥ P ∗ = γ−δ1
, P (t) → P ∗ = γ−δ1
.

Theorem 2.2. Ω is positively invariant under system (2.1).

14
2.2 Disease free equilibrium point E0 and R0

In this section, we derive a disease free equilibrium point and the basic reproduction number,
R0 . The basic reproduction number is defined as the average number of secondary infections
produced when one infected individual is introduced into a host population [50, 53]. The basic
reproduction number R0 is often considered as the threshold quantity that determines when an
infection can invade and persist in a new host population [53].
To find the disease free equilibrium point of (2.1) we set x2 = x3 = x4 = 0 and y2 = y3 =
0. Let E0 = (x∗1 , 0, 0, 0, y1∗ , 0, 0) be the disease free equilibrium point. From the state system
(2.1), we obtain the following system.
dx1
= ρ + Λ(x1 ) − µx1
dt
dx2 dx3 dx4
= = =0
dt dt dt (2.5)
dy1
= δ0 + δ1 (y1 ) − γy1
dt
dy2 dy3
= =0
dt dt
ρ
Then the disease free equilibrium point is E0 = (x∗1 , 0, 0, 0, y1∗ , 0, 0) where x∗1 = µ−Λ
, y1∗ =
δ0
γ−δ1
,γ > δ1 , µ > Λ.
To find the epidemiology threshold R0 , we use the Next-Generation Approach [53,66]. The
key concept is that we need to average the expected number of new infectious over all possible
infected types. Let G be a next generation matrix in which the ijth element of G, gij , is the
expected number of secondary infectious of type i caused by a single infected individual of type
j. That is, each element of the matrix G is a reproduction number, but one where who infects
whom is accounted for. R0 is the average of all the elements of G. [29]

1. We consider equations in the state system (2.1) which correspond to the infected com-
partments which are related to x2 , x3 , x4 , y2 , y3 and let y2 = x5 and y3 = x6 .

2. We split the right-hand side in the infected compartments in the following way.
dxi
= Fi (x) − Vi (x), i = 2, 3, 4, 5, 6 (2.6)
dt
where

15
 • Fi (x) is the rate of appearance of new infection compartment i.

• Vi (x) = Vi (x)+ − Vi (x)− , Vi (x)+ is the rate of transfer of individuals into com-
partment i by all other means, and Vi− is the rate of transfer of individuals out of
compartment i.

Note that this decomposition may not be unique. Different decompositions may occurred
to different interpretations of each terms and may lead to a different basic reproduction
number. The decomposition should satisfy the following properties.

• if xi ≥ 0, then Fi (x), Vi (x)+ , Vi (x)− ≥ 0 for i = 2, 3, 4, 5, 6.

• xi = 0 if and only if Vi (x)− = 0

• Fi (x) = 0 and Vi (x) = 0 for i = 2, 3, 4 and x5 , x6 ≥ 0

6
P
• Vi (x) ≥ 0 for all xi ≥ 0.
i=2

3. Determine the matrices F and V with components.

   
∂Fi (x) ∂Vi (x)
F = and V = , for i, j = 2, 3, 4, 5, 6 (2.7)
∂xj x=E0 ∂xj x=E0

where E0 is DEF.

4. The Next-Generation matrix, K is defined as

K = F V −1 (2.8)

5. The basic reproduction number is defines as

R0 = ρ(F V −1 ) (2.9)

where ρ(A) denotes the spectral radius of A.

Definition 2.1. The spectral radius of a matrix A is defined as the maximum of the absolute
values of the eigenvalues of A

ρ(A) = sup{|λ| : λ ∈ σ(A)} (2.10)

where σ(A) denotes the set of eigenvalues of A.

16
 By the next-generation approach, we can find Fi (x) and Vi (x) and let Fe(x) = (Fi (x)) and
Ve (x) = (Vi (x)) for i = 2, 3, 4, 5, 6. Then
 
βφy3 x1 φθ1 x3 y1 φθ2 x2 y1
F (x) =
e + Λζ1 x3 , 0, 0, + + δ1 ζ2 y3 , 0
N N N (2.11)
Ve (x) = (k1 x2 , −dx2 + k3 x3 , −rx3 + k2 x4 , k4 y2 , −y2 + γy3 )
where k1 = d + µ, k2 = ψ + µ, k3 = r + α + µ, k4 =  + γ.
Next we calculate F and V which are Jacobian matrices of Fe and Ve respectively evaluated at
ρ
DFE, E0 , with N = µ−Λ
since x∗1 = N at DFE.
 
  0 Λζ1 0 0 βφ
∂F1 ∂F1 ∂F1 ∂F1 ∂F1
 
 dx2 dx3 dx4 dy2 dy3 

 0

 ∂F  0 0 0 0 
∂F2 
 2 ···
 
 dx dy3 
 
F =  .2 .. .. =
..  
 0 0 0 0 0 
..

 .. . . . .  



   θ2xφy 1 θ1 φy1
0 0 ζ2 δ1 
 
∂F5 ∂F5 1 x1
dx2
··· dy3
 
0 0 0 0 0
  (2.12)
k1 0 0 0 0
 
 
−d k3 0 0 0 
 
 
V =
 0 −r k2 0 0 

 
0 0 0 k4 0 
 
 
0 0 0 − γ
Then we see that V is a nonsingular M -matrix and we have V −1 .
 
1
0 0 0 0
 k1 
 d 1


 k1 k3 k3
0 0 0 

−1
 
V = dr r 1 (2.13)
 k1 k2 k3 k2 k3 k2 0 0 

 
1
 0 0 0 k4 0 
 
 
 1
0 0 0 γk4 γ

17
Definition 2.2. Let A be a n × n real Z-matrix. That is, A = (aij ) where aij ≤ 0 for all i 6= j,
1 ≤ i, j ≤ n. Then matrix A is also a non-singular M -matrix if it can be expressed in the form
A = sI − B, where B = (bij ) with bij ≥ 0 for all 1 ≤ i, j ≤ n, where s > ρ(B), and I is an
identity matrix , and A is a singular M -matrix if s = ρ(P ).

Definition 2.3. A matrix A is called an M -matrix if

• A has Z-pattern, that is, the off-diagonal elements of A are nonpositive.

• The inverse A exists and has nonnegative elements: A−1 ≥ 0.

Lemma 2.1. Let A be a non-singular M -matrix and suppose B and BA−1 are Z-matrices.
Then B is a non-singular M -matrix if and only if BA−1 is a non-singular M -matrix.

Since F is nonnegative matrix, K = F V −1 is also nonnegative matrix. This implies

that the next-generation matrix has its spectral radius as an eigenvalue and there are no other
eigenvalues with lager modulus by the Perron-Frobenius theorem. This largest eigenvalues
gives R0 . The basic reproduction number R0 is defined as the spectral radius of
 
dζ1 Λ ζ1 Λ βφ βφ
k1 k3 k3
0 k4 γ γ 

 

 0 0 0 0 0 

−1
 
FV =   0 0 0 0 0   (2.14)
 
 dy1 θ1 φ y1 θ2 φ y1 θ1 φ δ1 ζ2 δ1 ζ 2 
 k1 k3 x1 + k1 x1 k3 x1 0 k4 γ γ 


0 0 0 0 0
Thus, the spectral radius of F V −1 is
  s 2
−1 1 ζ1 Λd ζ δ
2 1  1 ζ1 Λd ζ2 δ1 
R0 = R0 (ζ1 , ζ2 ) = ρ(F V ) = + + ( − ) + Rh2 (2.15)
2 k1 k3 γk4 2 k1 k3 γk4
q
β(θ1 d+θ2 k3 )y1∗
where Rh = φ γk1 k3 k4 x∗
.
1

2.3 Stability of disease free equilibrium point

In this section, we show that the disease free equilibrium point, E0 , is locally asymptotically
stale and also globally asymptotically stable. First, we show the local stability of E0 .
To prove the local stability, we use the Hartman-Grobman theorem in Misha Guysinsky
[36].

18
Theorem 2.3 (Hartman-Grobman theorem). Let U ⊂ Rn be a neighborhood of 0, f : U → Rn
continuously differentiable with 0 as a hyperbolic fixed point. Then there is a homeomorphism
h of a neighborhood of 0 with h ◦ f = Df0 ◦ h near 0.

The theorem guarantees that the stability of the steady state of the original system is the
same as the stability of the trivial steady state of the linearized system.
Also, we use the linear stability analysis. It is stated as the following theorem.

Theorem 2.4. Given the differentia equations on Rn

x0 = f (x)

Let x0 be an equilibrium point of the above equation and A = Df (x0 ) be the Jacobian matrix
of f at the equilibrium point x0 . If all eigenvalues of A have strictly negative real part, then x0
is locally asymptotically stable.

We can see the local stability of E0 from the theorem in P. van den Driessche [66] and the
following theorems.

Theorem 2.5. The disease-free equilibrium point E0 , of system (2.1) is locally asymptotically
stable if R0 < 1 and unstable if R0 > 1, where R0 is defined by (2.15)

Next, we show the global stability of the disease free equilibrium point. To show that, we
use the following theorem in Castillo-Chávez, C [22]. First, the System (2.1) must be written
in the form
dx
= F (x, I)
dt (2.16)
dI
= G(x, I), G(x, 0) = 0
dt
where x ∈ Rn denotes (its component) the number of uninfected individuals and I ∈ Rn
denotes (its component) the number of infected individuals including exposed, infectious, etc.
The conditions (H1) and (H2) below must be to guarantee local asymptotic stability.

(H1) For dx
dt
= F (x, 0), x∗ is globally asymptotically stable.

(H2) G(x, I) = AI − G(x, b i) ≥ 0 for (x, I) ∈ Ω.

b I), G(x, 0) = 0 and G(x,

19
where A = DI G(x∗ , 0) is M-matrix (the off diagonal elements of A are nonnegative) and Ω
is (2.4) If the state system (2.1) satisfies the above two condition then the following theorem
holds.

Theorem 2.6. The fixed point E0 = (x∗ , 0) is globally asymptotically stable equilibrium of the
state system (2.1) provided that R0 < 1 and that assumptions (H1) and (H2) are satisfied.

Proof. First, we define new variables and break the state system (2.1) into subsystems. Let
I = (x2 , x3 , y2 , y3 ) and x = (x1 , x4 , y1 ). Then the state system (2.1) can be written as
dx
= F (x, I)
dt (2.17)
dI
= G(x, I)
dt
where
 βx1 y3 φ
F (x, I) = ρ − µx1 − + Λ(x1 + x2 + x4 ) + (1 − ζ1 )Λx3 + ψx4 , rx3 − x4 (µ + ψ),
N
θ1 x3 y1 φ θ2 x2 y1 φ T
δ0 − − + γ(−y1 ) + δ1 (y1 + y2 ) + δ1 (1 − ζ2 )y3
N N
(2.18)

 βx1 y3 φ
G(x, I) = − dx2 + − µx2 + ζ1 Λx3 , dx2 − x3 (α + µ + r),
N (2.19)
θ1 x3 y1 φ θ2 x2 y1 φ T
+ + γy2 − y2  + δ1 ζ2 y3 , y2  − γy3
N N
where T is transpose. To show (2.17) satisfies the condition (H1), consider the system dx
dt
=
F (x, 0)
dx1
= (x1 + x4 )Λ − x1 µ + ρ + x4 ψ
dt
dx4
= −x4 (µ + ψ) (2.20)
dt
dy1
= −y1 γ + δ0 + y1 δ1
dt
ρ δ0
Then, x∗ = (x∗1 , x∗4 , y1∗ ) = ( µ−Λ , 0, γ−δ 1
). We see x∗ is globally asymptotically stable under the
system (2.20). To see that, we solve the second equation in (2.20) and obtain

x4 (t) = e−t(µ+ψ) x4 (0) (2.21)

We have that x4 (t) → 0 as t → ∞. Similarly from the last equation, we obtain

δ0
y1 (t) = + e−t(γ−δ1 ) y1 (0) (2.22)
γ − δ1

20
 δ0
Since γ > δ1 , we have that y1 (t) → γ−δ1
as t → ∞. By solving the first equation using x4 (t),
we obtain
ρ e−t(µ+ψ) e−t(µ+ψ)
x1 (t) = − x4 (0) − Λ x4 (0) + e−t(µ−Λ) x1 (0) (2.23)
µ−Λ Λ+ψ Λ+ψ
ρ
Since µ > Λ, we have that x1 (t) → µ−Λ
as t → ∞.
We see that the convergence are independent of initial condition. Hence, the convergence
of solutions of (2.20) is global and x∗ is globally asymptotically stable.
Next, we show the system (2.17) satisfies the condition (H2). We see easily that G(x, 0) =
0. We find A = DI G(x∗ , 0) and G(x,
b I).
 
 −d − µ − x1Ny32βφ ζ1 Λ − x1 y3 βφ
N2
0 x1 βφ
N 
 
 d −r − α − µ 0 0 
DI G(x, I) = 
 

 y1 θ2 φ − x3 y1 θ1 φ − x2 y1 θ2 φ y1 θ1 φ
− x3 y1 θ1 φ
− x2 y1 θ2 φ
−γ −  δ1 ζ2 
 N N2 N2 N N2 N2 
 
0 0  −γ
(2.24)
Then, A = Di G(x∗ , 0)
 
 −d − µ ζ1Λ 0 βφ 
 
 d −r − α − µ 0 0 
A= (2.25)
 

 δ0θ2(µ−Λ)φ δ0θ1(µ−Λ)φ −γ −  δ1ζ2 
 (γ−δ1)ρ (γ−δ1)ρ 
 
0 0  −γ
and Let
 
βy3 φ(x2 +x3 +x4 )
 N 
 
 0 
G(x, I) =  (2.26)
b  

 φ(θ2 x2 +θ1 x3 )(δ0 x1(µ−Λ)+δ0 x2(µ−Λ)+δ0 (µ−Λ)x3 +δ0 (µ−Λ)x4 +(γ−δ1 )ρy1 ) 
 ρ(γ−δ1 )N 
 
0
b I) ≥ 0 for (x, I) ∈ Ω. We obtain G(x, I) = AI − G(x,
Since µ > Λ and γ > δ1 , G(x, b I).

Therefore, the disease free equilibrium point E0 is globally asymptotically stable.

21
2.4 Sensitivity analysis of R0

In this section we perform sensitivity analysis of the basic reproduction number already obtai-
ned, (2.15), to identify the parameters which are important in contributing variability in the
outcome of the basic reproduction number. There are several methods to perform the sensiti-
vity analysis. We use the fixed point estimation used in Samsuzzoha [63]. Sensitivity analysis
using the fixed point estimations has been applied to determine the relative importance of dif-
ferent parameters responsible for the disease transmission related to the basic reproductions
number. Sensitivity indices for the basic reproduction number change with the change in para-
meters values. The normalized forward sensitivity index ( in Nakul Chitnis [26] ) of a variable
to a parameter is the ratio of the relative change in the variable to the relative change in the
parameter. For example, let u be a variable that depends on p a parameter and δ > 0 be a
small perturbation corresponding to p. We have
u(p + δ) − u(p) ∂u
δu = u(p + δ) − u(p) = δ≈δ (2.27)
δ ∂p
We define the normalized forward sensitivity index , γpu as
δu δ p ∂u
γpu = / = (2.28)
u p u ∂p

Definition 2.4. The normalized forward sensitivity index of a variable , u, that depends diffe-
rentiably on a parameter, p, is defined as
∂u p
γpu = × (2.29)
∂p u

We evaluate the normalized forward sensitivity index for each parameters in R0 using the
definitiom and the values of parameters in the table 2.2.

22
 Parameters Values Parameters Values
1 1
γ [ 17 , 10 ] Λ (0, µ)
1 1
µ [ 70(365) , 45(365) ] θ1 [0, 1)

φ ≥1 θ2 (0, θ1 )
1 1
 [ 14 , 7] α (0, 0.001)

r [0, 17 ] ψ [0, 1)

ζ2 [0, 1) β (0, 1)
1 1
d [ 14 , 3] ζ1 [0, 1)

δ0 [700, 10000] δ1 [0, 1)

Table 2.2: Parameter estimation for dengue fever in [12, 13]. Rate is
per day.

For the simplicity, we set we set ζ1 = 0 throughout the sensitivity. The sensitivity index
γζR10 ≈ 1.91002 × 10−6 . It means that decreasing ( or increasing) ζ1 by 100% decreases ( or
increases) R0 by 1.91002 × 10−5 %. Since the formula for sensitivity index of parameters are
complicate, we see the sensitivity index formula for ζ2 as an example.
δ1 ζ2 
γζR20 = r   (2.30)
4βγδ0 φ2 (γ+)(µ−Λ)(dθ1 +θ2 (α+µ+r)) 2 2
 ρ(γ−δ1 )(d+µ)(α+µ+r)
+ δ1 ζ2 

Using the definition 2.4 and the formula for each paramters, we evaluate the normalized
forward sensitivity index of R0 .

23
 Parameter Sensitivity index

γ −2.19232

δ1 1.05291

φ 0.975482

µ 0.780124

δ0 0.487741

β 0.487741

ρ −0.487741

Λ −0.292644

θ1 0.25891

r −0.258659

θ2 0.22883

d −0.22864

 0.163923

ζ2 0.0245184

α −0.000181061

Table 2.3: Sensitivity index for parameters in R0 . The parameters are

orders from most sensitivity to [Link] are : θ1 = 0.0083,
1 0.375
θ2 = 0.00513, µ = 70(365) , Λ = 70(365) , ρ = 205, δ0 = 1050, β = 0.2,
1 1
r = 7 , δ1 = 0.0399, d = 10 , ψ = 0.00233, ζ2 = 0.0087, α = 0.0001,
φ = 2,  = 18 , γ = 17
1

The most sensitive parameter is the host recovery rate, γ. Other important parameters
include the factor for density dependent maturation of mosquitos to adulthood, δ1 , the number
of contacts between a host and a vector, φ, and the host death rate, µ. Since γγR0 = −2.19232,
decreasing ( or increasing) γ by 10% increases ( or decreases) R0 by 21.92%. Similarly, as
γδR10 = 1.05291, increasing ( or decreasing) δ1 by 10% increases ( or decreases) R0 by 10.5%.
The least important parameter is the disease-induced host date rate, α.

24
2.5 Numerical Simulations

In this section we perform simulation for the state system model (2.1). In this simulation, we
consider initial conditions, x1 (0) = 100, x2 (0) = 20, x3 (0) = 20, x4 (0) = 10, y1 (0) = 1000,
y2 (0) = 20, and y3 (0) = 30. These simulations are performed for different γ, δ1 , φ, R0 . R0 is
the basic reproduction number given by the equation (2.15). The values of parameters, Table
2.4 are considered based on the dengue virus, Table 2.2.
It is clear from Theorem 2.6 and 2.5 that the disease is endemic for R0 > 1. The numerical
simulation shows that the number of exposed and infectious host increase when the number of
the contact between a host and a vector, φ, increases.

Parameters Values Parameters Values

θ1 0.0083 θ2 0.00513
1 0.375
µ 70(365)
Λ 70(365)

ρ 205 δ0 1050
1
β 0.2 r 7
1
δ1 0.0399 d 10

ψ 0.00233 ζ1 0.00001
ζ2 0.0087 α 0.0001
1
φ 2  8
1
γ 17

Table 2.4: Parameter values for a dengue virus. In this case, R0 ≈

0.0838.

25
Figure 2.2: Solution for hosts, x1 , x2 , x3 , and x4 , with the values of the
parameters in Table 2.4.

26
 Figure 2.3: Solution for vectors y1 , y2 , and y3 , with the values of the
parameters in Table 2.4.

The disease free equilibrium is stable for R0 < 1 as it is given by Theorem 2.6. Since
R0 ≈ 0.0838 < 1 with the values of the parameters in Table 2.4, we can see in Figure 2.2 and
2.3 that the size of all exposed and infected groups in each population die out.

27
 Figure 2.4: Solution for Exposed, x2 , and Infectious,x3 , host for γ =
1 1 1 1
10
, γ = 13 , γ = 15 , γ = 17 with the values of the parameters in
Table 2.4. In each cases, R0 ≈ 0.0328, R0 = 0.0502, R0 = 0.065,
R0 = 0.0838, respectively.

The decreased death rate of vectors, γ, is a factor for the increased size of exposed and
infectious hosts in Figure 2.4. We see that the basic reproduction number R0 is increasing as γ
decreasing, as we studied in the section 2.4.

28
 Figure 2.5: Solution for Exposed, x2 , and Infectious,x3 , host for δ1 =
1
0.04, δ1 = 0.06, δ1 = 0.08 with γ = 10 and the values of the parameters
in Table 2.4. In each cases, R0 ≈ 0.0328, R0 ≈ 0.0405, R0 ≈ 0.0571,
respectively.

The decreased maturation of mosquitoes to adulthood, δ1 , is a factor for the decreased

size of exposed and infectious hosts in Figure 2.5. We see that the basic reproduction number
R0 is increasing as δ1 increasing, as we studied in the section 2.4.

29
 Figure 2.6: Solution for Exposed, x2 , and Infectious,x3 , host for φ = 2,
1
φ = 3, φ = 4 with δ1 = 0.05, γ = 17 , and the values of the parameters
in Table 2.4. In each cases, R0 ≈ 0.1223, R0 ≈ 0.1822, R0 ≈ 0.2421,
respectively.

The increased number of contacts between a host and a vector, φ, is a factor for the
increased size of exposed and infectious hosts in Figure 2.6. We see that the basic reproduction
number R0 is increasing as δ1 increasing, as we studied in the section 2.4.

30
 Figure 2.7: Solution for Exposed, x2 , and Infectious,x3 , host for φ =
1
17, φ = 18, φ = 19 with δ1 = 0.05, γ = 17 , and the values of the
parameters in Table 2.4. In each cases, R0 ≈ 1.0209, R0 ≈ 1.0808,
R0 ≈ 1.1407, respectively.

Since R0 > 1 for each φ, the disease free equilibrium point is unstable. In Figures 2.4,
2.5, 2.6, 2.7, we see that the number of exposed and infectious host are changing as γ, δ1 , and
φ is changing. However, it does not necessarily result in the spread of the disease among hosts.
This is shown by the phase portrait for φ = 17 in Figure 2.8.

31
 Figure 2.8: A phase plane portrait for the Exposed and the Infectious
host,x2 + x3 , against the Exposed and the Infectious vector, y2 + y3 ,
where φ = 17.

From the figure 2.8, we see that an increase in exposed and infected vectors does not
necessarily in the spread of the disease among hosts. Also, we see that the opposite case is not
necessary. So, we can not say about the change of the number of exposed and infected hosts
just looking at the change of the number of exposed and infected vectors.

32
 Chapter 3

Optimal control for Vertically transmitted vector-borne disease:treament and prevention of

human-mosquito interaction

3.1 A Model for Optimal Control of vertically transmitted vector-borne disease

In this section, we formulate an optimal control model for vertically transmitted vector-borne
disease to drive an optimal prevention of the contact between host and vector and an optimal
treatment for host with the minimal implementation cost. The control functions, u1 and u2 ,
represent time dependent efforts of prevention and treatment respectively on a time interval
[0, T ]. For the prevention, we can do the followings;

• Aware of the areas that are heavily infested and avoid

• Wear long-sleeved shirts and pants in neutral colors

• Inspect vehicles before entering

• Avoid bushes

• Use insect repellent

For the treatment, we can do the followings;

• Screen patient

• Administer drug intake and patients’ condition

33
We will introduce u1 and u2 into the model (2.1). The transmission model of the vertically
transmitted vector-borne disease with the prevention and the treatment controls is given by
dx1 βφy3 x1 (1 − u1 (t))
= ρ + Λ(x1 + x2 + x4 ) + Λ(1 − ζ1 )x3 + ψx4 − − µx1
dt N
dx2 βφy3 x1 (1 − u1 (t))
= + Λζ1 x3 − dx2 − µx2
dt N
dx3
= dx2 − (r + α + µ + r0 u2 (t))x3
dt
dx4
= (r + r0 u2 (t))x3 − (ψ + µ)x4
dt
dy1 φθ1 x3 y1 (1 − u1 (t)) φθ2 x2 y1 (1 − u1 (t))
= δ0 + δ1 (y1 + y2 ) + δ1 (1 − ζ2 )y3 − − − γy1
dt N N
dy2 φθ1 x3 y1 (1 − u1 (t)) φθ2 x2 y1 (1 − u1 (t))
= + + δ1 ζ2 y3 − y2 − γy2
dt N N
dy3
= y2 − γy3
dt
(3.1)
with initial condition xi (0) ≥ 0, i = 1, 2, 3, 4 and yj (0) ≥ 0, j = 1, 2, 3 and t ∈ [0, T ].
In the model (3.1), 1 − u1 (t) describes the failure rate of prevention efforts. The per capita
recovery rate is r0 u2 (t), where 0 ≤ r0 ≤ 1 is the proportion of effective treatment.
We discuss the boundedness of the host and vector population. From the model (3.1), we have
that by adding the first four equation,
dN
= ρ − αx3 + (Λ − µ)N ≤ ρ + (Λ − µ)N (3.2)
dt
ρ ρ
By the theorem (2.1), N ≤ µ−Λ
for the initial value N (0) ≤ µ−Λ
. Similarly, adding the
last three equations, we have
dP
= δ0 − (γ − δ1 )P (3.3)
dt
δ0 δ0
Thus, for the initial value P (0) ≤ γ−δ1
, we have P ≤ γ−δ1
.
Let X = (x1 , x2 , x3 , x4 ) and Y = (y1 , y2 , y3 ) and define a set
ρ δ0
Ω = {(X, Y )| R4+ × R3+ , 0 ≤ N ≤ ,0 ≤ P ≤ } (3.4)
µ−Λ γ − δ1

Theorem 3.1. Ω is positively invariant under system (3.1)

Proof. First, we show that the solutions with an initial values in Ω remains nonnegative for all
t ≥ 0. Let C1 = −βφ − µ, C2 = −d − µ, C3 = −γ − α − µ, and C4 = −ψ − µ. Then we

34
have the following inequalities
dxi
≥ Ci xi for i = 1, 2, 3, 4 (3.5)
dt
Similarly, let D1 = −φθ1 − φθ2 − γ, D2 = − − γ, and D3 = −γ. Then we have the following
inequalities
dyi
≥ Di yi for i = 1, 2, 3 (3.6)
dt
It implies that the solutions with an initial values in Ω remains nonnegative for all t ≥ 0. From
ρ δ0
. Since N = 4i=0 xi and P = 3i=0 yi , Ω is positively
P P
(3.2), (3.3) we have N ≤ µ−Λ and γ−δ 1

invariant under the system (3.1).

3.2 Optimal Control Problem

We consider an optimal control problem with the objective (cost) functional given by
Z T
J(u1 , u2 ) = (A1 x2 (t) + A2 x3 (t) + B1 u21 (t) + B2 u22 (t))dt + l(x1 (T ), x4 (T )) (3.7)
0

where A1 and A2 are positive weight constants of the susceptible and infectious group,
respectively and B1 and B2 are positive weight constants for prevention and treatment efforts,
respectively. We choose a quadratic for the cost on the controls for the technical reason and
that is similar in other literature, that is, B1 u21 is the cost of prevention, and B2 y22 is the cost
of the treatment effort. l(x1 (T ), x4 (T )) is the fitness of the susceptible and treated group at
the end of the process as a result of the prevention and the treatment efforts and we want to
maximize while the cost function J is minimized. We seek an optimal control pair (u∗1 , u∗2 ) such
that

J(u∗1 , u∗2 ) = min{J(u1 , u2 )|(u1 , u2 ) ∈ Γ} (3.8)

where the control set is

Γ = {(u1 , u2 )|ui (t) is piecewise continuous on[0, T ], ai ≤ ui (t) ≤ bi , i = 1, 2} (3.9)

where ai and bi , i = 1, 2 are constants in [0, 1]. We discuss the existence of the optimal control
and then the optimal system.

35
3.3 Existence of an Optimal Control

In this section, we show that the optimal control exists by using a result from Fleming and
Rishel [33] and Carathodory’s existence theorem [49].

Theorem 3.2 (Carathodory’s existence theorem). Consider the differential equation

y 0 (t) = f (t, y(t)), y(t0 ) = y0 (3.10)

with f defined in the rectangular domain R = {(t, y)||t − t0 | ≤ a, |y − y0 | ≤ b}. If the function
f satisfies the following three conditions

• f (t, y) is continuous in y for each fixed t.

• f (t, y) is measurable in t for each fixed y.

• there is a Lebesgue-integrable function m(t), |t − t0 | ≤ a, such that |f (t, y)| ≤ m(t) for
all (t, y) ∈ R.

then, the differential equation has a solution in the extended sense in a neighborhood of the
initial condition.

Theorem 3.3. Consider the objective functional J(u1 , u2 ) given by (3.7) with (u1 , u2 ) ∈ Γ sub-
jected to the system (3.1). There exists (u∗1 , u∗2 ) ∈ Γ such that J(u∗1 , u∗2 ) = min{J(u1 , u2 )|(u1 , u2 ) ∈
Γ}.

Proof. By a result from Fleming and Rishel [33], if the following conditions are satisfied, then
there exist (u∗1 , u∗2 ) ∈ Γ.

1. The set of controls and corresponding state variables are non-empty.

2. Γ is convex and closed

3. The right hand side of (3.1) is bounded by a linear function in the state and control.

4. The integrand of the equation (3.7) is convex on Γ and is bounded below by c1 (|u1 |2 +
β
|u2 |2 ) 2 − c2 where c1 > 0, c2 > 0, and β > 1.

36
 5. The function l is continuous.

Carathodory’s existence theorem for the state system (3.1) with bounded coefficients gives the
first condition. The control set Γ is convex and closed by definition. The right hand side of the
system (3.1) satisfies the third condition as the state solutions are bounded from Theorem 3.1.
The integrand in the objective functional 3.7 is convex. Let L(t, x1 , x2 , u1 , u2 ) = A1 x2 (t) +
A2 x3 (t) + B1 u21 + B2 u22 . We show that

L(t, x1 , x2 , λu1 + (1 − λ)u01 , λu2 + (1 − λ)u02 )

(3.11)
≤ λL(t, x1 , x2 , u1 , u2 ) + (1 − λ)L(t, x1 , x2 , u01 , u02 )
for λ ∈ [0, 1] and (u1 , u2 ), (u01 , u02 ) ∈ Γ. By the definition, we have that

L(t, x1 , x2 , λu1 + (1 − λ)u01 , λu2 + (1 − λ)u02 )

(3.12)
= A1 x2 (t) + A2 x3 (t) + B1 (λu1 + (1 − λ)u01 )2 + B2 (λu2 + (1 − λ)u02 )2
and
λL(t, x1 , x2 , u1 , u2 ) + (1 − λ)L(t, x1 , x2 , u01 , u02 )
(3.13)
= A1 x2 (t) + A2 x3 (t) + λB1 u21 + λB2 u22 + B1 (1 − λ)u02
1 + B2 (1 − λ)u02
2

Then we get that

A1 x2 (t) + A2 x3 (t) + B1 (λu1 + (1 − λ)u01 )2 + B2 (λu2 + (1 − λ)u02 )2

(3.14)
≤ A1 x2 (t) + A2 x3 (t) + λB1 u21 + λB2 u22 + B1 (1 − λ)u02 02
1 + B2 (1 − λ)u2

It means that we need to show that

B1 (λu1 + (1 − λ)u01 )2 + B2 (λu2 + (1 − λ)u02 )2

≤ λB1 u21 + λB2 u22 + B1 (1 − λ)u02 02

1 + B2 (1 − λ)u2
(3.15)
0 0
= B1 (λu21 (1 − λ)u12 ) + B2 (λu22 + (1 − λ)u22 )

37
We show that (λu1 + (1 − λ)u01 )2 ≤ λu21 + (1 − λ)u02
1.

(λu1 + (1 − λ)u01 )2 = λ2 u21 + 2λ(1 − λ)u1 u01 + (1 − λ)2 u02

1

= λ(λ − 1 + 1)u21 + 2λ(1 − λ)u1 u1 + (1 − λ)(1 − λ)u02

1

= λu21 + λ(λ − 1)u21 + 2λ(1 − λ)u1 u01 + (1 − λ)u02 02

1 + λ(λ − 1)u1

= λu21 + (1 − λ)u02 2 0 02
1 + λ(λ − 1)(u1 − 2u1 u1 + u1 )

= λu21 + (1 − λ)u02 0 2
1 + λ(λ − 1)(u1 − u1 )

≤ λu21 + (1 − λ)u02
1
(3.16)
Similarly, we can show that (λu2 + (1 − λ)u02 )2 ≤ λu22 + (1 − λ)u02
2 . Thus, the integrand in

the objective functional , A1 x2 (t) + A2 x3 (t) + B1 u21 + B2 u22 , is convex on Γ. Since the state
variables are bounded, there are c1 > 0 , c2 > 0 and β > 1 satisfying
β
A1 x2 (t) + A2 x3 (t) + B1 u21 + B2 u22 ≥ c1 (|u1 |2 + |u2 |2 ) 2 − c2

We want to maximize the function l(x1 (T ), x4 (T )) while the cost functional J is minimized.
If we define the function l as follows,

l(x1 (T ), x4 (T )) = −Q1 x1 (T ) − Q2 x4 (T ), Q1 ≥ 0, Q2 ≥ 0 (3.17)

Then the function l is clearly continuous. Finally there exists an optimal control pair (u∗1 , u∗2 )
that minimizes the objective functional J(u1 , u2 ).

3.4 Optimality System

We present the optimality system using a result from Lewis and Syrmos [48] and Pontryagin’s
Maximum Principle. From the theorems in Lenhart and Workman [47] and Clarke [28], the
optimality system can be used to compute candidates for the optimal control pair.

Theorem 3.4. If u∗ (t) and x∗ (t) are optimal for

Z t1
minu f (t, x(t), u(t))dt
t0

subjected to x0 (t) = g(t, x(t), u(t)) (3.18)

x(t0 ) = x(0) and x(t1 ) free.

38
,then there exists a piecewise differentiable adjoint variable λ(t) such that

H(t, x∗ (t), u∗ (t), λ(t)) ≤ H(t, x∗ (t), u(t), λ(t)) (3.19)

for all controls u at each time t, where Hamiltonian H is

H = f (t, x(t), u(t)) + λ(t)g(t, x(t), u(t)) (3.20)

and
∂H(t, x∗ (t), u∗ (t), λ(t))
λ0 (t) = −
∂x (3.21)
λ(t1 ) = 0

Theorem 3.5. Suppose that f (t, x, u) and g(t, x, u) are noth continuously differentiable functi-
ons in their three arguments and concave in u. Suppose u∗ is an optimal control for problem
(3.18), with associated state x∗ , and λ a piecewise differentiable function with λ(t) ≥ 0 for all
t. Suppose for all t0 ≤ t ≤ t1

0 = Hu (t, x∗ t(t), u∗ (t), λ(t)). (3.22)

Then for all controls u and each t + 0 ≤ t ≤ t1 , we have

H(t, x∗ (t), u∗ (t), λ(t)) ≤ H(t, x∗ (t), u(t), λ(t)) (3.23)

We have introduced an adjoint variable λ, which is similar to a Lagrange multiplier. It

attaches the differential equations information onto the minimization of the objective functional
(3.7). To apply the theory, we follow an outline.

1. Form the Hamiltonian for the problem.

2. Write the adjoint differential equation, termianl boundary condition, and the optimiality
condition. Now there are three unknowns, u∗ , x∗ , and λ.

3. Try to eliminate u∗ by using the optimality equation HU = 0, i.e., solve for u∗ in terms of
x∗ and λ.

4. Solve the two differential equations for x∗ and λ with two boundary conditions, substitu-
ting u∗ in the differential equations with the expression for the optimal control from the
previous step.

39
 5. After finding the optimal state and adjoint, solve for the optimal control.

We define a Lagrangian which is the Hamiltonian augmented with penalty terms for the
control constraints. Let Z = (X, Y ) ∈ Ω and U = (u1 , u2 ) ∈ Γ, where X = (x1 , x2 , x3 , x4 )
and Y = (y1 , y2 , y3 ), where Ω is defined by (3.4) and Γ is defined by (3.9).

H(Z, U, Π) =A1 x2 (t) + A2 x3 (t) + B1 u21 + B2 u22

 
βφy3 x1 (1 − u1 (t))
+ λ1 ρ + Λ(x1 + x2 + x4 ) + Λ(1 − ζ1 )x3 + ψx4 − − µx1
N
 
βφy3 x1 (1 − u1 (t))
+ λ2 + Λζ1 x3 − dx2 − µx2
N
+ λ3 (dx2 − (r + α + µ + r0 u2 (t))x3 )

+ λ4 ((r + r0 u2 (t))x3 − (ψ + µ)x4 )

 
φθ1 x3 y1 (1 − u1 (t)) φθ2 x2 y1 (1 − u1 (t))
+ λ5 δ0 + δ1 (y1 + y2 ) + δ1 (1 − ζ2 )y3 − − − γy1
N N
 
φθ1 x3 y1 (1 − u1 (t)) φθ2 x2 y1 (1 − u1 (t))
+ λ6 + + δ1 ζ2 y3 − y2 − γy2
N N
+ λ7 (y2 − γy3 )

− ω11 (u1 − a1 ) − ω12 (b1 − u1 ) − ω21 (u2 − a2 ) − ω22 (b2 − u2 )

(3.24)
where the adjoint variable Π = (λ1 , λ2 , λ3 , λ4 , λ5 , λ6 , λ7 ) and the penalty multipliers ωij (t) ≥
0 satisfying

ω11 (u1 (t) − a1 ) = ω12 (t)(b1 − u1 (t)) = 0 at optimal control u∗1

and

ω21 (u2 (t) − a2 ) = ω22 (t)(b2 − u2 (t)) = 0 at optimal control u∗2

We can check the concavity condition of H that minimizes the objective functional;
∂ 2H
= 2Bi > 0 at u∗i for i = 1, 2. (3.25)
∂u2i

40
 Theorem 3.6. Given an optimal control pair, (u∗1 , u∗2 ),and solutions, x1 , x2 , x3 , x4 , y1 , y2 and
y3 of the corresponding state system (3.1), there exist adjoint variables Π satisfying
    
β(1 − u 1 )x 1 y3 φ β(1 − u 1 )y 3 φ β(1 − u1 )y 3 φ β(1 − u 1 )x 1 y3 φ
λ˙1 = − λ1 Λ − µ + − + λ2 −
N2 N N N2
   
θ2 (1 − u1 )x2 y1 φ θ1 (1 − u1 )x3 y1 φ θ2 (1 − u1 )x2 y1 φ θ1 (1 − u1 )x3 y1 φ
+ λ5 + + λ6 − −
N2 N2 N2 N2
    
β(1 − u1 )x1 y3 φ β(1 − u1 )x1 y3 φ
λ˙2 = − A1 + dλ3 + λ2 −d − µ − 2
+ λ1 Λ +
N N2
 
θ2 (1 − u1 )x2 y1 φ θ1 (1 − u1 )x3 y1 φ θ2 (1 − u1 )y1 φ
+ λ5 + −
N2 N2 N
 
θ2 (1 − u1 )x2 y1 φ θ1 (1 − u1 )x3 y1 φ θ2 (1 − u1 )y1 φ
+ λ6 − − +
N2 N2 N
  
β(1 − u1 )x 1 y 3 φ
λ˙3 = − A2 + λ3 (−α − µ − r0 u2 − r) + λ4 (r0 u2 + r) + λ1 (1 − ζ1 )Λ +
N2
   
β(1 − u1 )x1 y3 φ θ2 (1 − u1 )x2 y1 φ θ1 (1 − u1 )x3 y1 φ θ1 (1 − u1 )y1 φ
+ λ2 ζ1 Λ − + λ5 + −
N2 N2 N2 N
 
θ2 (1 − u1 )x2 y1 φ θ1 (1 − u1 )x3 y1 φ θ1 (1 − u1 )y1 φ
+ λ6 − − +
N2 N2 N
  
β(1 − u1 )x1 y3 φ βλ2 (1 − u1 )x1 y3 φ
λ˙4 = − λ4 (−µ − ψ) + λ1 Λ + 2
+ψ −
N N2
   
θ2 (1 − u1 )x2 y1 φ θ1 (1 − u1 )x3 y1 φ θ2 (1 − u1 )x2 y1 φ θ1 (1 − u1 )x3 y1 φ
+ λ5 + + λ6 − −
N2 N2 N2 N2
    
θ2 (1 − u 1 )x 2 φ θ1 (1 − u 1 )x 3 φ θ 2 (1 − u 1 )x 2 φ θ1 (1 − u1 )x 3 φ
λ˙5 = − λ5 −γ + δ1 − − + λ6 +
N N N N
 
λ˙6 = − λ6 (−γ − ) + δ1 λ5 + λ7 
 
βλ 1 (1 − u 1 )x 1 φ βλ 2 (1 − u 1 )x 1 φ
λ˙7 = − − γλ7 + δ1 (1 − ζ2 ) λ5 + δ1 ζ2 λ6 − +
N N
(3.26)
with the terminal conditions,
∂l ∂l
λ1 (T ) = , λ4 (T ) = , λi (T ) = 0, for i = 2, 3, 5, 6, 7. (3.27)
∂x1 T ∂x4 T

Furthermore, u∗1 and u∗2 are represented by

  
∗ βλ1 (−x1 )y3 φ + βλ2 x1 y3 φ + (λ6 − λ5 ) y1 φ (θ2 x2 + θ1 x3 )
u1 = max a1 , min b1 ,
2B1 N
  
(λ3 − λ4 ) r0 x3
u∗2 = max a2 , min b2 ,
2B2
(3.28)

41
Proof. We differentiate the Lagrangian H with respect to states, Z = (x1 , x2 , x3 , x4 , y1 , y2 , y3 ).
Then the adjoint system can be written as
∂L ∂L ∂L
λ̇1 = − , λ̇2 = − , λ̇3 = − ,
∂x1 ∂x2 ∂x3
∂L ∂L ∂L ∂L
λ̇3 = − , λ̇5 = − , λ̇6 = − , λ̇7 = −
∂x4 ∂y1 ∂y2 ∂y3
The terminal condition of the adjoint equations can be given by
∂l
− Π = 0, at t = T.
∂Z
To obtain the optimality conditions, we differentiate the Lagrangian H with respect to U =
(u1 , u2 ) and set it equal to zero.
 
∂L βλ1 x1 y3 φ βλ2 x1 y3 φ θ2 x2 y1 φ θ1 x3 y1 φ
= 2B1 u1 + − + λ5 +
∂u1 N N N N
 
θ2 x2 y1 φ θ1 x3 y1 φ
+ λ6 − − − ω11 + ω12 = 0 (3.29)
N N
∂L
= 2B2 u2 − λ3 r0 x3 + λ4 r0 x3 − ω21 + ω22 = 0
∂u2
Solving for the optimal control, we obtain
N (ω11 − ω12 ) + βλ1 (−x1 )y3 φ + βλ2 x1 y3 φ + (λ6 − λ5 ) y1 φ (θ2 x2 + θ1 x3 )
u∗1 =
2B1 N
(3.30)
∗ (λ3 − λ4 ) r0 x3 + ω21 − ω22
u2 =
2B2
We consider the following three cases to have an explicit expression for the optimal control.
For the optimal control u∗1 ,

1. On the set {t|a1 < u∗1 (t) < b1 }, we have ω11 (t) = ω12 = 0. Hence the optimal control
is
βλ1 (−x1 )y3 φ + βλ2 x1 y3 φ + (λ6 − λ5 ) y1 φ (θ2 x2 + θ1 x3 )
u∗1 =
2B1 N

2. On the set {t|u∗1 (t) = b1 }, we have ω11 (t) = 0. Hence

N (−ω12 ) + βλ1 (−x1 )y3 φ + βλ2 x1 y3 φ + (λ6 − λ5 ) y1 φ (θ2 x2 + θ1 x3 )
b1 = u∗1 =
2B1 N
Since ω12 (t) ≥ 0, we have that
βλ1 (−x1 )y3 φ + βλ2 x1 y3 φ + (λ6 − λ5 ) y1 φ (θ2 x2 + θ1 x3 )
≥ b1
2B1 N

42
 3. On the set {t|u∗1 (t) = a1 }, we have ω12 (t) = 0. Hence
N (ω11 ) + βλ1 (−x1 )y3 φ + βλ2 x1 y3 φ + (λ6 − λ5 ) y1 φ (θ2 x2 + θ1 x3 )
a1 = u∗1 =
2B1 N
Since ω11 (t) ≥ 0, we have that
βλ1 (−x1 )y3 φ + βλ2 x1 y3 φ + (λ6 − λ5 ) y1 φ (θ2 x2 + θ1 x3 )
≤ a1
2B1 N

Combining these three cases, the optimal control u∗1 is characterized as

  
∗ βλ1 (−x1 )y3 φ + βλ2 x1 y3 φ + (λ6 − λ5 ) y1 φ (θ2 x2 + θ1 x3 )
u1 = max a1 , min b1 ,
2B1 N
For the optimal control u∗2 ,

1. On the set {t|a2 < u∗2 (t) < b2 }, we have ω21 (t) = ω22 = 0. Hence the optimal control
is
(λ3 − λ4 ) r0 x3
u∗2 =
2B2

2. On the set {t|u∗2 (t) = b2 }, we have ω21 (t) = 0. Hence

(λ3 − λ4 ) r0 x3 − ω22
b2 = u∗2 =
2B2
Since ω22 (t) ≥ 0, we have that
(λ3 − λ4 ) r0 x3
≥ b2
2B2

3. On the set {t|u∗2 (t) = a2 }, we have ω22 (t) = 0. Hence

(λ3 − λ4 ) r0 x3 + ω21
u∗2 =
2B2
Since ω21 (t) ≥ 0, we have that
(λ3 − λ4 ) r0 x3
≥ a2
2B2

Combining these three cases, the optimal control u∗2 is characterized as

  
∗ (λ3 − λ4 ) r0 x3
u2 = max a2 , min b2 ,
2B2

43
3.5 The forward-backward sweep method Algorithm

In this section , we introduce a numerical algorithm to find an optimal control, u∗ . We break

the time interval [t0 , t1 ] into pieces with specific points of interest t0 = b1 , b2 , . . . , bN , bN +1 =
t1 . These points will usually be equally spaced. The approximation will be a vector →
−
u =
(u1 , u2 , . . . , uN , uN +1 ) where ui = (u1 (bi ), u2 (bi ))T . There are various methods to solve op-
timal control problem. For example, Rodrigues and Monteiro [62], Wang [69], or Betts [10].
We use the forward-backwardsweep method introduced in Anita and Capasso [3] or Lenhart
and Workman [47]. We consider a control problem.
Z t1
minu J(u) = f (t, x(t), u(t))dt (3.31)
t0

subjected to

x0 (t) = g(t, u(t), x(t)), t ∈ (t0 , t1 )


(3.32)

x(t0 ) = x0


and

∗ ∗
 λ0 (t) = − ∂H(t,x (t),u (t),λ(t))


∂x
(3.33)

 λ(t1 ) = 0


where G(t, u, x) is the integrand of the cost functional. Then the forward-backward sweep
method is

S0: Choose u(0) ∈ Γ. Set k = 0

S1: Compute x(k) the solution to (3.32) corresponding to u = u(k)

S2: Compute λ(k) the solution to (3.33) corresponding to u = u(k)

S3: Compute u(k+1) the solution to the equation

Gu (t, u(t), x(k) (t)) + fu∗ (t, u(t), x(k) (t))λ(k) = 0

S4: (The stopping criterion)

If min(δku(k+1) k1 −ku(k+1) −u(k) k1 , δkx(k+1) k1 −kx(k+1) −x(k) k1 , δkλ(k+1) k1 −kλ(k+1) −
λ(k) k1 ) > 0

44
 then Stop
else k = k + 1 go to S0.

Here, k · k1 refers to the l1 norm for vectors, i.e.,

N
X +1
ku − olduk = |ui − oldui |. (3.34)
i=1
Many types of convergence tests exist for S5. We use the stopping criterion explained in Lenhart
and Workman [47]. It requires the relative errors to be small,
ku − olduk1
≤δ
kuk1
kx − oldxk1
≤δ (3.35)
kxk1
kλ − oldλk1
≤δ
kλk1
where δ is the accepted tolerance. To allow the zero for the controls, multiply both sides by
kuk1 in (3.35) to remove it from the denominator. Then we have

δkuk1 − ku − olduk ≥ 0

or
N
X +1 N
X +1
δ |ui | − |ui − oldui | ≥ 0. (3.36)
i=1 i=1
In the same reason, we have

δkxk1 − kx − oldxk ≥ 0

δkλk1 − kλ − oldλk ≥ 0

This method has two restrictions as explained in Lenhart and Workman [47], 1) the Lipschitz
constants for the state, adjoint, and control is small enough and 2) the time interval is small.
Because of these restrictions, we choose the parameters and t1 very carefully. The convergence
and stability of the forward-backward sweep algorithm can be found in Mcasey [54].

3.6 The Runge-Kutta method in 3-dimension

In this section, we introduce a numerical method required to solve in S1 and S2 in the previous
section. To solve the state (3.1) , the Runge-Kutta method is applied and we consider a modified
Runge-Kutta to solve the adjoint (3.26). We solve the state (2.1) forward in time and the adjoint

45
(3.26) backward in time. We consider 4th order Runge-Kutter methods. Let x0 (t) = f (t, x(t)).
Then the 4th order Runge-Kutta method is
h
xn+1 = xn + (k1 + 2k2 + 2k3 + k4 ) for n = 0, 1, 2, 3, . . . .
6
where
k1 = f (tn , xn )
h k1
k2 = f (tn + , xn + h)
2 2 (3.37)
h k2
k3 = f (tn + , xn + h)
2 2
k4 = f (tn + h, xn + hk3 )
k1
To find k2 in (3.37), xn is replaced with xn + 2
h and tn is replaced with tn + h2 . So, to
calculate a control, u, we should consider un + h2 k2 . However, there is no explicit dependence
on t in the differential equation for u. So this value is not assigned by our vector. There are
many ways to approximate this value. For example, an interpolating polynomial or spline of
u could be generated. In most of the literature, it usually suffices to approximate it as the
following

un (1 − cn ) + un+1 cn (3.38)

where n is the current iteration and 0 < cn < 1. This is weighted average, where the weight
shifts each iteration towards the current iteration. In our numerical experiment, we use the
average
un + un+1
. (3.39)
2
We consider the 4th order Runge-Kutta for 3 inputs, so we can solve the states forward in
time,

k1 = f (tn , xn , un )
h h 1
k2 = f (tn + , xn + k1 , (un + un+1 ))
2 2 2 (3.40)
h h 1
k3 = f (tn + , xn + k2 , (un + un+1 ))
2 2 2
k4 = f (tn + h, xn + hk3 , un+1 )

46
To solve the adjoints backward in time,

l =N +2−n

k1 = f (tl , λl , xl , ul )
h 1 1
k2 = f (tl − k1 , (xl + xl−1 ), (ul + ul−1 ))
2 2 2
(3.41)
h 1 1
k3 = f (tl − k2 , (xl + xl−1 ), (ul + ul−1 ))
2 2 2
k4 = f (tl − h, λl − hk3 , xl−1 , ul−1 )
h
λl−1 = λl − (k1 + 2k2 + 3k3 + k4 )
6
where h is the step size between time, t, N is the total number of time steps, n = 1, 2, 3 . . . , N ,
and l = N + 1, N, . . . , 2.
The error for the 4th order Runge-Kutta method is O(h4 ). The stability and accuracy of
the 4th order Runge-Kutta method is found in Butcher [17, 18].

3.7 Numerical Results

In this section, we perform a simulation for the state system (3.1), the adjoint system (3.26),
and the optimal control (3.30). The model considered in the experiment is tested with data
taken from the dengue virus. The optimality system is a two-point boundary problem because
of the initial condition Z(0) of the state system (3.1) and the terminal condition Π(0) (3.26).
First, we make an initial guess for the control functions. Second, we solve the initial valued
state system forward in time. Then, using the same guess for the control functions, we solve the
adjoint system with the terminal conditions backeard in time. The controls are updated in each
iteration using the optimality conditions (3.30). To focus on the controls, we choose weight
constant values A1 = A2 = 1, B1 = B2 = 50, and Q1 = Q2 = 0.1 in the objective functional
(3.7) and the Hamitonian (3.24).
We consider the initial conditions x1 (0) = 100, x2 (0) = 20, x3 (0) = 20, x4 (0) = 10, y1 (0) =
1000, y2 (0) = 20 and y3 (0) = 30. For the boundary of prevention and treatment efficiency we
choose a1 = a2 = 0 and b1 = b2 = 1.

47
 Parameter Value Parameter Value
1
 12
δ0 1050
1 1
γ 15
r 7

β 0.2 r0 0.04

ζ2 0.67 δ1 0.0399
1
θ1 0.0082 d 10

θ2 0.0289 ψ 0.0014
1.01
µ 70∗365
ζ1 0
0.379
Λ 70∗365
α 0.0238

ρ 205 φ 3

Table 3.1: Parameter values are estimated based on a dengue virus.

Figure 3.1: Comparison between Controlled(left) and Uncontrol-

led(right) for susceptible host

48
Figure 3.2: Comparison between Controlled(left) and Uncontrol-
led(right) for exposed host

Figure 3.3: Comparison between Controlled(left) and Uncontrol-

led(right) for infected host

49
Figure 3.4: Comparison between Controlled(left) and Uncontrol-
led(right) for treated host

Figure 3.5: Comparison between Controlled(left) and Uncontrol-

led(right) for susceptible vector

50
Figure 3.6: Comparison between Controlled(left) and Uncontrol-
led(right) for exposed vector

Figure 3.7: Comparison between Controlled(left) and Uncontrol-

led(right) for infected vector

51
 Figure 3.8: Optimal Control u1 (left) and u2 (right)

From 3.3 and 3.4, we see that the number of exposed and infected host are reduced in
very short time compare to the result without controls. Also, From 3.6 and 3.7, we see that the
number of exposed and infected vector are reduced in vary short time compare to the result
without controls. We see that from the figures 3.1 - 3.7, the result with the optimal controls, u1
and u2 are better than the result without the optimal control.
From 3.8, we see that the effectiveness of the prevention is higher than the effectiveness of
the treatment. So, we see the numerical result with only with the prevention u1 .

52
 Chapter 4

Summary

In this dissertation we studied a vertically transmitted deterministic vector-borne disease model

and the optimal control of the deterministic model. The deterministic model is a compartment
model. We consider SEIR model for the host and SEI model for the host. We divide the total
host population into susceptible, exposed, infectious, and treated groups. Similarly, the total
vector population is divided into susceptible, exposed, and infectious groups. By studying which
elements affect the transmission between a host and a vector and how the disease is transmitted
from one compartment to another, we develop a system of nonlinear differential equations that
describes the epidemiology of vector-borne disease.
In chapter 2, we analyzed the vertically transmitted deterministic vector-borne disease
model. We find the disease free equilibrium point E0 by setting the number of exposed, infecti-
ous, and treated groups is equal to zero. Then, we calculated the basic reproduction number
( or the epidemiology threshold) R0 using the Next-Generation Approach which is finding the
spectral radius of the next-generation matrix. We proved that if R0 < 1, the disease free equi-
librium points is locally and globally asymptotically stable and the disease is extinguished. If
R0 > 1, then the disease free equilibrium point is unstable. We also studied the sensitivity
of R0 . We studied that which parameter make the value of R0 change most and least using
the normalized forward sensitivity index γpu with a variable u that depends differentiably on a
parameter p. We found that the change of host recovery rate γ, a factor for density dependent
maturation of mosquitoes to adulthood δ1 , and the number of contacts between a host and a
vector α affect the value of R0 most and the change of disease-induced host death rate α and
vertical disease transmission rate in the vector ζ2 least. We do the numerical experiment for the

53
different value of parameters to explain the analytic results. We had the numerical experiments
for R0 < 1 and R0 > 1 and we see that for R0 < 1, the disease die out and for R0 > 1, the
disease persist.
In chapter 3, we studied the optimal control for the vertically transmitted deterministic
vector-borne disease model. We consider the two optimal control functions ,The prevention u1
and the treatment u2 , which are piecewise continuous on a certain time interval. We introduce
these two controls into the deterministic model in chapter 2, which u1 is related to susceptible,
exposed host groups, susceptible, and exposed vector groups and u2 is related to infectious and
treated host groups. We find u1 and u2 by minimizing the cost J(u1 , u2 ) which minimizes the
number of exposed and infected host groups and the cost for the prevention and the treatment
and maximizes the number of susceptible and treatment host groups at the final time step. We
shoe the existence of the optimal control functions using Carathodory’s existence theorem. We
calculate the optimal controls by finding the optimal system. We formulate Hamiltonian to
find the optimal system with the deterministic model, the integrand of the cost functional with
the adjoint variables which is similar to a Lagrange multiplier using Pontry-yagins Maximum
Principle. We find the explicit formula for u1 and u2 in terms of the status variables and the
adjoint variables. For the numerical simulation, we use the forward-backward sweep method
since the status system which is the deterministic model with the optimal controls has the ini-
tial condition and the system of adjoint variables has the final condition. We solve the status
system and the adjoint system using the Runge-Kutta method in 3-dimensions. In the numerical
simulation, we focus on the effectiveness of the controls u1 and u2 . We choose the parameters
based on the dengue virus. We compared the result between when we control the disease and
the case when we don’t. We see that when the controls are included, the disease die out faster
than the result without the controls and we archive more the number of susceptible and treated
host groups faster than without controls. We also see the effectiveness of each controls. The
effectiveness of the prevention is higher than the treatment.
In the future, we will study for the periodic model and the stochastic model of the deter-
ministic model in chapter 2. The transmission of a disease is affected by the season and the
temperature which is periodic. We expect that we can get the better result which fits the result

54
from the real life. Also, we can study the transmission in the small number of host and vector
by considering the stochastic model.

55
 References

[1] Ben Adams and Michael Boots. How important is vertical transmission in mosquitoes for
the persistence of dengue? insights from a mathematical model. Epidemics, 2(1):1–10,
2010.

[2] John F Anderson and Andy J Main. Importance of vertical and horizontal transmission of
west nile virus by culex pipiens in the northeastern united states. The Journal of infectious
diseases, 194(11):1577–1579, 2006.

[3] Sebastian Anita, Vincenzo Capasso, and Viorel Arnautu. An Introduction to Optimal Con-
trol Problems in Life Sciences and Economics: From Mathematical Models to Numerical
Simulation with MATLAB®. Springer, 2011.

[4] Kendall E Atkinson. An introduction to numerical analysis. John Wiley & Sons, 2008.

[5] Vadim Azhmyakov, Ruthber Rodriguez Serrezuela, Angela Magnolia Rios Gallardo, and
Winston Gerardo Vargas. An approximations based approach to optimal control of swit-
ched dynamic systems. Mathematical Problems in Engineering, 2014, 2014.

[6] Nicolas Bacaër and Souad Guernaoui. The epidemic threshold of vector-borne diseases
with seasonality. Journal of mathematical biology, 53(3):421–436, 2006.

[7] Shahida Baqar, Curtis G Hayes, James R Murphy, and Douglas M Watts. Vertical trans-
mission of west nile virus by culex and aedes species mosquitoes. The American journal
of tropical medicine and hygiene, 48(6):757–762, 1993.

[8] B.J. Beaty and W.C. Marquardt. The Biology of Disease Vectors. University Press of
Colorado, 1996.

56
 [9] Abraham Berman and Robert J Plemmons. Nonnegative matrices in the mathematical
sciences. SIAM, 1994.

[10] John T Betts. Practical methods for optimal control and estimation using nonlinear pro-
gramming, volume 19. Siam, 2010.

[11] Kbenesh Blayneh. Vertically transmitted vector-borne diseases and the effects of extreme
temperature. International Journal of Applied Mathematics, 30(2):177–209, 2017.

[12] Kbenesh Blayneh, Yanzhao Cao, and Hee-Dae Kwon. Optimal control of vector-borne
diseases: treatment and prevention. Discrete and Continuous Dynamical Systems Series
B, 11(3):587–611, 2009.

[13] Kbenesh W Blayneh, Abba B Gumel, Suzanne Lenhart, and Tim Clayton. Backward
bifurcation and optimal control in transmission dynamics of west nile virus. Bulletin of
mathematical biology, 72(4):1006–1028, 2010.

[14] Kbenesh W Blayneh and Jemal Mohammed-Awel. Insecticide-resistant mosquitoes and

malaria control. Mathematical biosciences, 252:14–26, 2014.

[15] Francesco Borrelli. Constrained optimal control of linear and hybrid systems, volume
290. Springer, 2003.

[16] CF Bosio, RE Thomas, PR Grimstad, and KS Rai. Variation in the efficiency of verti-
cal transmission of dengue-1 virus by strains of aedes albopictus (diptera: Culicidae).
Journal of medical entomology, 29(6):985–989, 1992.

[17] JC Butcher. On the convergence of numerical solutions to ordinary differential equations.

Mathematics of Computation, 20(93):1–10, 1966.

[18] John C Butcher. A multistep generalization of runge-kutta methods with four or five
stages. Journal of the ACM (JACM), 14(1):84–99, 1967.

[19] Liming Cai and Xuezhi Li. Analysis of a simple vector-host epidemic model with direct
transmission. Discrete Dynamics in Nature and Society, 2010, 2010.

57
[20] Dean A Carlson, Alain B Haurie, and Arie Leizarowitz. Infinite horizon optimal control:
deterministic and stochastic systems. Springer Science & Business Media, 2012.

[21] Jack Carr. Applications of centre manifold theory, volume 35. Springer Science & Busi-
ness Media, 2012.

[22] C Castillo-Chávez, Z Feng, and W Huang. On the computation of ro and its role on global
stability in mathematical approaches for emerging and re-emerging infectious diseases,
part i, ima, 125.

[23] Carlos Castillo-Chavez and Baojun Song. Dynamical models of tuberculosis and their
applications. Mathematical biosciences and engineering, 1(2):361–404, 2004.

[24] Benoit Chachuat. Nonlinear and dynamic optimization: From theory to practice. Techni-
cal report, 2007.

[25] Carmen Chicone. Ordinary differential equations with applications, volume 34. Springer
Science & Business Media, 2006.

[26] Nakul Chitnis, James M Hyman, and Jim M Cushing. Determining important parameters
in the spread of malaria through the sensitivity analysis of a mathematical model. Bulletin
of mathematical biology, 70(5):1272, 2008.

[27] Gerardo Chowell and Fred Brauer. The basic reproduction number of infectious diseases:
computation and estimation using compartmental epidemic models. In Mathematical and
statistical estimation approaches in epidemiology, pages 1–30. Springer, 2009.

[28] Frank H Clarke. Optimization and nonsmooth analysis, volume 5. Siam, 1990.

[29] Odo Diekmann, Hans Heesterbeek, and Tom Britton. Mathematical tools for understan-
ding infectious disease dynamics. Princeton University Press, 2012.

[30] Odo Diekmann and Johan Andre Peter Heesterbeek. Mathematical epidemiology of in-
fectious diseases: model building, analysis and interpretation, volume 5. John Wiley &
Sons, 2000.

58
[31] Lourdes Esteva and Cristobal Vargas. A model for dengue disease with variable human
population. Journal of mathematical biology, 38(3):220–240, 1999.

[32] Zhilan Feng and Jorge X Velasco-Hernández. Competitive exclusion in a vector-host

model for the dengue fever. Journal of mathematical biology, 35(5):523–544, 1997.

[33] Wendell H Fleming and Raymond W Rishel. Deterministic and stochastic optimal control,
volume 1. Springer Science & Business Media, 2012.

[34] Martin Grunnill and Michael Boots. How important is vertical transmission of dengue
viruses by mosquitoes (diptera: Culicidae)? Journal of medical entomology, 53(1):1–19,
2015.

[35] AB Gumel. Causes of backward bifurcations in some epidemiological models. Journal of

Mathematical Analysis and Applications, 395(1):355–365, 2012.

[36] Misha Guysinsky, Boris Hasselblatt, and Victoria Rayskin. Differentiability of the
hartman-grobman linearization. Discrete and Continuous Dynamical Systems, 9(4):979–
984, 2003.

[37] J.K. Hale. Ordinary Differential Equations: Pure and Applied Mathematics. (Wiley-
Interscience) 21. Pure and applied mathematics, 21. Wiley-Interscience, 1969.

[38] Jane M Heffernan, Robert J Smith, and Lindi M Wahl. Perspectives on the basic repro-
ductive ratio. Journal of the Royal Society Interface, 2(4):281–293, 2005.

[39] Herbert W Hethcote. The mathematics of infectious diseases. SIAM review, 42(4):599–
653, 2000.

[40] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press,
2012.

[41] James Holland Jones. Notes on r0 . Department of Anthropology Science, 2007.

[42] Hem Raj Joshi. Optimal control of an hiv immunology model. Optimal control applicati-
ons and methods, 23(4):199–213, 2002.

59
[43] E Jung, Suzanne Lenhart, and Z Feng. Optimal control of treatments in a two-strain
tuberculosis model. Discrete and Continuous Dynamical Systems Series B, 2(4):473–
482, 2002.

[44] M Keeling and Pejman Rohani. Modeling infectious diseases in humans and animals.
Clinical Infectious Diseases, 47:864–6, 2008.

[45] Denise Kirschner, Suzanne Lenhart, and Steve Serbin. Optimal control of the chemother-
apy of hiv. Journal of mathematical biology, 35(7):775–792, 1997.

[46] Joseph P La Salle. The stability of dynamical systems. SIAM, 1976.

[47] Suzanne Lenhart and John T Workman. Optimal control applied to biological models.
Crc Press, 2007.

[48] Frank L Lewis, Draguna Vrabie, and Vassilis L Syrmos. Optimal control. John Wiley &
Sons, 2012.

[49] D. L. Lukas. Differential equation classical to controlled mathematical in science and

engineering, 1982.

[50] Stefan Ma and Yingcun Xia. Mathematical understanding of infectious disease dynamics,
volume 16. World Scientific, 2009.

[51] Xuerong Mao and Chenggui Yuan. Stochastic differential equations with Markovian swit-
ching. Imperial College Press, 2006.

[52] William H Marquardt. Biology of disease vectors. Elsevier, 2004.

[53] Maia Martcheva. An introduction to mathematical epidemiology, volume 61. Springer,

2015.

[54] Michael McAsey, Libin Mou, and Weimin Han. Convergence of the forward-backward
sweep method in optimal control. Computational Optimization and Applications,
53(1):207–226, 2012.

60
[55] Yukihiko Nakata and Toshikazu Kuniya. Global dynamics of a class of seirs epidemic
models in a periodic environment. Journal of Mathematical Analysis and Applications,
363(1):230–237, 2010.

[56] Calistus N Ngonghala, Sara Y Del Valle, Ruijun Zhao, and Jemal Mohammed-Awel.
Quantifying the impact of decay in bed-net efficacy on malaria transmission. Journal
of theoretical biology, 363:247–261, 2014.

[57] Gideon A Ngwa and William S Shu. A mathematical model for endemic malaria with va-
riable human and mosquito populations. Mathematical and Computer Modelling, 32(7-
8):747–763, 2000.

[58] Muhammad Ozair, Abid Ali Lashari, Il Hyo Jung, and Kazeem Oare Okosun. Stability
analysis and optimal control of a vector-borne disease with nonlinear incidence. Discrete
Dynamics in Nature and Society, 2012, 2012.

[59] Lev Semenovich Pontryagin. Mathematical theory of optimal processes. CRC Press,
1987.

[60] Vadrevu Sree Hari Rao and Ravi Durvasula. Dynamic models of infectious diseases,
volume 1. Springer, 2013.

[61] Robert C Reiner, T Alex Perkins, Christopher M Barker, Tianchan Niu, Luis Fernando
Chaves, Alicia M Ellis, Dylan B George, Arnaud Le Menach, Juliet RC Pulliam, Donal
Bisanzio, et al. A systematic review of mathematical models of mosquito-borne pathogen
transmission: 1970–2010. Journal of The Royal Society Interface, 10(81):20120921,
2013.

[62] Helena Sofia Rodrigues, M Teresa T Monteiro, and Delfim FM Torres. Optimal control
and numerical software: an overview. arXiv preprint arXiv:1401.7279, 2014.

[63] Md Samsuzzoha, Manmohan Singh, and David Lucy. Uncertainty and sensitivity analysis
of the basic reproduction number of a vaccinated epidemic model of influenza. Applied
Mathematical Modelling, 37(3):903–915, 2013.

61
[64] Lisa Sattenspiel. Modeling the spread of infectious disease in human populations. Ame-
rican Journal of Physical Anthropology, 33(S11):245–276, 1990.

[65] M. W. Service. Blood-sucking insects, vectors of disease / Michael W. Service. E. Arnold

London ; Baltimore, Md., USA, 1986.

[66] Pauline Van den Driessche and James Watmough. Reproduction numbers and sub-
threshold endemic equilibria for compartmental models of disease transmission. Mat-
hematical biosciences, 180(1):29–48, 2002.

[67] Paul Waltman. A second course in elementary differential equations. Courier Corpora-
tion, 2004.

[68] Wendi Wang and Xiao-Qiang Zhao. Threshold dynamics for compartmental epide-
mic models in periodic environments. Journal of Dynamics and Differential Equations,
20(3):699–717, 2008.

[69] Xuezhong Wang. Solving optimal control problems with matlab: Indirect methods. ISE
Dept., NCSU, Raleigh, NC, 27695, 2009.

[70] Hui-Ming Wei, Xue-Zhi Li, and Maia Martcheva. An epidemic model of a vector-borne
disease with direct transmission and time delay. Journal of Mathematical Analysis and
Applications, 342(2):895–908, 2008.

[71] WHO. Dengue and severe dengue. [Link]

fact-sheets/detail/dengue-and-severe-dengue.

[72] WHO. Vector-borne diseases. [Link]

fact-sheets/detail/vector-borne-diseases.

[73] Jerzy Zabczyk. Mathematical control theory: an introduction. Springer Science & Busi-
ness Media, 2009.

62
 Appendices

.1 Matlab Code

Matlab code for the numerical simulation of the optimal control. main.m

1 clc
2 clear
3
4 epsilon = 1/12;
5 gamma = 1/15;
6 beta = 0.2;
7 theta1 = 0.0082;
8 theta2 = 0.0289;
9 mu = 1.01/(70*365);
10 Lambda = 0.379/(70*365);
11 rho = 205;
12 delta0 = 1050;
13 r= 1/7;
14 r0 = 0.04;
15 delta1 = 0.0399;
16 d = 1/10;
17 psi = 0.0014;
18 zeta1 = 0;
19 zeta2 = 0.67;
20 alpha = 0.0238;

63
21 phi = 3;
22 A1 = 1;
23 A2 = 1;
24 B1 = 50;
25 B2 = 50;
26 T = 365; %final time
27 N = 200;
28 x0 = [100; 20; 20; 10; 1000; 20; 30]; % initial value of x
29 u0 = [0; 0]; % initial value of u
30 lambdafinal =[−0.1;0;0;−0.1;0;0;0]; %final values of lambda
31 t = linspace(0,T,N+1);
32 h = T/N;
33 x = zeros(7,N+1); %initialize x
34 u = zeros(2,N+1); %initialize u
35 lambda = zeros(7,N+1); %initialize x
36 x (:,1) = x0; %assign the x0
37 u (:,1) = u0; %assign the u0
38 lambda(:,N+1) = lambdafinal; %assign the lambdafinal
39
40 %R0
41 %1/2*(zeta1 * Lambda*d/(k1*k2) + zeta2 * delta1*epsilon/(gamma*k4))...
42 %+sqrt((1/2 *(zeta1*Lambda*d/(k1*k3)−zeta2*delta1*epsilon/(gamma*k4)))ˆ2...
43 %+phiˆ2 * (beta*epsilon*(theta1*d+theta2*k3)*(delta0/(gamma−delta1)))/(gamma*
k1*k3*k4*(rho/(mu−Lambda))))
44
45 parameters=[epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0, ...
46 delta1 , d, psi , zeta1, zeta2,alpha,phi,A1,A2,B1,B2];
47

64
48 k = 1;%counter of the iteration
49 delta = 0.001; %error bound
50 test = −1; %error
51 while ( test<0 && k<1000)
52 oldx = x;
53 oldu = u;
54 oldlambda = lambda;
55
56 %forward Runge−Kutta 4th with 3 input algorithm for state
57 for i = 1:N
58 k1(1:7,1) = state(t( i ) ,x (:, i ) ,u (:, i ) ,parameters);
59 k2(1:7,1) = state(t( i )+h/2,x(:,i)+h*k1/2,(u(:,i)+u(:,i+1))/2,parameters);
60 k3(1:7,1) = state(t( i )+h/2,x(:,i)+h*k2/2,(u(:,i)+u(:,i+1))/2,parameters);
61 k4(1:7,1) = state(t( i )+h,x(:, i )+h*k3,u(:,i+1),parameters);
62 x (:, i+1) = x(:,i)+ (h/6) * (k1+2*k2+2*k3+k4);
63 end
64
65 %backward Runge−Kutta 4th with 3 input algorithm for adjoint
66 for i = 1:N
67 j = N+2−i;
68 k1(1:7,1) = adjoint(t(j ) , lambda(:,j) , x (:, j) ,u (:, j ) ,parameters);
69 k2(1:7,1) = adjoint(t(j )−h/2, lambda(:,j)−h*k1/2,(x(:,j)+x(:,j−1))/2 ...
70 ,( u (:, j)+u(:,j−1))/2,parameters);
71 k3(1:7,1) = adjoint(t(j )−h/2,lambda(:,j)−h*k2/2,(x(:,j)+x(:,j−1))/2 ...
72 ,( u (:, j)+u(:,j−1))/2,parameters);
73 k4(1:7,1) = adjoint(t(j )−h/2, lambda(:,j)−h*k3/2,x(:,j−1),u(:,j−1),
parameters);
74 lambda(:,j−1) = lambda(:,j) − (h/6) *(k1+2*k2+2*k3+k4);

65
 75 end
76
77 %Find controls
78 for i = 1:N+1
79 u(1, i ) = max(0,min(1,(beta*lambda(1,i)*(−x(1,i))*x(7,i)*phi...
80 +beta*lambda(2,i)*x(1,i)*x(7,i)*phi ...
81 +(lambda(6,i)−lambda(5,i))*x(5,i)*phi*(theta2*x(2,i) ...
82 +theta1*x(3,i)))/(2*B1*(x(1,i)+x(2,i)+x(3,i)+x(4,i))))) ;
83 u(2, i ) = max(0,min(1,(lambda(3,i)−lambda(4,i))*r0*x(3,i)/(2*B2)));
84 end
85
86 %updates Control
87 c = 0.8;
88 u = (1−c)*u + c*oldu;
89 % for i=1:N+1
90 % if (u(1, i ) > oldu(1,i))
91 % u(1, i ) = (1 − c) + oldu(1,i)*c;
92 % else
93 % u(1, i ) = oldu(1,i)*c;
94 % end
95 % if (u(2, i ) > oldu(2,i))
96 % u(2, i ) = (1 − c) + oldu(2,i)*c;
97 % else
98 % u(2, i ) = oldu(2,i)*c;
99 % end
100 % end
101
102 %error bewteen old and new

66
103 tempu = min(delta*sum(abs(u),2)−sum(abs(oldu−u),2));
104 tempx = min(delta*sum(abs(x),2)−sum(abs(oldx−x),2));
105 templambda = min(delta*sum(abs(lambda),2)−sum(abs(oldlambda−lambda),2));
106 test = min(tempu,min(tempx,templambda))
107
108 %The cost
109 %trapz(A1*oldx(2,:)+A2*oldx(3,:)+B1*oldu(1,:).ˆ2+B2*oldu(2,:).ˆ2)...
110 % −0.1*oldx(1,N+1)− 0.1*oldx(4,N+1)
111 %trapz(A1*x(2,:)+A2*x(3,:)+B1*u(1,:).ˆ2+B2*u(2,:).ˆ2) ...
112 % −0.1*x(1,N+1)−0.1*x(4,N+1)
113
114 k=k+1;
115 end
116
117
118 plot(t ,x (3,:) , t ,x (7,:) )
119 title (' infectious ')
120 legend('host','vector')

67
 The Matlab code for the system of the status (state.m).

1 function dxdt = state(t,x,u,parameters)

2 parameters = num2cell(parameters);
3 [epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0,
...
4 delta1, d, psi, zeta1, zeta2,alpha,phi,A1,A2,B1,B2] = deal(
parameters{:});
5
6 dxdt = zeros(7,1);
7 dxdt(1) = rho + Lambda*(x(1)+x(2)+x(4)) + Lambda*(1-zeta1)*x
(3)+psi*x(4)...
8 -beta * phi*x(7)*x(1)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) - mu *
x(1);
9 dxdt(2) = beta*phi*x(7)*x(1)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) +
Lambda*zeta1*x(3)-d*x(2)-mu*x(2);
10 dxdt(3) = d*x(2) - (r + alpha + mu+r0*u(2))*x(3);
11 dxdt(4) = (r+r0*u(2))* x(3)- (psi + mu)*x(4);
12 dxdt(5) = delta0 + delta1*(x(5)+x(6)) + delta1*(1-zeta2)*x(7)
...
13 - phi*theta1*x(3)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) ...
14 - phi*theta2*x(2)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4)) -
gamma*x(5);
15 dxdt(6) = phi*theta1*x(3)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4))
...
16 + phi*theta2*x(2)*x(5)*(1-u(1))/(x(1)+x(2)+x(3)+x(4))...
17 + delta1*zeta2*x(7)-epsilon*x(6)-gamma*x(6);
18 dxdt(7) = epsilon*x(6)-gamma*x(7);

68
 The Matlab code for the system of the adjoint (adjoint.m).

1 function dlambdadt = adjoint(t,lambda,x,u,parameters)

2 parameters = num2cell(parameters);
3 [epsilon,gamma,beta,theta1,theta2,mu,Lambda,rho,delta0,r, r0,
...
4 delta1, d, psi, zeta1, zeta2,alpha,phi,A1,A2,B1,B2] = deal(
parameters{:});
5
6 dlambdadt = zeros(7,1);
7 dlambdadt(1) = -(lambda(1)*(Lambda-mu+beta*(1-u(1))*x(1)*x(7)*
phi/(x(1)+x(2)+x(3)+x(4))ˆ2 ...
8 -beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4))) ...
9 +lambda(2)*(beta*(1-u(1))*x(7)*phi/(x(1)+x(2)+x(3)+x(4))
...
10 -beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2)...
11 +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x
(4))ˆ2 ...
12 +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2)
...
13 +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+
x(4))ˆ2 ...
14 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2))
;
15 dlambdadt(2) = -(A1+d*lambda(3) ...
16 +lambda(2)*(-d-mu-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x
(3)+x(4))ˆ2) ...
17 +lambda(1)*(Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x
(3)+x(4))ˆ2) ...

69
18 +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x
(4))ˆ2 ...
19 +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2
...
20 -theta2*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))) ...
21 +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+
x(4))ˆ2 ...
22 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2
...
23 +theta2*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))));
24 dlambdadt(3) = -(A2+lambda(3)*(-alpha-mu-r0*u(2)-r)+lambda(4)

*(r0*u(2)+r)...
25 +lambda(1)*((1-zeta1)*Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x
(1)+x(2)+x(3)+x(4))ˆ2)...
26 +lambda(2)*(zeta1*Lambda-beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+
x(2)+x(3)+x(4))ˆ2)...
27 +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x
(4))ˆ2 ...
28 +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2
...
29 -theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4)))...
30 +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+
x(4))ˆ2 ...
31 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2
...
32 +theta1*(1-u(1))*x(5)*phi/(x(1)+x(2)+x(3)+x(4))));
33 dlambdadt(4) = -(lambda(4)*(-mu-psi)+...

70
34 lambda(1)*(Lambda+beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x
(3)+x(4))ˆ2+psi)...
35 -lambda(2)*(beta*(1-u(1))*x(1)*x(7)*phi/(x(1)+x(2)+x(3)+x
(4))ˆ2) ...
36 +lambda(5)*(theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+x
(4))ˆ2 ...
37 +theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2)
...
38 +lambda(6)*(-theta2*(1-u(1))*x(2)*x(5)*phi/(x(1)+x(2)+x(3)+
x(4))ˆ2)...
39 -theta1*(1-u(1))*x(3)*x(5)*phi/(x(1)+x(2)+x(3)+x(4))ˆ2);
40 dlambdadt(5) = -(lambda(5)*(-gamma+delta1-theta2*(1-u(1))*x(2)

*phi/(x(1)+x(2)+x(3)+x(4)) ...
41 -theta1*(1-u(1))*x(3)*phi/(x(1)+x(2)+x(3)+x(4))) ...
42 +lambda(6)*(theta2*(1-u(1))*x(2)*phi/(x(1)+x(2)+x(3)+x(4))
...
43 +theta1*(1-u(1))*x(3)*phi/(x(1)+x(2)+x(3)+x(4))));
44 dlambdadt(6) = -(lambda(6)*(-gamma-epsilon)+delta1*lambda(5)+
lambda(7)*epsilon);
45 dlambdadt(7) = -(-gamma*lambda(7) +delta1*(1-zeta2)*lambda(5)+
delta1*zeta2*lambda(6)...
46 -beta*lambda(1)*(1-u(1))*x(1)*phi/(x(1)+x(2)+x(3)+x(4))...
47 +beta*lambda(2)*(1-u(1))*x(1)*phi/(x(1)+x(2)+x(3)+x(4)));

71