A relação entre a epidemiologia, a modelação matemática e as ferramentas computacionais permite construir e testar teorias sobre o desenvolvimento e combate de uma doença. Esta tese tem como motivação o estudo de modelos epidemiológicos aplicados a doenças infeciosas numa perspetiva de ControloÓtimo, dando particular relevância ao Dengue. Sendo uma doença tropical e subtropical transmitida por mosquitos, afecta cerca de 100 milhões de pessoas por ano, eé considerada pela Organização Mundial de Saúde como uma grande preocupação para a saúde pública. Os modelos matemáticos desenvolvidos e testados neste trabalho, baseiamse em equações diferenciais ordinárias que descrevem a dinâmica subjacentè a doença nomeadamente a interação entre humanos e mosquitos.É feito um estudo analítico dos mesmos relativamente aos pontos de equilíbrio, sua estabilidade e número básico de reprodução. A propagação do Dengue pode ser atenuada através de medidas de controlo do vetor transmissor, tais como o uso de inseticidas específicos e campanhas educacionais. Como o desenvolvimento de uma potencial vacina tem sido uma aposta mundial recente, são propostos modelos baseados na simulação de um hipotético processo de vacinação numa população.
Tendo por base a teoria de ControloÓtimo, são analisadas as estratégiasótimas para o uso destes controlos e respetivas repercussões na redução/erradicação da doença aquando de um surto na população, considerando uma abordagem bioeconómica. Os problemas formulados são resolvidos numericamente usando métodos diretos e indiretos. Os primeiros discretizam o problema reformulando-o num problema de optimização não linear. Os métodos indiretos usam o Princípio do Máximo de Pontryagin como condição necessária para encontrar a curvá otima para o respetivo controlo. Nestas duas estratégias utilizam-se vários pacotes de software numérico. Ao longo deste trabalho, houve sempre um compromisso entre o realismo dos modelos epidemiológicos e a sua tratabilidade em termos matemáticos. Introduction "Mathematical biology is a fast-growing, well-recognized, albeit not clearly defined, subject and is, to my mind, the most exciting modern application of mathematics."
-J. D. Murray, Mathematical Biology, 2002 Epidemiology has become an important issue for modern society. The relationship between mathematics and epidemiology has been increasing. For the mathematician, epidemiology provides new and exciting branches, while for the epidemiologist, mathematical modeling offers an important research tool in the study of the evolution of diseases.
In 1760, a smallpox model was proposed by Daniel Bernoulli and is considered by many authors the first epidemiological mathematical model. Theoretical papers by Kermack and McKendrinck, between 1927 and 1933 about infectious disease models, have had a great influence in the development of mathematical epidemiology models [93] . Most of the basic theory had been developed during that time, but the theoretical progress has been steady since then [14] . Mathematical models are being increasingly used to elucidate the transmission of several diseases. These models, usually based on compartment models, may be rather simple, but studying them is crucial in gaining important knowledge of the underlying aspects of the infectious diseases spread out [63] , and to evaluate the potential impact of control programs in reducing morbidity and mortality.
After the Second World War, the strategy of public health has been focusing on the control and elimination of the organisms that cause the diseases. The appearance of new antibiotics and vaccines brought a positive perspective of the diseases eradication. However, factors such as resistance to the medicine by the microorganisms, demographic evolution, accelerated urbanization, increased travelling and climate change, led to new diseases and the resurgence of old ones. In 1981, the human immunodeficiency virus (HIV) appears and since then, become as important sexually transmitted disease throughout the world [64] . Futhermore, malaria, tuberculosis, dengue and yellow fever have re-emerged and, as a result of climate changes, has been spreading into new regions [64] .
Recent years have seen an increasing trend in the representation of mathematical models in publications in the epidemiological literature, from specialist journals of medicine, biology and mathematics to the highest impact generalist journals [50] , showing the importance of interdisciplinary. Their role in comparing, planning, implementing and evaluating various control programs is of major importance for public health decision makers. This interest has been reinforced by the recent examples of SARS -Severe Acute Respiratory Syndrome -epidemic in 2003 and Influenza pandemic in 2009.
Although chronic diseases, such as cancer and heart diseases have been receiving more attention in developed countries, infectious diseases are still important and cause suffering and mortality in developing countries. These, remain a serious medical burden all around the world with 15 million deaths per year estimated to be directly related to infectious diseases [71] .
The successful containment of the emerging diseases is not just linked to medical infrastructure but also on the capacity to recognize its transmission characteristics and apply optimal medical and logistic policies. Public health often asks information such as [64] : how many people will be infected, how many require hospitalization, what is the maximum number of people ill at a given time and how long will the epidemic last. As a result, it is necessary an ever-increasing capacity for a rapid response.
Education, vaccination campaigns, preventive drugs administration and surveillance programs, are all examples of prevention methods that authorities must consider for disease prevention. Whenever the disease declares itself, the emergency interventions such as disinfectants, insecticide application, mechanical controls and quarantine measures must be considered. Intervention strategies can be modelled with the goal of understanding how they will influence the disease's battle.
As financial resources are limited, there is a pressing need to optimize investments for disease prevention and fight. Traditionally, the study of disease dynamics has been focused on identifying the mechanisms responsible for epidemics but has taken little into account economic constraints in analyzing control strategies. On the other hand, economic models have given insight into optimal control under constraints imposed by limited resources, but they are frequently ignored by the spatial and temporal dynamics of the disease. Therefore, progress requires a combination of epidemiological and economic factors for modelling what until here tended to remain separate. More recently, bioeconomic approaches to disease management have been advocated, since infectious diseases can be modelled thinking that the limited resources involved require trade-offs. Finding the optimal strategy depends on the balance of economic and epidemiological parameters that reflect the nature of the host-pathogen system and the efficiency of the control method.
The main goal of this thesis is to formulate epidemiological models, giving a special importance to Dengue disease. Moreover, it is our aim to frame the disease management question into an optimal control problem requiring the maximization/minimization of some objective function that depends on the infected individuals (biological issues) and control costs (economic issues), given some initial conditions. This way, will allow us to propose practical control measures to the authorities to assess and forecast the disease burden, such as an attack rate, morbidity, hospitalization and mortality.
The thesis is composed by two parts. The First Part, comprising chapters 1 to 3, gives a mathematical background to support the original results presented in the Second Part, that is composed by chapters 4 to 7. In Chapter 1, the definition of Optimal Control Problem, its possible versions and the adapted first order necessary conditions based on the Pontryagin Maximum Principle are introduced. Simple examples are chosen to exemplify the mathematical concepts.
With the increasing of variables and complexity, Optimal Control problems can no longer be solved analytically and numerical methods are required. For this purpose, in Chapter 2, direct and indirect methods are presented for their resolution. Direct methods consist in the discretization of the Optimal xii Control problem, reducing it to a nonlinear constrained optimization problem. Indirect methods are based on the Pontryagin Maximum Principle, which in turn reduces the problem to a boundary value problem. For each approach, the software packages used in this thesis are described.
In Chapter 3, the basic building blocks of most epidemiological models are reviewed: SIR (composed by Susceptible-Infected-Recovered) and SIS models (Susceptible-Infected-Susceptible). For these, it is possible to develop some analytical results which can be useful in the understanding of simple epidemics. Taking this as the basis, we discuss the dynamics of other compartmental models, bringing more complex realities, such as those with exposed or carrier classes.
The Second Part of the thesis contains the original results and is focused on Dengue Fever disease. Dengue is a vector borne disease, caused by a mosquito from the Aedes family. It is mostly found in tropical and sub-tropical climates, mostly in urban areas. It can provoke a severe flu-like illness, and sometimes, in severe cases can be lethal. According to the World Health Organization about 40% of the world's population is now at risk [141] .
The main reasons for the choice of this particular disease are:
• the importance of this disease around the world, as well as the challenges of its transmission features, prevention and control measures;
• two Portuguese-speaking countries (Brazil and Cape Verde) have already experience with Dengue, and in the last one, a first outbreak occurred during the development of this thesis, which allowed the development of a groundbreaking work;
• the mosquito Aedes Aegypti, the main vector that transmits Dengue, is already in Portugal, on Madeira island [5] , which without carrying the disease, is considered a potential threat to public health and has been followed by the health authorities.
In Chapter 4 information about the mosquito, disease symptoms, and measures to fight Dengue are reported. An old problem related to Dengue is revisited and solved by different approaches. Finally, a numerical study is performed to compare different discretization schemes in order to analyze the best strategies for future implementations.
In Chapter 5, a SEIR+ASEI model is studied. The basic reproduction number and the equilibrium points are computed as well as their relationship with the local stability of the Disease Free Equilibrium. This model implements a control measure: adulticide. A study to find the best strategies available to apply the insecticide is made. Continuous and piecewise constant strategies are used involving the system of ordinary differential equations only or resorting to the Optimal Control theory.
Chapter 6 is concerned with a SIR+ASI model that incorporates three controls: adulticide, larvicide and mechanical control. A detailed discussion on the effects of each control, individually or together, on the development of the disease is given. An analysis of the importance of each control in the decreasing of the basic reproduction number is given. These results are strengthened when the optimal strategy for the model is calculated. Bioeconomic approaches, using distinct weights for the respective control costs and treatments for infected individuals have also been provided.
In Chapter 7 simulations for a hypothetical vaccine for Dengue are carried out. The features of the vaccine are unknown because the clinical trials are ongoing. Using models with a new compartment for vaccinated individuals, perfect and imperfect vaccines are studied aiming the analysis of the repercussions of the vaccination process in the morbidity and/or eradication of the disease. Then, an optimal control approach is studied, considering the vaccination not as a new compartment, but as a measure control in fighting the disease.
Finally, the main conclusions are reported and future directions of research are pointed out.
xiv Part I State of the art 1 Chapter 1
The optimal control definition and its possible formulations are introduced, followed by some examples related to epidemiological models. The Pontryagin Maximum Principle is presented with the aim of finding the best control policy.
Optimal Control (OC) is the process of determining control and state trajectories for a dynamic system over a period of time in order to minimize a performance index [16] .
Historically, OC is an extension of the calculus of variations. In the seventeenth century, the first formal results of calculus of variations can be found. Johann Bernoulli challenged other famous contemporary mathematicians -such as Newton, Leibniz, Jacob Bernoulli, L'Hôpital and von Tschirnhaus -with the Brachistochrone problem: "if a small object moves under the influence of gravity, which part between two fixed points enables it to make the trip in the shortest time?"
Other specific problems were solved and a general mathematical theory was developed by Euler and Lagrange. The most fruitful applications of the calculus of variations have been to theoretical physics, particularly in connection with Hamilton's principle or the Principle of Least Action. Early applications to economics appeared in the late 1920s and early 1930s by Ross, Evans, Hottelling and Ramsey, with further applications published occasionally thereafter [129] .
The generalization of the calculus of variations to optimal control theory was strongly motivated by military applications and has developed rapidly since 1950. The decisive breakthrough was achieved by the Russian mathematician Lev S. Pontryagin (1908 Pontryagin ( -1988 and his co-workers (V. G. Boltyanskii, R. V. Gamkrelidz and E. F. Misshchenko) with the formulation and demonstration of the Pontryagin Maximum Principle [106] . This principle has provided research with suitable conditions for optimization problems with differential equations as constraints. The Russian team generalized variational problems by separating control and state variables and admitting control constraints. In such problems, OC gives equivalents results, as one would have expected. However, the two approaches differ and the OC approach sometimes affords insight into a problem that might be less readily apparent through the calculus of variations. OC is also applied to problems where the calculus of variations is not convenient, such as those involving constraints on the derivatives of functions [80] .
The theory of OC brought new approaches to Mathematics with Dynamic Programming. Introduced by R. E. Bellman, Dynamic Programming makes use of the principle of optimality and it is suitable for solving discrete problems, allowing for a significant reduction in the computation of the optimal controls (see [78] ). It is also possible to obtain a continuous approach to the principle of optimality that leads to the solution of a partial differential equation called the Hamilton-Jacobi-Bellman equation. This result allowed to bring new connections between the OC problem and the Lyapunov stability theory.
Before the arrival of the computer, only fairly simple the OC problems could be solved. The arrival of the computer age enabled the application of OC theory and some methods to many complex problems. Selected examples are as follows:
• Physical systems, such as stable performance of motors and machinery, robotics, optimal guidance of rockets [59, 90] ;
• Aerospace, including driven problems, orbits transfers, development of satellite launchers and recoverable problems of atmospheric reentry [12, 62] ;
• Economics and management, such as optimal exploitation of natural resources, energy policies, optimal investment of production strategies [92, 127] ;
• Biology and medicine, as regulation of physiological functions, plants growth, infectious diseases, oncology, radiotherapy [72, 73, 81, 95] .
Today, the OC theory is extensive and with several approaches. One can adjust controls in a system to achieve a goal, where the underlying system can include: ordinary differential equations, partial differential equations, discrete equations, stochastic differential equations, integro-difference equations, combination of discrete and continuous systems. In this work the goal is the OC theory of ordinary differential equations with time fixed.
A typical OC problem requires a performance index or cost functional (J[x(t), u(t)]), a set of state variables (x(t) ∈ X), a set of control variables (u(t) ∈ U ) in a time t, with t 0 ≤ t ≤ t f . The main goal consists in finding a piecewise continuous control u(t) and the associated state variable x(t) to maximize a given objective functional. The development of this chapter will be closely structured from Lenhart and Workman work [81] .
x(t f ) could be free, which means that the value of x(t f ) is unrestricted, or could be fixed, i.e, x(t f ) = x f .
For our purposes, f and g will always be continuously differentiable functions in all three arguments. We assume that the control set U is a Lebesgue measurable function. Thus, as the control(s) will always be piecewise continuous, the associated states will always be piecewise differentiable.
We have been focused on finding the maximum of a function. We can switch back and forth between maximization and minimization by simply negating the cost functional:
An OC problem can be presented in different ways, but equivalent, depending on the purpose or the software to be used.
There are three well known equivalent formulations to describe the OC problem, which are Lagrange (already presented in previous section), Mayer and Bolza forms [25, 144] .
Definition 2 (Bolza formulation). The Bolza formulation of the OC problem can be defined as
where φ is a continuously differentiable function.
Definition 3 (Mayer formulation). The Mayer formulation of the OC problem can be defined as Proof.
(2) ⇒ (1) To get the proof, we formulate the Bolza problem as one of Lagrange, using an extended state vector.
Let (x(·), u(·)) an admissible pair for the problem (1.2) and let z(t) = (x a (t),
. So, (z(·), u(·)) is an admissible pair for the Lagrange problem
Thus, the value of the functionals in both formulations matches.
(1) ⇒ (2) Conversely, each admissible pair (z(·), u(·)) for the problem (1.4) corresponds the pair (x(·), u(·)), where x(·) is composed by the last component of z, admissible for the problem (1.2) and matching the respective values of the functionals.
(2) ⇒ (3) For this statement we also need to use an extended state vector. Let be
Thus we have (z(·), u(·)) an admissible pair for the following Mayer problem:
This way, the values of the functional for both formulations are the same.
(3) ⇒ (2) Conversely, to all admissible pair (z(·), u(·)) for the Mayer problem (1.5) corresponds to an admissible pair (x(·), u(·)) for the Bolza problem (1.2) , where x(·) consists in the last component of z(·).
For the proof of the previous theorem it was not necessary to show that the Lagrange problem is equivalent to the Mayer formulation. However, in the second part of the thesis, and due to computational issues, some of the OC problems (usually presented in the Lagrange form) will be converted into the equivalent Mayer one. Hence, using a standard procedure is possible to rewrite the cost functional (cf. [82] ), augmenting the state vector with an extra component. So, the Lagrange formulation (1.1) can be rewritten as max
The necessary first order conditions to find the optimal control were developed by Pontryagin and his co-workers. This result is considered as one of the most important results of Mathematics in the 20th century.
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 function of several variables to be maximized or minimized.
Definition 4 (Hamiltonian). Let the previous OC problem considered in (1.1). The function H(t, x(t), u(t), λ(t)) = f (t, x(t), u(t)) + λ(t)g(t, x(t), u(t)) If u * (t) and x * (t) are optimal for problem (1.1), then there exists a piecewise differentiable adjoint variable λ(t) such that
for all controls u at each time t, where H is the Hamiltonian previously defined and
Proof. The proof of this theorem is quite technical and we opted to omit it. The original Pontryagin's text [106] or Clarke's book [29] are good references to find the proof.
Remark 1. The last condition, λ(t f ) = 0, called transversality condition, is only used when the OC problem does not have terminal value in the state variable, i.e., x(t f ) is free.
This principle converted the problem of finding a control which maximizes the objective functional subject to the state ODE and initial condition into the problem of optimizing the Hamiltonian pointwise. As consequence, with this adjoint equation and Hamiltonian, we have
at u * for each t, namely, the Hamiltonian has a critical point; usually this condition is called optimality condition. Thus to find the necessary conditions, we do not need to calculate the integral in the objective functional, but only use the Hamiltonian.
Here is presented a simple example to illustrate this principle.
Example 1 (from [97] ). Consider the OC problem:
The calculus of this OC problem can be done by steps.
Step 1 -Form the Hamiltonian for the problem. The Hamiltonian can be written as:
Step 2 -Write the adjoint differential equation, the optimality condition and transversality boundary condition (if necessary). Try to eliminate u * by using the optimality equation H u = 0, i.e., solve for u * in terms of x * and λ.
Using the Hamiltonian to find the differential equation of the adjoint λ, we obtained
The optimality condition is given by
In this way we obtain an expression for the OC:
As the problem has just an initial condition for the state variable, it is necessary to calculate the transversality condition: λ(2) = 0.
Step 3 -Solve the set of two differential equations for x * and λ with the boundary conditions, replacing u * in the differential equations by the expression for the optimal control from the previous step.
By the adjoint equation λ = −1 − λ and the transversality condition λ(2) = 0 we have
Hence, the optimality condition leads to u * = −λ ⇔ u * = 1 − e 2−t and the associated state is
If the Hamiltonian is linear in the control variable u, it can be difficult to calculate u * from the optimality equation, since ∂H ∂u would not contain u. Specific ways of solving these kind of problems can be found in [81] .
Until here we have showed necessary conditions to solve basic optimal control problems. Now, it is important to study some conditions that can guarantee the existence of a finite objective functional value at the optimal control and state variables, based on [53, 75, 81, 86] . The following is an example of a sufficient condition result.
Suppose that f (t, x, u) and g(t, x, u) are both continuously differentiable functions in their three arguments and concave in x and u. Suppose u * is a control with associated state x * , and λ a piecewise differentiable function, such that u * , x * and λ together satisfy on t 0 ≤ t ≤ t f :
Then for all controls u, we have
Proof. The proof of this theorem is available on [81] .
This result is not strong enough to guarantee that J(u * ) is finite. Such results usually require some conditions on f and/or g. Next theorem is an example of an existence result from [53] . Theorem 1.3.3. Let the set of controls for problem (1.1) be Lebesgue integrable functions on t 0 ≤ t ≤ t f in R. Suppose that f (t, x, u) is convex in u, and there exist constants C 1 , C 2 , C 3 > 0, C 4 and β > 1 such that
Then there exists an optimal control u * maximizing J(u), with J(u * ) finite.
Proof. The proof of this theorem is available on [53] .
For a minimization problem, g would have a concave property and the inequality on f would be reversed. Note that the necessary conditions developed to this point deal with piecewise continuous optimal controls, while this existence theorem guarantees an optimal control which is only Lebesgue integrable. This disconnection can be overcome by extending the necessary conditions to Lebesgue integrable functions [81, 86] , but we did not expose this idea in the thesis. See the existence of OC results in [52] .
In some cases it is necessary, not only minimize (or maximize) terms over the entire time interval, but also minimize (or maximize) a function value at one particular point in time, specifically, the end of the time interval. There are some situations where the objective function must take into account the value of the state at the terminal time, e.g., the number of infected individuals at the final time in an epidemic model [81] .
Definition 5 (OC problem with payoff term). An OC problem with payoff term is in the form
where φ(x(t f )) is a goal with respect to the final position or population level x(t f ). The term φ(x(t f )) is called payoff or salvage.
Using the PMP, adapted necessary conditions can be derived for this problem.
Proposition 1 (Necessary conditions). If u * (t) and x * (t) are optimal for problem (1.8), then there exists a piecewise differentiable adjoint variable λ(t) such that
for all controls u at each time t, where H is the Hamiltonian previously defined and
Proof. The proof of this result can be found in [75] .
A new example is given to illustrate this proposition.
Example 2 (from [97] ). Let x(t) represent the number of tumor cells at time t, with exponential growth factor α, and u(t) the drug concentration. The aim is to minimize the number of tumor cells at the end of the treatment period and the accumulated harmful effects of the drug on the body. This problem is formulated as
Let us consider the Hamiltonian
The optimality condition is given by
The adjoint condition is given by
Using the transversality condition λ(t f ) = 1 (note that φ(s) = s, so φ (s) = 1), we obtain
The optimally state trajectory is (usingẋ = αx − u and x(0) = x 0 ):
Many problems require bounds on the control to achieve a realistic solution. For example, the amount of drugs in the organism must be non-negative and it is necessary to impose a limit. For the last example, despite being simplistic, makes more sense to constraint the control as 0 ≤ u ≤ 1.
Definition 6 (OC with bounded control). An OC with bounded control can be written in the form
where a, b are fixed real constants and a < b.
To solve problems with bounds on the control, it is necessary to develop alternative necessary conditions. Proposition 2 (Necessary conditions). If u * (t) and x * (t) are optimal for problem (1.9), then there exists a piecewise differentiable adjoint variable λ(t) such that
for all controls u at each time t, where H is the Hamiltonian previously defined and
λ(t f ) = 0 (transversality condition).
By an adaptation of the PMP, the OC must satisfy (optimality condition):
i.e., the maximization is over all admissible controls, andũ is obtained by the expression ∂H ∂u = 0. In particular, the optimal control u * maximizes H pointwise with respect to a ≤ u ≤ b.
Proof. The proof of this result can be found in [75] .
If we have a minimization problem, then u * is instead chosen to minimize H pointwise. This has the effect of reversing < and > in the first and third lines of optimality condition.
Remark 3. In some software packages there are no specific characterization for the bounds of the control. In those cases, and when the implementation allows, we can write in a compact way the optimal controlũ obtained without truncation, bounded by a and b: u * (t) = min{a, max{b,ũ}}.
So far, we have only examined problems with one control and with one dependent state variable. Often, it is necessary to consider more variables.
Definition 7 (OC with several variables and several controls). An OC with n state variables, m control variables and a payoff function φ can be written in the form max u1,...,um
x i (t 0 ) = x i0 , i = 1, 2, . . . , n (1.10)
where the functions f , g i are continuously differentiable in all variables.
From now on, to simplify the notation, let x(t) = [x 1 (t), . . . , x n (t)], u(t) = [u 1 (t), . . . , u m (t)], x 0 = [x 10 , . . . , x n0 ], and g(t, x, u) = [g 1 (t, x, u), . . . , g n (t, x, u)].
So, the previous problem can be rewritten in a compact way as
x(t 0 ) = x 0 , i = 1, 2, . . . , n (1.11)
Using the same approach of the previous subsections, it is possible to derive generalized necessary conditions. Proposition 3 (Necessary conditions). Let u * be a vector of optimal control functions and x * be the vector of corresponding optimal state variables. With n states, we will need n adjoints, one for each state. There is a piecewise differentiable vector-valued function λ(t) = [λ 1 (t), . . . , λ n (t)] , where each λ i is the adjoint variable corresponding to x i , and the Hamiltonian is
It is possible to find the variables satisfying identical optimality, adjoint and transversality conditions in each vector component. Namely, u * maximizes H(t, x * , u, λ) with respect to u at each t, and u * , x * and λ satisfy
x(t f )) for j = 1, . . . , n (transversality conditions) ∂H ∂u k = 0 at u * k for k = 1, . . . , m (optimality conditions) By φ xj , it is meant the partial derivative in the x j component. Note, if φ ≡ 0, then λ j (t f ) = 0 for all j, as usual.
Similarly to the previous section, if bounds are placed on a control variable, a k ≤ u k ≤ b k (for k = 1, . . . , m), then the optimality condition is changed from ∂H ∂u k = 0 to
Below, an optimal control problem related to rubella is presented.
Rubella, commonly known as German measles, is a most common in child age, caused by the rubella virus. Children recover more quickly than adults, and can be very serious in pregnancy. The virus is contracted through the respiratory tract and has an incubation period of 2 to 3 weeks. The primary symptom of rubella virus infection is the appearance of a rash on the face which spreads to the trunk and limbs and usually fades after three days. Other symptoms include low grade fever, swollen glands, joint pains, headache and conjunctivitis.
It is presented now an optimal control problem to study the dynamics of rubella in China over three years, using a vaccination process (u) as a measure to control the disease (more details can be found in [17] ). Let x 1 represent the susceptible population, x 2 the proportion of population that is in the incubation period, x 3 the proportion of population that is infected with rubella and x 4 the rule that remains the population constant. The optimal control problem can be defined as: It is very difficult to solve analytically this problem. For most of the epidemiologic problems it is necessary to employ numerical methods. Some of them will be described in the next chapter.
In this chapter, some numerical approaches to solve a system of ordinary differential equations, such as shooting methods and multi-steps methods, are introduced. Then, two distinct philosophies to solve OC problems are presented: indirect methods centered in the PMP and the direct ones focus on problem discretization solved by nonlinear optimization codes. A set of software packages used all over the thesis is summarily exposed.
In the last decades the computational world has been developed in an amazing way. Not only in hardware issues such as efficiency, memory capacity, speed, but also in terms of the software robustness. Groundbreaking achievements in the field of numerical solution techniques for differential and integral equations have enabled the simulation of highly complex real world scenarios. This way, OC also won with these improvements and numerical methods and algorithms have evolved significantly.
The next section concerns on the resolution of differential equations systems.
A dynamic system is mathematically characterized by a set of ordinary differential equations (ODE). Specifically, the dynamics are described for t 0 ≤ t ≤ t f , by a system of n ODEṡ
. . , y n (t), t) f 2 (y 1 (t), . . . , y n (t), t)
. . .
The problems of solving an ODE are classified into initial value problems (IVP) and boundary value problems (BVP), depending on how the conditions at the endpoints of the domain are specified. All the conditions of an initial-value problem are specified at the initial point. On the other hand, the problem becomes a boundary-value problem if the conditions are needed for both initial and final points.
There exist many numerical methods to solve initial value problems -such as Euler, Runge-Kutta or adaptive methods -and boundary value problems, such as shooting methods.
One can visualize the shooting method as the simplest technique for solving BVP. Supposing it is desired to determine the initial angle of a cannon so that, when a cannonball is fired, it strikes a desired target. An initial guess is made for the angle of the cannon, and the cannon is fired. If the cannon does not hit the target, the angle is adjusted based on the amount of the miss and another cannon is fired. The process is repeated until the target is hit [108] .
Suppose we want to find y(t 0 ) = y 0 such that y(t f ) = b. The shooting method can be summarized as follows [9] :
Step 1. guess initial conditions x = y(t 0 );
Step 2. propagate the differential system equations from t 0 to t f , i.e., shoot;
Step 3. evaluate the error in the boundary conditions c(x) = y(t f ) − b;
Step 4. use a nonlinear program to adjust the variables x to satisfy the constraints c(x) = 0, i.e., repeat steps 1-3. Despite its simplicity, from a practical standpoint, the shooting method is used only when the problem has a small number of variables. This method has a major disadvantage: a small change in the initial condition can produce a very large change in the final conditions. In order to overcome the numerical difficulties of the simple method, the multiple shooting method is presented.
In a multiple shooting method, the time-interval [t 0 , t f ] is divided into M − 1 subintervals. Then is applied over each subinterval [t i , t i+1 ] with the initial values of the differential equations in the interior intervals being unknown that need to be determined. In order to enforce continuity, the following conditions are enforced at the interface of each subinterval:
A scheme of the multiple shooting method is shown in Figure 2 With the multiple shooting approach the problem size is increased: additional variables and constraints are introduced for each shooting segment. In particular, the number of nonlinear variables and constraints for a multiple shooting application is n = n y (M − 1), where n y is the number of dynamic variables y and M − 1 is the number of segments [9] .
Both shooting and multiple shooting methods require a good guess for initial conditions and the propagation of the shoots for problem of high dimension is not feasible. For this reason, other methods can be implemented, using initial value problems.
The numerical solution of the IVP is fundamental to most optimal control methods. The problem can stand as follows: compute the value of y(t f ) for some value of t < t f that satisfies (2.1) with the known initial value y(t 0 ) = y 0 .
Numerical methods for solving the ODE IVP are relatively mature in comparison to other fields in Optimal Control. It will be considered two methods: single-step and multiple-step methods. In both, the solution of the differential system at each step t k is sequentially obtained using current and/or previous information about the solution. In both cases, it is assumed that the time t = nh moves ahead in uniform steps of length h [9, 51] .
The most common single-step method is Euler method. In this discretization scheme, if a differential equation is written likeẋ = f (x(t), t), is possible to make a convenient approximation of this:
This approximation x n+1 of x(t) at the point t n+1 has an error of order h 2 . Clearly, there is a trade-off between accuracy and complexity of calculation which depends heavily on the chosen value for h. In general, as h is decreased the calculation takes longer but is more accurate.
For many higher order systems it is very difficult to make Euler approximation effective. For this reason more accurate and elaborate techniques were developed. One of these methods is the Runge-Kutta method.
A Runge-Kutta method is a multiple-step method, where the solution at time t k+1 is obtained from a defined set of previous values t j−k , . . . , t k and j is the number of steps.
If a differential equation is written likeẋ = f (x(t), t), it is possible to make a convenient approximation of this, using the second order Runge-Kutta method
or the fourth order Runge-Kutta method
x n+1 x n + h 6 (k 1 + 2k 2 + 2k 3 + k4)
where
This approximation x n+1 of x(t) at the point t n+1 has an error depending on h 3 and h 5 , for the Runge-Kutta methods of second and fourth order, respectively.
Numerical methods for solving OC problems date back to the 1950s with Bellman investigation. From that time to present, the complexity of methods and corresponding complexity and variety of applications has substantially increased [108] .
There are two major classes of numerical methods for solving OC problems: indirect methods and direct methods. The first ones, indirectly solve the problem by converting the optimal control problem to a boundary-value problem, using the PMP. On the other hand, in a direct method, the optimal solution is found by transcribing an infinite-dimensional optimization problem to a finite-dimensional optimization problem.
In an indirect method, the PMP is used to determine the first-order optimality conditions of the original OC problem. The indirect approach leads to a multiple-point boundary-value problem that is solved to determine candidate optimal trajectories called extremals.
For an indirect method it is necessary to explicitly get the adjoint equations, the control equations and all the transversality conditions, if there exist. Notice that there is no correlation between the method used to solve the problem and its formulation: one may consider applying a multiple shooting method solution technique to either an indirect or a direct formulation. In the following subsection a numerical approach using the indirect method is presented.
This method is described in a recent book by Lenhart and Workman [81] and it is known as forwardbackward sweep method. The process begins with an initial guess on the control variable. Then, the state equations are simultaneously solved forward in time and the adjoint equations are solved backward in time. The control is updated by inserting the new values of states and adjoints into its characterization, and the process is repeated until convergence occurs.
Considering x = (x 1 , . . . , x N + 1) and λ = (λ 1 , . . . , λ N + 1) the vector approximations for the state and the adjoint. The main idea of the algorithm is described as follows:
Step 1. Make an initial guess for u over the interval ( u ≡ 0 is almost always sufficient);
Step 2. Using the initial condition x 1 = x(t 0 ) = a and the values for u, solve x forward in time according to its differential equation in the optimality system;
Step 3. Using the transversality condition λ N +1 = λ(t f ) = 0 and the values for u and x, solve λ backward in time according to its differential equation in the optimality system;
Step 4. Update u by entering the new x and λ values into the characterization of the optimal control;
Step 5. Verify convergence: if the variables are sufficiently close to the corresponding in the previous iteration, then output the current values as solutions, else return to Step 2.
For Steps 2 and 3, Lenhart and Workman used for the state and adjoint systems the Runge-Kutta fourth order procedure to make the discretization process.
On the other hand, Wang [136] , applied the same philosophy but solving the differential equations with the solver ode45 for Matlab. This solver is based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair. That means the numerical solver ode45 combines a fourth and a fifth order methods, both of which are similar to the classical fourth order Runge-Kutta method discussed above. These vary the step size, choosing it at each step an attempt to achieve the desired accuracy. Therefore, the solver ode45 is suitable for a wide variety of initial value problems in practical applications. In general, ode45 is the best method to apply as a first attempt for most problems [70] .
Let consider the open problem defined in Chapter 1 (Example 3) about rubella disease. With x(t) = (x 1 (t), x 2 (t), x 3 (t), x 4 (t)) and λ(t) = (λ 1 (t), λ 2 (t), λ 3 (t), λ 4 (t)), the Hamiltonian of this problem can be written as
Using the PMP the optimal control problem can be studied with the state variableṡ
with initial conditions x 1 (0) = 0.0555, x 2 (0) = 0.0003, x 3 (0) = 0.0004 and x 4 (0) = 1 and the adjoint variables:
The optimal control is
Here it is only presented the main part of the code using the backward-forward sweep method with fourth order Runge-Kutta. The completed one can be found in the website [110] .
.. beta*(x1(i)+h2*m11)*(x3(i)+h2*m13)-(0.5*(u(i) + u(i+1)))*(x1(i)+h2*m11); m22 = b*p*(x2(i)+h2*m12)+beta*(x1(i)+h2*m11)*(x3(i)+h2*m13)-(e+b)*(x2(i)+h2*m12); m23 = e*(x2(i)+h2*m12)-(g+b)*(x3(i)+h2*m13); m24 = b-b*(x4(i)+h2*m14); m31 = b-b*(p*(x2(i)+h2*m22)+q*(x3(i)+h2*m23))-b*(x1(i)+h2*m21)-... beta*(x1(i)+h2*m21)*(x3(i)+h2*m23)-(0.5*(u(i) + u(i+1)))*(x1(i)+h2*m21); m32 = b*p*(x2(i)+h2*m22)+beta*(x1(i)+h2*m21)*(x3(i)+h2*m23)-(e+b)*(x2(i)+h2*m22); m33 = e*(x2(i)+h2*m22)-(g+b)*(x3(i)+h2*m23); m34 = b-b*(x4(i)+h2*m24); m41 = b-b*(p*(x2(i)+h2*m32)+q*(x3(i)+h2*m33))-b*(x1(i)+h2*m31)-... beta*(x1(i)+h2*m31)*(x3(i)+h2*m33)-u(i+1)*(x1(i)+h2*m31); m42 = b*p*(x2(i)+h2*m32)+beta*(x1(i)+h2*m31)*(x3(i)+h2*m33)-(e+b)*(x2(i)+h2*m32); m43 = e*(x2(i)+h2*m32)-(g+b)*(x3(i)+h2*m33); m44 = b-b*(x4(i)+h2*m34); x1(i+1) = x1(i) + (h/6)*(m11 + 2*m21 + 2*m31 + m41); x2(i+1) = x2(i) + (h/6)*(m12 + 2*m22 + 2*m32 + m42); x3(i+1) = x3(i) + (h/6)*(m13 + 2*m23 + 2*m33 + m43); x4(i+1) = x4(i) + (h/6)*(m14 + 2*m24 + 2*m34 + m44); end for i = 1:M j = M + 2 -i; n11 = lambda1(j)*(b+u(j)+beta*x3(j))-lambda2(j)*beta*x3(j); n12 = lambda1(j)*b*p+lambda2(j)*(e+b-p*b)-lambda3(j)*e; n13 = -A+lambda1(j)*(b*q+beta*x1(j))-lambda2(j)*beta*x1(j)+lambda3(j)*(g+b); n14 = b*lambda4(j); n21 = (lambda1(j) -h2*n11)*(b+u(j)+beta*(0.5*(x3(j)+x3(j-1))))-... (lambda2(j) -h2*n12)*beta*(0.5*(x3(j)+x3(j-1))); n22 = (lambda1(j) -h2*n11)*b*p+(lambda2(j) -h2*n12)*(e+b-p*b)-(lambda3(j) -h2*n13)*e; n23 = -A+(lambda1(j) -h2*n11)*(b*q+beta*(0.5*(x1(j)+x1(j-1))))-...
(lambda2(j) -h2*n12)*beta*(0.5*(x1(j)+x1(j-1)))+(lambda3(j) -h2*n13)*(g+b); n24 = b*(lambda4(j) -h2*n14); n31 = (lambda1(j) -h2*n21)*(b+u(j)+beta*(0.5*(x3(j)+x3(j-1))))-...
(lambda2(j) -h2*n22)*beta*(0.5*(x3(j)+x3(j-1))); n32 = (lambda1(j) -h2*n21)*b*p+(lambda2(j) -h2*n22)*(e+b-p*b)-(lambda3(j) -h2*n23)*e; n33 = -A+(lambda1(j) -h2*n21)*(b*q+beta*(0.5*(x1(j)+x1(j-1))))-...
(lambda2(j) -h2*n22)*beta*(0.5*(x1(j)+x1(j-1)))+(lambda3(j) -h2*n23)*(g+b); n34 = b*(lambda4(j) -h2*n24); n41 = (lambda1(j) -h2*n31)*(b+u(j)+beta*x3(j-1))-(lambda2(j) -h2*n32)*beta*x3(j-1); n42 = (lambda1(j) -h2*n31)*b*p+(lambda2(j) -h2*n32)*(e+b-p*b)-(lambda3(j) -h2*n33)*e; n43 = -A+(lambda1(j) -h2*n31)*(b*q+beta*x1(j-1))-...
(lambda2(j) -h2*n32)*beta*x1(j-1)+(lambda3(j) -h2*n33)*(g+b); n44 = b*(lambda4(j) -h2*n34); lambda1(j-1) = lambda1(j) -h/6*(n11 + 2*n21 + 2*n31 + n41); lambda2(j-1) = lambda2(j) -h/6*(n12 + 2*n22 + 2*n32 + n42); lambda3(j-1) = lambda3(j) -h/6*(n13 + 2*n23 + 2*n33 + n43); lambda4(j-1) = lambda4(j) -h/6*(n14 + 2*n24 + 2*n34 + n44); end u1 = min(0.9,max(0,lambda1.*x1/2));
The optimal curves for the states variables and optimal control are shown in Figure 2 There are several difficulties to overcome when an optimal control problem is solved by indirect methods. Firstly, is necessary to calculate the hamiltonian, adjoint equations, optimality condition and transversality conditions. Besides, the approach is not flexible, since each time a new problem is formulated, a new derivation is required. In contrast, a direct method does not require explicit derivation neither the necessary conditions.
Due to these practical difficulties with the indirect formulation, the main focus will be centered on the direct methods. This approach has been gaining popularity in numerical optimal control over the past three decades [9] .
A new family of numerical methods for dynamic optimization has emerged, referred to as direct methods. This development has been driven by the industrial need to solve large-scale optimization problems and it has also been supported by the rapidly increasing computational power.
A direct method constructs a sequence of points x 1 , x 2 , . . . , x * such that the objective function in minimized and typical F (x 1 ) > F (x 2 ) > · · · > F (x * ). Here the state and/or control are approximated using an appropriate function approximation (e.g., polynomial approximation or piecewise constant parameterization). Simultaneously, the cost functional is approximated as a cost function. Then, the coefficients of the function approximations are treated as optimization variables and the problem is reformulated to a standard nonlinear optimization problem (NLP) in the form:
where c i , i ∈ E e c j , j ∈ I are the set of equality and inequality constraints, respectively. In fact, the NLP is easier to solve than the boundary-value problem, mainly due to the sparsity of the NLP and the many well-known software programs that can handle with this feature. As a result, the range of problems that can be solved via direct methods is significantly larger than the range of problems that can be solved via indirect methods. Direct methods have become so popular these days that many people have written sophisticated software programs that employ these methods. Here we present two types of codes/packages: specific solvers for OC problems and standard NLP solvers used after a discretization process.
The OC-ODE [57] , Optimal Control of Ordinary-Differential Equations, by Matthias Gerdts, is a collection of Fortran 77 routines for optimal control problems subject to ordinary differential equations. It uses an automatic direct discretization method for the transformation of the OC problem into a finite-dimensional NLP. OC-ODE includes procedures for numerical adjoint estimation and sensitivity analysis.
Considering the same problem (Example 3), here is the main part of the code in OC-ODE. The completed one can be found in the website [110] . The achieved solution is similar to the indirect approach, therefore we will not present it. 
The DOTcvp [67] , Dynamic Optimization Toolbox with Vector Control Parametrization is a dynamic optimization toolbox for Matlab. The toolbox provides environment for a Fortran compiler to create the '.dll' files of the ODE, Jacobian, and sensitivities. However, a Fortran compiler has to be installed in a Matlab environment. The toolbox uses the control vector parametrization approach for the calculation of the optimal control profiles, giving a piecewise solution for the control. The OC problem has to be defined in Mayer form. For solving the NLP, the user can choose several deterministic solvers -Ipopt, Fmincon, FSQP -or stochastic solvers -DE, SRES.
The modified SUNDIALS tool [66] is used for solving the IVP and for the gradients and Jacobian automatic generation. Forward integration of the ODE system is ensured by CVODES, a part of SUNDIALS, which is able to perform the simultaneous or staggered sensitivity analysis too. The IVP problem can be solved with the Newton or Functional iteration module and with the Adams or BDF linear multistep method. Note that the sensitivity equations are analytically provided and the error control strategy for the sensitivity variables could be enabled. DOTcvp has a user friendly graphical interface (GUI).
Considering the same problem (Example 3), here is a part of the code used in DOTcvp. The completed one can be found in the website [110] . The solution, despite being piecewise continuous, follows the curve obtained by the previous programs. Functional for non-stiff problems data.odes.LinearSolver With the GUI interface this method was the preferred to be tested, due to its simple way to implement the code.
In NEOS [98] platform there is a large set of software packages. NEOS is considered as the state of the art in optimization. One recent solver is Muscod-II [79] (Multiple Shooting CODe for Optimal Control) for the solution of mixed integer nonlinear ODE or DAE constrained optimal control problems in an extended AMPL format.
AMPL [55] is a modelling language for mathematical programming and was created by Fourer, Gay and Kernighan. The modelling languages organize and automate the tasks of modelling, which can handle a large volume of data and, moreover, can be used in machines and independent solvers, allowing the user to concentrate on the model instead of the methodology to reach solution. However, the AMPL modelling language itself does not allow the formulation of differential equations. Hence, the TACO Toolkit has been designed to implement a small set of extensions for easy and convenient modeling of optimal control problems in AMPL, without the need for explicit encoding of discretization schemes. Both the TACO Toolkit and the NEOS interface to Muscod-II are still under development.
Probably for this reason, the Example 3 could not be solved by this software which crashed after some runs. Anyway, we opted to also put the code, for the same example that is being used, to show the differences of modelling language used in each program.
include OptimalControl.mod; var t ; var x1, >=0 <=1; var x2, >=0 <=1; var x3, >=0 <=1; var x4, >=0 <=1; var u >=0, <=0.9 suffix type "u0";
The three nonlinear optimization software packages presented here, were used through the NEOS platform with codes formulated in AMPL.
The Ipopt [135] , Interior Point OPTimizer, is a software package for large-scale nonlinear optimization. It is written in Fortran and C. Ipopt implements a primal-dual interior point method and uses a line search strategy based on filter method. Ipopt can be used from various modeling environments.
Ipopt is designed to exploit 1st and 2nd derivative information if provided, usually via automatic differentiation routines in modeling environments such as AMPL. If no Hessians are provided, Ipopt will approximate them using a quasi-Newton methods, specifically a BFGS update.
Continuing with Example 3 the AMPL code is shown for Ipopt. The Euler discretization was the selected. Indeed, this code can also be implemented in other nonlinear software packages available in NEOS platform, reason why the code for the next two software packages will not be shown. The full version can be found on the website [110] .
, short for "Nonlinear Interior point Trust Region Optimization", was created primarily by Richard Waltz, Jorge Nocedal, Todd Plantenga and Richard Byrd. It was introduced in 2001 as a derivative of academic research at Northwestern, and has undergone continual improvement since then.
Knitro is also a software for solving large scale mathematical optimization problems based mainly on the two Interior Point (IP) methods and one active set algorithm. Knitro is specialized for nonlinear optimization, but also solves linear programming problems, quadratic programming problems, and systems of nonlinear equations. The unknowns in these problems must be continuous variables in continuous functions; however, functions can be convex or nonconvex. The code also provides a multistart option for promoting the computation of the global minimum. This software was tested through the NEOS platform.
Snopt [58] , by Philip Gill, Walter Murray and Michael Saunders, is a software package for solving large-scale optimization problems (linear and nonlinear programs). It is specially effective for nonlinear problems whose functions and gradients are expensive to evaluate. The functions should be smooth but do not need to be convex. Snopt is implemented in Fortran 77 and distributed as source code. It uses the SQP (Sequential Quadratic Programming) philosophy, with an augmented Lagrangian approach combining a trust region approach adapted to handle the bound constraints. Snopt is also available in NEOS platform.
Choosing a method for solving an OC problem depends largely on the type of problem to be solved and the amount of time that can be invested in coding. An indirect shooting method has the advantage of being simple to understand and produces highly accurate solutions when it converges [108] . The accuracy and robustness of a direct method is highly dependent upon the method used. Nevertheless, it is easier to formulate highly complex problems in a direct way and can be used standard NLP solvers as an extra advantage. This last feature has the benefit of converging with poor initial guesses and being extremely computationally efficient since most of the solvers exploit the sparsity of the derivatives in the constraints and objective function.
In the next chapter the basic concepts from epidemiology are provided, in order to formulate/implement more complex OC problems in the health area.
In this chapter, the simplest epidemiologic models composed by mutually exclusive compartments are introduced. Based on the models SIS (susceptible-infected-susceptible) and SIR (susceptible-infected-recovered), other models are presented introducing new issues related to maternal immunity or the latent period, fitting the features to distinct diseases. Illustrative examples are presented, with diseases that can be described by each model. The basic reproduction number is calculated and presented as a threshold value for the eradication or the persistence of the disease in a population.
In the 14th century, occurred one of most famous epidemic events: the Black Death. It killed approximately one third of the European population. From 1918-19, twenty to forty percent of the world's population suffered from the Spanish Flu, the most severe pandemic in history. In 1978, the United Nations promoted an ambitious agreement between the countries forecasting that in the year 2000 infectious diseases would be eradicated. This conjecture failed, mainly due to the assumption that the microorganisms were biologically stationary and consequently they were not modified and became resistant to the medicines. Besides, the improvements in the transportation allowing for a faster movement of individuals and the population growing especially in developing countries, led to the appearance of new diseases and the resurgence of old ones in distinct places. Nowadays, AIDS is the most scrutinized. In 2007 there were an estimated 33.2 million sufferers worldwide, and 2.1 million deaths with over three quarters of these occurring in sub-Saharan Africa [93] .
Epidemiology -the study of patterns of diseases including those which are non-communicable of infections in population -has become more relevant and indispensable in the development of new models and explanations for the outbreaks, namely due to their propagation and causes. In epidemiology, an infection is said to be endemic in a population when it is maintained in the population without the need for external inputs. An epidemic occurs when new cases of a certain disease appears, in a given human population during a given period, and then essentially disappears.
There are several types of diseases, depending on their type of transmission mechanism, of which stand out:
• bacteria, which do not confer immunity against the reinfection and frequently produce harmful toxins to the host; in case of infection, the antibiotics are usually efficient (examples: tuberculosis, meningitis, gonorrhea, syphilis, tetanus);
• viral agents, that confer immunity against reinfection; here the antibiotics do not produce effects and usually it is hoped that the immune system of the host responds to an infection by the virus or it will be necessary to take antiviral drugs that retard the multiplication of the virus (examples influenza, chicken pox, measles, rubella, mumps, HIV / AIDS, smallpox);
• vectors, that are usually mosquitoes or ticks and are infected by humans and then transmit the disease to other humans (examples: malaria, yellow fever, dengue, chikungunya).
The transmission can happen in a direct or indirect way. The direct transmission of a disease can happen by physical proximity (such as sneezing, coughing, kissing, sexual contact) or even by a specific parasite that penetrates the host through ingestion or the skin. The indirect transmission involves the vectors that are intermediaries or carriers of the infection.
In most of the cases, the direct and indirect transmission of the disease happens between the member that coexists in the host population; this is called the horizontal transmission. When the direct transmission occurs from one ascendent to a descendent not yet born (egg or embryo) it is said that vertical transmission happens [76] . When formulating a model for a particular disease, we should make a trade-off between simple models -that omit several details and generally are used for specific situations in a short time, but have the disadvantage of possibly being naive and unrealistic -and more complex models, with more details and more realistic, but generally more difficult to solve or could contain parameters which their estimates cannot be obtained.
Choosing the most appropriated model depends on the precision or generality required, the available data, and the time frame in which the results are needed. By definition, all models are "wrong", in the sense that even the most complex will make some simplifying assumptions. It is, therefore, difficult to definitively express which model is right, though naturally we are interested in developing models that capture the essential features of a system. The art of epidemiological modelling is to make suitable choices in the model formulation making it as simple as possible and yet suitable for the question being considered [65] .
Mathematical models are a simplified representation of how an infection spreads across a population over time, and generally come in two forms: stochastic and deterministic models. The first ones, employ randomness, with variables being described by probability distributions. Deterministic models split the population into subclasses, and an ODE with respect to time is formulated for each. The state variable are determined using parameters and initial conditions. The main focus in this chapter will be the deterministic models, neglecting the others.
Most epidemic models are based on dividing the population into a small number of compartments. Each containing individuals that are identical in terms of their status with respect to the disease in question.
Here are some of the main compartments that a model can contain.
• Passive immune (M ): is composed by newborns that are temporary passively immune due to antibodies transferred by their mothers;
• Susceptible (S): is the class of individuals who are susceptible to infection; this can include the passively immune 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 referred to the individuals that despite infected, do not exhibit obvious signs of infection and the abundance of the pathogen may be too low to allow further transmission;
• Infected (I): in this class, the level of parasite is sufficiently large within the host and there is potential in transmitting the infection to other susceptible individuals;
• Recovered or Resistant (R): includes all individuals who have been infected and have recovered.
The choice of which compartments to include in a model depends on the characteristics of the particular disease being modelled and the purpose of the model. The exposed compartment is sometimes neglected, when the latent period is considered very short. Besides, the compartment of the recovered individuals cannot always be considered since there are diseases where the host has never became resistent. Acronyms for epidemiology models are often based on the flow patterns between the compartments such as MSEIR, MSEIRS, SEIR, SEIRS, SIR, SIRS, SEI, SEIS, SI, SIS.
There are three commonly used threshold values in epidemiology: R 0 , σ and R. The most common and probably the most important is the basic reproduction number [61, 64, 65] .
Definition 8 (Basic reproduction number). The basic reproduction number, denoted by R 0 , is defined as the average number of secondary infections that occurs when one infective is introduced into a completely susceptible population.
This threshold, R 0 , is a famous result due to Kermack and McKendrick [77] and is referred to as the "threshold phenomenon", giving a borderline between a persistence or a disease death. R 0 it is also called the basic reproduction ratio or basic reproductive rate.
Definition 9 (Contact number). The contact number, σ is the average number of adequate contacts of a typical infective during the infectious period.
An adequate contact is one that is sufficient for transmission, if the individual contacted by the susceptible is an infective. It is implicitly assumed that the infected outsider is in the host population for the entire infectious period and mixes with the host population in exactly the same way that a population native would mix.
Definition 10 (Replacement number). The replacement number, R, is the average number of secondary infections produced by a typical infective during the entire period of infectiousness.
Note that the replacement number R changes as a function of time t as the disease evolves after the initial invasion.
These three quantities R 0 , σ and R are all equal at the beginning of the spreading of an infectious disease when the entire population (except the infective invader) is susceptible. R 0 is only defined at the time of invasion, whereas σ and R are defined at all times.
The replacement number R is the actual number of secondary cases from a typical infective, so that after the infection has invaded a population and everyone is no longer susceptible, R is always less than the basic reproduction number R 0 . Also after the invasion, the susceptible fraction is less than one, and as such not all adequate contacts result in a new case. Thus the replacement number R is always less than the contact number σ after the invasion [64] . Combining these results leads to
Note that R 0 = σ for most models, and σ > R after the invasion for all models.
For the models throughout this study the basic reproduction number, R 0 , will be applied. When R 0 < 1 the disease cannot invade the population and the infection will die out over a period of time. The amount of time this will take generally depends on how small R 0 is. When
invasion is possible and infection can spread through the population. Generally, the larger the value of R 0 the more severe, and possibly widespread, the epidemic will be [42] .
In Table 3 .1 are some example diseases with their estimated basic reproduction number. Due to differences in demographic rates, rural-urban gradients, and contact structure, different human populations may be associated with different values of R 0 for the same disease [7] . In the next section some of the epidemiologic models will be presented.
Numerous infectious diseases confer no long-lasting immunity. The SIS models are suitable for some bacterial agent diseases like meningitis, sexually transmitted diseases such as gonorrhea and for protozoan agent diseases where malaria and the sleeping sickness are good examples. For these diseases, an individual can be infected multiple times throughout their lives, with no apparent immunity. Here, recovery from infection is followed by an instant return to the susceptible compartment.
Throughout this chapter we will consider that population is constant, neglecting the tourism and immigration factors. Also it is considered that the population is homogeneously mixed, which means that every individual interacts with another at the same level and therefore all individuals have the same risk of contracting the disease.
The number of individuals in each compartment must be integer, but if the population size N is sufficiently large, it is possible to treat S and I as continuous variables. Calculating the proportion of these compartments, varying from 0 to 1, it is considered that the total population is constant over time, i.e., 1 = S + I. The compartment changes are expressed by a system of differential equations.
The SIS model can be mathematically represented as follows.
Definition 11 (SIS model). The SIS model can be formulated as
subject to initial conditions S(0) > 0 and I(0) ≥ 0.
It is considered β the transmission rate (per capita) and γ the recovery rate, so the mean infectious period is 1/γ. The vital dynamics (births and deaths) were not considered, but a similar model can be constructed with these effects [65] . Despite this lack of susceptible births, the disease can still persist because the recovery of infected individuals replenishes the susceptible class and will guarantee the long-term persistence of the disease.
Remark 5. The SI model is a particular case of the SIS model when the recovery rate (γ) is null.
A new infected individual is expected to infect another at a transmission rate β, during 1 γ that is the infectious period, so the expected basic reproduction number is
Trachoma is an infectious disease causing a characteristic roughening of the inner surface of the eyelids. Also called granular conjunctivitis or Egyptian ophthalmia, it is the leading cause of infectious blindness in the world. Adapting a model from [109] , and using β = 0.047 as transmission rate and the recovery rate γ = 0.017, we have the basic reproduction number R 0 approximately 2.76 and the following representation of the state variables (Figure 3 .2). 
The SIR model was initially studied in depth by Kermack and McKendrick and categorizes hosts within a population as Susceptible, Infected and Recovered [77] . It captures the dynamics of acute infections that confers lifelong immunity once recovered. Diseases where individuals acquire permanent immunity, and for which this model may be applied, include measles, smallpox, chickenpox, mumps, typhoid fever and diphtheria.
Once again we will consider the total population size constant, i.e., 1 = S + I + R. Two cases will be studied, distinguished by the inclusion or exclusion of demographic factors.
Having compartmentalized the population, we now need a set of equations that specify how the sizes of compartments change over time.
Definition 12 (SIR model without demography). The SIR model, excluding births and deaths, can be defined as
In addition, the transmission rate, per capita, is β and the recovery rate is γ. Since population is constant and as R does not appear in the first two differential equations, most of the times the last equation is omitted, indeed, R(t) = 1 − S(t) − I(t). By assuming equations from (3.2) is possible to notice that dS dt
Then a newly introduced infected individual can be expected to infect other people at the rate β during the expected infectious period 1/γ. Thus, this first infective individual can be expected to infect
Consider an epidemic of influenza in a British boarding school in early 1978 [76] . Three boys were reported to the school infirmary with the typical symptoms of influenza. Over the next few days, a very large fraction of the 763 boys in the school had contact with the infection. Within two weeks, the infection had become extinguished. The best fit parameters yield an estimated active infectious period of 1/γ = 2.2 days and a mean transmission rate β = 1.66 per day. Therefore, the estimated R 0 is 3.652. Figure 3 .4 represents the dynamics of the three state variables.
It was presented a SIR model given the assumption that the time scale of disease spread was sufficiently fast that births and deaths can be neglected. In the next subsection we explore the long-term persistence and endemic dynamics of an infectious disease. 
The simplest and most common way of introducing demography into the SIR model is to assume there is a natural host lifespan, 1/µ years. Then, the rate at which individuals, at any epidemiological compartment, suffer natural mortality is given by µ. It is important to emphasize that this factor is independent of the disease and is not intended to reflect the pathogenicity of the infectious agent. Historically, it has been assumed that µ also represents the population's crude birth rate, thus ensuring that total population size does not change through time, or in other words, dS dt + dI dt + dR dt = 0. Putting all these assumptions together, we have a new definition.
Definition 13 (SIR model with demography). The SIR model, including births and deaths, can be defined as
The epidemiological scheme is in Figure 3 .5.
It is important to introduce the R 0 expression for this model. The parameter β represents the trans-mission rate per infective and the negative terms in the equation tell us that each individual spends an average 1 γ+µ time units in this class. Therefore, if we assume the entire population is susceptible, then the average number of new infectious per infectious individual is determined by
The inclusion of demographic dynamics may allow a disease to die out or persist in a population in the long term. For this reason it is important to explore what happens when the system is at equilibrium.
Definition 14 (Equilibrium points). A model defined SIR has an equilibrium point, if a triple E * = (S * , I * , R * ) satisfies the following system:
If the equilibrium point has the infectious component equal to zero (I * = 0), this means that the pathogen suffered extinction and E * is called Disease Free equilibrium (DFE).
If I * > 0 the disease persist in the population and E * is called Endemic Equilibrium (EE).
With some calculations and algebraic manipulations, it is possible to obtain two equilibria for the system (3.3):
DFE:
When R 0 < 1, each infected individual produces, on average, less than one new infected individual, and therefore, predictable that the infection will be cleared from the population. If R 0 > 1, the pathogen is able to invade the susceptible population [61, 64] . It is possible to prove that for the Endemic Equilibrium to be stable, R 0 must be greater than one, otherwise the Disease Free Equilibrium is stable. More detailed information about local and global stability of the equilibrium point can be found in [22, 74, 83, 91] .
This threshold behavior is very useful, once we can determine which control measures, and at what magnitude, would be most effective in reducing R 0 below one, providing important guidance for public health initiatives.
The SIR model below, Figure 3 .6, is plotted using similar parameters and initial conditions, except the transmission rate β (adapted from [76] ).
It is shown that in case (a), using β = 520, the basic reproduction number is greater than one. With demographic effects could have damped oscillations with a decreasing amplitude, in the EE direction. In case (b), with β = 10, we obtain R 0 < 1 and the system tends to go to a DFE.
The next three sections present a brief description of other possible refinements of the basic models SIS and SIR. 
In the SEIR case, the transmission process occurs with an initial inoculation with a very small number of pathogens. Then, during a period of time the pathogen reproduces rapidly within the host, relatively unchallenged by the immune system. During this stage, pathogen abundance is too low for active transmission to other susceptible host, and yet the pathogen is present. The time in this stage is very difficult to quantify, since there is no symptomatic features of the disease. It is called the latent or exposed period. Assuming the average duration of the latent period is given by 1 ν , the SEIR model can be described as follow.
Definition 15 (SEIR model). The SEIR model is formulated as
The parameter β and γ were defined in the previous section. The epidemiological scheme for SEIR model is presented in Figure 3 
This threshold is the product of the contact rate β per unit time, the average infectious period adjusts to the population growth of 1 γ+µ , and the fraction ν ν+µ of exposed people surviving the latent class E. Finding steady states of the system, we obtain the following equilibrium points:
Although the SIR and SEIR models behave similarly at equilibrium, when the parameters are suitably adapted, the SEIR model has a slower growth rate after pathogen invasion. This is due to the fact that individuals need to stay some time in the exposed class before their contribution in the disease transmission chain [76] .
An infected or vaccinated mother transfers some antibodies across the placenta to her fetus, so that the newborn infant has temporary passive immunity to an infection. Since the infant can not produce new antibodies, when these passive antibodies are gone at a rate δ, the baby passes from the immune state M to the susceptible state S. Some childhood diseases have this feature. The birth rate µS into the susceptible class of size S corresponds to newborns whose mothers are susceptible, and the other newborns µ(1 − S) enter passively immune class of size M , since their mothers were infected or had some type of immunity [65] .
The transfer diagram for the MSEIR model is shown in The MSEIR is also composed by a system of differential equations.
Definition 16 (MSEIR model). The SEIR model can be described as
Thus, the basic reproduction number is equal to previous model SEIR, once the M compartment does not affect the transmission chain of the disease:
The equations (3.5) always have a DFE given by
There are several models with different epidemiological states, such as MSEIRS, SEIRS, SEI, SI, SIRS, depending on the specific features and the level of detail that is wished to introduce in the model. These models are similar to the previous ones presented. Other epidemiological models can be studied using more compartments such as Quarantine-Isolation (Q), Treatment (T ), Carrier (C) or Vaccination (V ). Besides, most of the populations can be subdivided into different groups (sex, age, health weaknesses, ...), depending upon characteristics that may influence the risk of catching and/or transmitting an infection. More details and examples can be found in [4, 14, 65, 134] . Other diseases can be caught and transmitted by numerous hosts or even is required a second different population to complete the transmission cycle, such as vector borne-diseases, using double models. This last case will be explored in the second part of the thesis.
In cases that include multiple compartments of infected individuals, in which vital and epidemiological parameters depend on factors as stage of the disease, spatial position, age, behavior, multigroups, the next generation method is the more generalized approach to calculate R 0 .
The definition of R 0 has more than one possible interpretation, depending on the field (ecology, demography or epidemiology) and there exist distinct methods and estimations to calculate this threshold.
The next generation method, introduced by Diekmann et al. [39] , defines R 0 as the spectral radius of the next generator operator. The formation of the operator involves the determination of two compartments, infected and non-infected from the model. Recent examples of this method are given in [38, 42, 61, 142] .
Let us assume that there are n compartments of which m are infected. We defined the vector x = (x 1 , . . . , x n ) T , where x i ≥ 0 denotes the number of proportion of individuals in the ith compartment. For clarity we sort the compartments so that the first m compartments correspond to infected individuals. The distinction between infected and uninfected compartments must be determined from the epidemiological interpretation, and not by the mathematical expressions.
It is necessary to define the set The disease transmission model consists of nonnegative initial conditions together with the following system of equations:ẋ
Note that F i should include only infections that are newly arising, but does not include terms which describe the transfer of infectious individuals from one infected compartment to another.
Let us consider the following assumptions:
(A5) If F(x) is set to zero, then all eigenvalues of Df (x 0 ) have negative real parts and Df (x 0 ) is the derivative ∂fi ∂xj evaluated at the DFE x 0 .
Assuming that F i and V i meet the assumptions above, we can form the next generation matrix F V −1 from matrices of partial derivatives of F i and V i . Specially,
where i, j = 1, . . . , m and x 0 is the DFE. The entries of F V −1 give the rate at which infected individuals in x j produce new infections in x i , times the average length of time an individual spends on a single visit to compartment j.
Definition 17 (Basic reproduction number using the next generator operator). The basic reproduction number is given by
where ρ denotes the spectral radius (dominant eigenvalue) of the matrix F V −1 .
The following theorem states that R 0 is a threshold parameter for the stability of the DFE. 
Proof. The proof of this theorem can be found in [42] .
Consider a simple model SEIT for tuberculosis with a treated compartment (adapted from [11] ). Tuberculosis, is a common, and in many cases lethal, infectious disease caused by various strains of mycobacteria. Tuberculosis typically attacks the lungs but can also affect other parts of the body. Most infections are asymptomatic and latent, but about one in ten latent infections eventually progresses to active disease which, if left untreated, kills more than 50% of those so infected.
In this model a constant population is considered where N = S + E + I + T . Exposed individuals progress to the infectious compartment at a rate ν. The treatment rates are r 1 for exposed individuals and r 2 for infectious individuals. However only a fraction q of the treatments of infectious individuals are successful. Unsuccessfully treated infectious individuals re-enter the exposed compartment (1 − q). The dynamics are illustrated in Figure 3 .9 and the differential equations are the following: E and I, which gives m = 2, and we will construct the matrices F and V only related to these compartments. Let
This way,
Epidemiological models can be understood as a framework to explain the mechanisms of the disease procedures and test ideas to implement control measures. R 0 is a key concept, used as a threshold parameter to predict where an infection will spread and how these controls can be effective.
In the second part of the thesis, the epidemiological models will be used for studying Dengue disease. Applying OC theory, the repercussions of several controls in the development of the disease will be analyzed.
During the last decades, the global prevalence of Dengue progressed dramatically. In this chapter, Dengue details are given, such as disease symptoms, transmission and epidemiological trends. An old OC model for Dengue is revisited, as a first approach to the disease. Due to the software robustness improvements and the higher computational capacity, a better solution for this problem is proposed. In order to study different discretization schemes for an OC problem, taking into account time performance and resources used, some numerical simulations are made using Euler and Runge-Kutta methods.
The origins of the word Dengue are not clear. Some researchers think that it is derived from the Swahili phrase "Ka-dinga pepo", meaning "cramp-like seizure caused by an evil spirit". The first recognized Dengue epidemics occurred almost simultaneously in Asia, Africa, and North America in the 1780s, shortly after the identification and naming of the disease in 1779 [33] .
Dengue transcends international borders and can be found in tropical and subtropical regions around the world, predominantly in urban and semi-urban areas. Dengue is a disease which is now endemic in more than one hundred countries of Africa, America, Asia and the Western Pacific. In Figure 4 .1 it is possible to see the areas that in 2008 have more surveillance.
Nevertheless, some studies have indicated that countries with a mild climate, such as in the Mediterranean, are at risk due to future climate conditions that may be favorable to this kind of disease [69] . In Europe there are no registered cases, but the main vector of the disease is already in the old continent and has been followed on Madeira island [5, 124] . This risk may be aggravated further due to climate changes and to the globalization, as a consequence of the huge volume of international tourism and trade [122] . Travelers play an essential role in the global epidemiology: they act as viremic travelers, carrying the disease into areas where mosquitos can transmit the infection. 
Dengue is a vector-borne disease transmitted from an infected human to a female Aedes mosquito by a bite. Then, the mosquito, that needs regular meals of blood to feed their eggs, bites a potential healthy human and transmits the disease making it a cycle.
There are two forms of Dengue: Dengue Fever (DF) and Dengue Hemorrhagic Fever (DHF). The first one is characterized by a sudden high fever without respiratory symptoms, accompanied by intense headaches, painful joints and muscles and lasts between three to seven days. Humans may only transmit the virus during the febrile stage [33] . DHF initially exhibits a similar, if more severe pathology as DF, but deviates from the classic pattern at the end of the febrile stage [54] . The hemorrhagic form has an additionally bleeding from the nose, mouth and gums or skin bruising, nausea, vomiting and fainting due to low blood pressure by fluid leakage. It usually lasts between two to three days and can lead to death [36] . Nowadays, Dengue is the mosquito-borne infection that has become a major international public health concern. According to the World Health Organization (WHO), 50 to 100 million Dengue Fever infections occur yearly, including 500000 Dengue Hemorrhagic Fever cases and 22000 deaths, mostly among children [141] .
There are four distinct, but closely related, viruses that cause Dengue. The four serotypes, named DEN-1 to DEN-4, belong to the Flavivirus family, but they are antigenically distinct. Recovery from infection by one virus provides lifelong immunity against that virus but confers only partial and transient protection against subsequent infection by the other three viruses. There is good evidence that a sequential infection increases the risk of developing DHF [137] .
Unfortunately, there is no specific effective treatment for Dengue. Activities, such as triage and management, are critical in determining the clinical outcome of Dengue. A rapid and efficient front-line response not only reduces the number of unnecessary hospital admissions but also saves lives.
Although up until now there is no effective and safe vaccine for Dengue, a number of candidates are undergoing various phases of clinical trials [140] . With four closely related viruses that can cause the disease, there is a need for a vaccine that would immunize against all four types to be effective. The main difficulty in the vaccine production is that there is a limited understanding of how the disease typically behaves and how the virus interacts with the immune system. Another challenge is that some studies show that some secondary Dengue infection can leave to DHF, and theoretically a vaccine could be a potential cause of severe disease if a solid immunity is not established against the four serotypes. Research to develop a vaccine is ongoing and the incentives to study the mechanism of protective immunity are gaining more support, now that the number of outbreaks around the world is increasing [33] .
The spread of Dengue is attributed to the geographic expansion of the mosquitoes responsible for the disease: Aedes aegypti and Aedes albopictus [23] . Due to its higher interaction with humans and its urban behavior, the first mosquito is considered the major responsible for Dengue transmission and, our attention will be focused on it.
Aedes aegypti, in Figure 4 .2, is an insect species closely associated with humans and their dwellings, thriving in crowded cities and biting primarily during the day. Humans not only provide blood meals for mosquitoes, but also nutrients needed for them to reproduce through water-holding containers, in and around their homes. In urban areas, Aedes mosquitoes breed on water collections in artificial containers such as cans, plastic cups, used tires, broken bottles and flower pots. With increasing urbanization, crowded cities, poor sanitation and lack of hygiene, environmental conditions foster the spread of the disease that, even in the absence of fatal forms, breed significant economic and social costs (absenteeism, immobilization, debilitation and medication) [35] . Dengue is spread only by adult females that require a blood meal for the development of their eggs, whereas male mosquitoes feed on fruit nectar and other sources of sugar. In this process the female acquire the virus while feeding on the blood of an infected person. After virus incubation from eight to twelve days (extrinsic period), an infected mosquito is capable, during probing and blood feeding, of transmitting the virus for the rest of its life to susceptible humans, and the intrinsic period for humans varies from 3 to 15 days.
The life cycle of a mosquito has four distinct stages: egg, larva, pupa and adult, as it is possible to see in Figure 4 .3. In the case of Aedes aegypti, the first three stages take place in or near water whilst air is the medium for the adult stage [103] . Female mosquitoes lay their eggs, but usually do not lay them all at once: it releases them in different places, increasing the probability of new births [139] . The eggs of Aedes aegypti can resist to droughts and low temperatures for up to one year. Although the hatching of mature eggs may occur spontaneously at any time, this is greatly stimulated by flooding. Larvae hatch when water inundates the eggs as a result of rains or an addition of water by people. In the following days, the larvae will feed on microorganisms and particulate organic matter. When the larva has acquired enough energy and size, metamorphosis is done, changing the larva into pupa. Pupae do not feed: they just change in form until the adult body is formed. The newly formed adult emerges from the water after breaking the pupal skin. The entire life cycle, from the aquatic phase (eggs, larvae, pupae) to the adult phase, lasts from 8 to 10 days at room temperature, depending on the level of feeding [27] . The adult stage of the mosquito is considered to last an average of eleven days in an urban environment, reaching up to 30 days in laboratory environment.
Studies suggest that most female mosquitoes may spend their lifetime in or around the houses where they emerge as adults. This means that people, rather than mosquitoes, rapidly move the virus within and between communities.
Aedes aegypti is one of the most efficient vectors for arboviruses because it is highly anthropophilic, frequently bites several times before complete oogenesis [141] . The extent of Dengue transmission is determined by a wide variety of factors: the level of herd immunity in the human population to circulating virus serotype(s); virulence characteristics of the viral strain; survival, feeding behavior, and abundance of aedes aegypti; climate; and human density, distribution, and movement [121] .
It is very difficult to control or eliminate aedes aegypti mosquitoes because they are highly resilient, quickly adapting to changes in the environment and they have the ability to rapidly bounce back to initial numbers after disturbances resulting from natural phenomena (e.g., droughts) or human interventions (e.g., control measures). We can safely expect that transmission thresholds will vary depending on a range of factors.
Primary prevention of Dengue resides mainly in mosquito control. There are two primary methods: larval control and adult mosquito control, depending on the intended target [96] . Larvicide treatment is done through long-lasting chemical in order to kill larvae and preferably have WHO clearance for use in drinking water [36] . The application of adulticides can have a powerful impact on the abundance of adult mosquito vector. However, the efficacy is often constrained by the difficulty in achieving sufficiently high coverage of resting surfaces [37] . This is the most common measure. However the long term use of adulticide has several risks: the resistance of the mosquito to the product reducing its efficacy, the killing of other species that live in the same habitat and has also been linked to numerous adverse health effects including the worsening of asthma and respiratory problems.
Larvicide treatment is an effective way to control the vector larvae, together with mechanical control, which is related to educational campaigns. The mechanical control must be done by both public health officials and residents in affected area. The participation of the entire population is essential in removing still water from domestic recipients and eliminating possible breeding sites [140] .
The most recent approach for fighting the disease is biological control. It is a natural process of population regulation through natural enemies. There are techniques that combine some parasites that kill partially the larval population; however there are some operational resistance, because there is lack of expertise in producing these type of parasites and there are some cultural objections in introducing something in the water for human consumption [23] .
Another way of insect control is to change the reproduction process, releasing sterile insects. This technique, named as Sterile Insect Technique, consists in releasing sterile insects in natural environment, so that the result of mating produces the non-viability of the eggs, and thus can lead to drastic reduction of the specie. This way of control, shows two types of inconvenience: it is expensive to produce and release insects, and can be confronted with social objections, because an uninformed population could not correctly understand the addition of insects as a good solution [8, 47, 131] .
Mathematical modeling became an interesting tool for understanding epidemiological diseases and for proposing effective strategies in fighting them. A set of mathematical models have been developed in literature to gain insights into the transmission dynamics of Dengue in a community. While Zeng and Velasco-Hérnandez [49] investigate the competitive exclusion principle in a two-strain Dengue model, Chowell et al. [26] estimates the basic reproduction number for Dengue using spatial epidemic data. In [100] the author studies the spread of Dengue thought statistical analysis, while in Tewa et al. [130] global asymptotic stability of the equilibrium of a single-strain Dengue model is established. The control of the mosquito by the introduction of a sterile insect technique is analyzed in Thomé et al. [132] . More recently, a study in disease persistence was made [88] in Brazil and Otero et al. [102] studied Dengue outbreaks. All these studies were made with the aim of providing a better understanding of the nature and dynamics of Dengue infection transmission. In the next section a temporal mathematical model that explores the dynamics between hosts (humans) and vectors (mosquitoes) is analyzed.
The aim of this section is to present a mathematical model to study the dynamic of the Dengue epidemics, in order to minimize the investments in disease's control, since financial resources are always scarce. Quantitative methods are applied to the optimization of investments in the control of the epidemiologic disease, in order to obtain a maximum of benefits from a fixed amount of financial resources. The used model depends on the dynamic of the mosquito growing, but also on the efforts of public management to motivate the population to break the reproduction cycle of the mosquitoes by avoiding the accumulation of still water in open-air recipients and spraying potential zones of reproduction.
The Dengue epidemic model described in this paper is based on the one proposed in [20] . It has four state variables and two control variables as follows.
State Variables: To describe the model it is also necessary to introduce some parameters:
Parameters: α R : average reproduction rate of mosquitoes α M : mortality rate of mosquitoes β : probability of contact between non-carrier mosquitoes and infected individuals η :
rate of treatment of infected individuals µ :
amplitude of seasonal oscillation in the reproduction rate of mosquitoes ρ :
probability of individuals becoming infected θ :
fear factor, reflecting the increase in the population willingness to take actions to combat the mosquitoes as a consequence of the high prevalence of the disease in the specific social environment τ :
forgetting rate for goodwill of the target population ϕ : phase angle to adjust the peak season for mosquitoes ω :
angular frequency of the mosquitoes proliferation cycle, corresponding to a 52 weeks period P :
population in the risk area (usually normalized to yield P = 1) γ D : the instantaneous costs due to the existence of infected individuals γ S : the costs of each operation of spraying insecticides γ E :
the cost associated to the instructive campaigns
The model consists in minimizing
subject to the following four nonlinear time-varying state equations [20] :
. . , 4. Equation (4.2a) represents the variation of the mosquitoes density per unit time to the natural cycle of reproduction and mortality (α R and α M ), due to seasonal effects µ sin(ωt + ϕ) and to human interference −x 4 (t) and u 1 (t). Equation (4.2b) expresses the variation of the mosquitoes density carrying the virus
x 2 (t) denotes the rate of the infected mosquitoes and β [x 1 (t) − x 2 (t)] x 3 (t) represents the increase rate of the infected mosquitoes due to the possible contact between the uninfected mosquitoes x 1 (t) − x 2 (t) and infected individuals denoted by x 3 (t). The dynamics of the infectious transmission is presented in equation (4.2c). The term −ηx 3 (t) is related to the rate of cure and ρx 2 (t) [P − x 3 (t)] describes the rate at which new cases spring up. The factor [P − x 3 (t)] is the number of individuals in the area, that are not infected. Equation (4.2d) is a model for the level of popular motivation (or goodwill) to combat the reproductive cycle of mosquitoes. Over time, the level of people motivated will have changed. As a consequence, it is necessary to invest in educational campaigns designed to increase consciousness of the population under risk. The expression −τ x 4 (t) represents the decay of the people's motivation over time, due to forgetfulness. The term θx 3 (t) describes the natural sensibilities of the public due to increase in the prevalence of the disease.
The goal is to minimize the cost functional (4.1). This functional includes social costs related to the existence of ill individuals -like absenteeism, hospital admission, treatments -, γ D x 2 3 (t), the recourses needed for the spraying of insecticide operations, γ S u 2 1 (t), and for educational campaigns, γ E u 2 2 (t). The model for the social cost is based on the concept of goodwill explored by Nerlove and Arrow [99] .
Due to computational issues, the optimal control problem (4.1)-(4.2d), that was written in the Lagrange form, was converted into an equivalent Mayer problem. Hence, using a standard procedure (cf., Section 1.2) to rewrite the cost functional, the state vector was augmented by an extra component x 5 ,
leading to the following equivalent terminal cost problem:
with given t f , subject to the control system (4.2a)-(4.2d) and (4.3).
Two different implementations were considered. In a first approach, the OC problem is solved by a specific Optimal Control package, OC-ODE, already described in Section 2.3. The OC problem considers (4.2a)-(4.4) and the code is available in [110] .
A second approach uses the nonlinear solver Ipopt [135] , also described in Section 2.3. In order to use this software, it was necessary to discretize the problem. The Euler discretization scheme was chosen (see Section 2.1 for more details). The discretization step length was h = 1/4, because it is a good compromise between precision and efficiency. Thus, the optimal control problem was discretized into the following nonlinear programming problem:
The error tolerance value was 10 −8 using the Ipopt solver. The discretized problem, after a presolve done by the software, has 1455 variables, 1243 of which are nonlinear; and 1039 constraints, 828 of which are nonlinear (see the AMPL code for this problem in [110] ). [104, 133] that was used by the authors of the paper [20] . It is important to salient that, at the time of the initial paper [20] , the authors had not the same computational resources that exist nowadays. The results with OC-ODE and Ipopt are better since the cost to fight the Dengue disease and the number of infected individuals are smaller. [20] . packages were tested, and they could not reach a solution -some crashed at the middle or some bad scaling issues were observed. Until some years ago, due to computational limitations, most of the models were run using codes made by the authors themselves as [20] . Nowadays, one can choose between several proper software packages "out of the box", that already take into account specific features of stiff problems, scaling problems, etc. With this work it is possible to realize that "old" problems can again be taken into account and be better analyzed with new technology and approaches, with the goal of finding global optimal solutions, instead of local ones.
For this purpose, and at an initial research stage, it is important to understand if different kinds of discretization for an OC problem have influence on the problem resolution.
This section aims to study the costs of different discretization processes, in terms of time performance, number of variables and iterations used. For the purpose of this analysis two discretization schemes, Euler and second order Runge-Kutta's scheme [9] (cf. Section 2.1), are considered to solve the problem described in the previous section.
This discretization process transforms the Dengue epidemics problem into a standard nonlinear optimization problem, with an objective function and a set of nonlinear constraints. This NLP problem was codified, for both discretization schemes, in the AMPL modeling language [55] and can be checked in [110] .
Two nonlinear solvers with distinct features were selected to solve the NLP problem: the Knitro (IP method) and the Snopt (SQP method). The NEOS Server [98] platform was used as interface with both solvers. The Ipopt (IP method), used in the previous section, was the first choice for our research. However, at the time of this investigation, the NEOS platform was moved to another Center of Research and some software packages were unavailable for long periods of time. So, we had to choose another Interior Point robust software. Table 4 .1 reports the results for both solvers, for each discretization method using three different discretization steps (h = 0.5, 0.25, 0.125), rising twelve numerical experiences. The columns # var. and # const. mean the number of variables and constraints, respectively. The next columns refer to the performance measures -number of iterations and total CPU time in seconds (time for solving the problem, for evaluate the objective and the constraints functions and for input/output). As the computational experiments were made in the NEOS server platform, the selected machine to run the program remains unknown as well as its technical specifications.
The optimal cost reached was ≈ 3 × 10 −3 for all tests. Comparing the general behavior of the solvers one can conclude that the IP based method (Knitro) presents much better performance than the SQP method (Snopt) in terms of the measures used. Regarding the Knitro results, one realizes that the Euler's discretization scheme has better times for h = 0.25 and h = 0.125 and similar time for h = 0.5, when compared to Runge-Kutta's method. Another obvious finding, for both solvers, is that the CPU time increases as far as the problem dimension increases (number of variables and constraints). With respect to the number of iterations, Snopt presents more iterations as the problem dimension increases. However this conclusion cannot be taken for Knitro -in fact, there does not exist a relation between the problem dimension and the number of iterations. The best version tested was Knitro using Runge-Kutta with h = 0.5 (best CPU time and fewer iterations), and the second one was Knitro with Euler's method using h = 0.25. An important evidence of this numerical experience is that it is not worth the reduction of the discretization step size because no significative advantages are obtained.
At this moment, as a result of major demographic changes, rapid urbanization on a massive scale, global travel and environmental change, the world faces enormous future challenges from emerging infectious diseases. Dengue illustrates well these challenges [141] . In this work we investigated an optimal control model for Dengue epidemics proposed in [20] , that includes the mosquitoes dynamics, the effect of educational campaigns. The cost functional reflects a compromise between financial spending in insecticides and educational campaigns and the population health. For comparison reasons, the same choice of data/parameters in [20] was considered.
The results obtained from OC-ODE and Ipopt are similar, improving the ones previously reported in [20] (cf. Section 4.2). Indeed, the obtained control policy in this work presents an important progress with respect to the previous best policy: the percentage of infected mosquitoes vanishes just after four weeks, while mosquitoes are completely eradicated after 30 weeks (Figures 4.4 and 4 .5); the number of infected individuals begin to decrease after four weeks while with the previous policy this only happens after 23 weeks (Figure 4.6) . Despite the fact that our results are better, they are accomplished with a much smaller cost with insecticides and educational campaigns (Figure 4.8) . The general improvement, which explains why the results are so successful, rely on an effective control policy of insecticides. The proposed strategy for insecticide application seems to explain the discrepancies between the results here obtained and the best policy of [20] . Our results show that applying insecticides in the first four weeks yields a substantial reduction in the cost of fighting Dengue, in terms of the functional proposed in [20] . The main conclusion is that health authorities should pay attention to the epidemic from the very beginning: effective control decisions in the first four weeks have a decisive role in the Dengue battle, and population and governments will both profit from it.
We successfully solved an OC problem by direct methods using nonlinear optimization software based on IP and SQP approaches. The problem was discretized through Euler and Runge-Kutta schemes. The implementation efforts of higher order discretization methods bring no advantages. The reduction of the discretization step and consequently the increase of the number of variables and constraints do not improve the performance with respect to the CPU time and to the number of iterations. We can point out the robustness of both solvers in spite of the dimension problem increase. The conclusions drawn in Section 4.3 were helpful for decision-making future processes of discretization all over the work.
As future work we intend to analyze how different parameters/weights associated to the variables in the objective function can influence the spread of the disease. This chapter was based on work available in the peer reviewed journal [112] and the peer reviewed conference proceedings [111] .
A model for the Dengue disease transmission is presented. It consists of eight mutually-exclusive compartments representing the human and vector dynamics. It also includes a control parameter, adulticide spray, as a measure to fight the disease. The model presents three possible equilibria: two Disease Free Equilibria (DFE) and an Endemic Equilibrium (EE). It has been proved that a DFE is locally asymptotically stable, whenever a certain epidemiological threshold, known as the basic reproduction number, is less than one. In this work we try to understand which is the best way to apply the control in order to effectively reduce the number of infected humans and mosquitoes. A case study, using outbreak data in 2009 in Cape Verde, is reported.
In Chapter 4 the Dengue epidemic was studied, mostly centered in people, specially in the goodwill of the individuals and spraying campaigns. However, the virus transmission scheme was overlooked: it was only considered two compartments for people and two compartments for adult mosquitoes. Here, the aim is to deepen the relationship between human and mosquitoes, creating a better framework to explain the development and transmission of the disease.
The mathematical model is based on [43, 44] , that describes the chikungunya disease transmitted by Aedes albopictus.
The notation used in our mathematical model includes four epidemiological states for humans:
It is assumed that the total human population (N h ) is constant, so, N h = S h + E h + I h + R h . There are also four other state variables related to the female mosquitoes (the male mosquitoes are not considered in this study because they do not bite humans and consequently they do not influence the dynamics of the disease):
Similarly, it is assumed that the total adult mosquito population is constant, which means N m = S m + E m + I m . In this way, we put our model more complex and reliable to the reality of Dengue epidemics. For this study we introduced a control variable:
The control variable, c(t), varies from 0 to 1. However, the model does not fit completely the reality. Epidemiologist and policy makers need to be aware of both strengths and weakness of the epidemiological modeling approach. An epidemiological model is always a simplification of reality. So, some assumptions were made to built this model:
• the total human population (N h ) is constant;
• there is no immigration of infected individuals into the human population;
• the population is homogeneous, which means that every individual of a compartment is homogenously mixed with the other individuals;
• the coefficient of transmission of the disease is fixed and does not vary seasonally;
• both human and mosquitoes are assumed to be born susceptible, i.e., there is no natural protection;
• for the mosquito there is no resistant phase, due to its short lifetime.
To completely describe the model it is necessary to use parameters, which are:
total population B :
average daily biting (per day) β mh :
transmission probability from I m (per bite) β hm :
transmission probability from I h (per bite) 1/µ h : average lifespan of humans (in days) 1/η h : mean viremic period (in days) 1/µ m : average lifespan of adult mosquitoes (in days) ϕ :
number of eggs at each deposit per capita (per day) µ A : natural mortality of larvae (per day) η A : maturation rate from larvae to adult (per day) 1/η m : extrinsic incubation period (in days) 1/ν h :
intrinsic incubation period (in days) m :
female mosquitoes per human k :
number of larvae per human K :
maximal capacity of larvae For notation simplicity, the independent variable t will be omitted when writing the dependent variables, example given, will be written S h instead of S h (t). The Dengue epidemic can be modelled by the following nonlinear time-varying state equations:
and vector population
with the initial conditions For the previous set of differential equations is now analyzed the equilibrium points of the system and is determined the threshold phenomena.
Let the set Proof. System (5.1)-(5.2) can be rewritten in the following way:
As M (X) has all off-diagonal entries nonnegative, M (X) is a Metzler matrix.
Using the fact that F ≥ 0, the system (5.4) is positively invariant in R 7 + [1] , which means that any trajectory of the system starting from an initial state in the positive orthant R 7 + remains forever in R 7 + .
Theorem 5.1.1. Let Ω be defined as above. Consider also
The system (5.1)-(5.2) admits at most two Disease Free Equilibrium points:
• if M ≤ 0, there is a Disease Free Equilibrium (DFE),
called Trivial Equilibrium;
• if M > 0, there is a Biologically Realistic Disease Free Equilibrium (BRDFE),
Proof. The equilibrium points are reached when the following equations hold:
Using the Mathematica software to solve the system (5.5), we obtained four solutions.
The first one is known as the Trivial Equilibrium, since the mosquitoes do not exist, so there is no disease:
In the second one, mosquitoes and humans interact, but there is only one outbreak of the disease, i.e., over time the disease goes away without killing all the mosquitoes. We have called this equilibrium point a Biologically Realistic Disease Free Equilibrium (BRDFE), since it is a more reasonable situation to find in nature than the previous one:
which is equivalent to E * 2 = N h , 0, 0, kN hM η A µ b , kN hM µ b µm , 0, 0 . This is biologically interesting only if M is greater than 0.
The third solution corresponds to a situation where humans and mosquitoes live together but the disease persists in both populations, which means that it is not a DFE. This equilibrium will be explained after (see Theorem 5.1.3). Thus the disease is not anymore an epidemic episode, but transforms into a endemic one. With some algebraic manipulations we obtained the following point:
With the Mathematica software we obtained a fourth solution. But some of the components are negative, which means that they do not belong to the Ω set. Remark 6. The condition M > 0 is equivalent, by algebraic manipulation, to the condition
which corresponds to the basic offspring number for mosquitoes. Thus, if M < 0, then the mosquito population will collapse and the only equilibrium for the whole system is the trivial equilibrium. If M ≥ 0, then the mosquito population is sustainable.
The amount of mosquitoes is also related to an epidemic threshold: the basic reproduction number of the disease, R 0 . Following [42] , we prove: 
.
The equilibrium point BRDFE is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.
Proof. To derive the basic reproduction number, we use the next generator approach. The basic reproduction number is calculated in a Disease Free Equilibrium. In this case we consider the most realistic one, BRDFE.
Following [3, 42] , let consider the vector x T = (E h , I h , E m , I m ) which corresponds to the components related to the progression of the disease.
The subsystem used is:
This subsystem can be partitioned,
represents the components related to new cases of disease (in this situation in the exposed compartments) and V(x) represents the other components. Thus the subsystem (5.6) can be rewritten as
Let us consider the Jacobian matrices associated with F and V:
According to [42] the basic reproduction number is R 0 = ρ(J F (x0) J V −1 (x0) ), where x 0 is a Disease Free Equilibrium (BRDFE) and ρ(A) defines the spectral radius of a matrix A. Using Mathematica
and we obtain the value for the threshold parameter, with M > 0.
In the model we have two different populations (human and vectors), so the expected basic reproduction number reflects the infection human-vector and also vector-human, that is, R 2 0 = R hm × R mh . The term Bβ hm Sm0 N h represents the product between the transmission probability of the disease from humans to vectors and the number of susceptible mosquitoes per human; 1 η h +µ h is related to the human's viremic period; and ηm c+ηm+µm represents the proportion of mosquitoes that survive to the incubation period. Analogously, the term Bβ mh S h0 N h is related to the transmission probability of the disease between mosquitos and human, in a susceptible population; When R 0 < 1, each infected individual produces, on average, less than one new infected individual, and therefore it is predictable that the infection will be cleared from the population. If R 0 > 1, the disease is able to invade the susceptible population [42, 83] .
Theorem 5.1.3. If M > 0 and R 0 > 1, then the system (5.1)-(5.2) also admits an endemic equilibrium (EE):
.
Proof. See Proof of Theorem 5.1.1.
From a biological point of view, it is desirable that humans and mosquitoes coexist without the disease reaching a level of endemicity. Good estimates of Dengue transmission intensity are therefore necessary to compare and interpret Dengue interventions conducted in different places and times and to evaluate options for Dengue control. The basic reproduction number has played a central role in the epidemiological theory for Dengue and other infectious diseases because it provides an index of transmission intensity and establishes a threshold criteria. We claim that proper use of control c, can result in the basic reproduction number remaining below unity and, therefore, making BRDFE stable. In order to make effective use of achievable insecticide control, and simultaneously to explain more easily to the competent authorities its effectiveness, we assume that c is constant. The goal is to find c such that R 2 0 < 1. For this purpose we have studied the reality of Cape Verde.
An unprecedented outbreak was detected in the Cape Verde archipelago in September 2009. This is the first report of Dengue virus activity in that country. As the population had never had contact with the virus, the herd immunity was very low. Dengue type 3 spread throughout the islands of the archipelago, reaching four of the nine islands. The worst outbreak occurred on the Santiago island, where most people live. The number of cases increased sharply since the beginning of November, reaching 1000 cases per day. The Cape Verde Ministry of Health reported more than 20000 cases of Dengue Fever within the archipelago between October and December 2009, which is about 5% of the total population of the country. From 173 cases of Dengue hemorrhagic fever reported, six people died [31, 40] .
It represented a challenge to the performance of the National Health Care System. Government officials launched a plan to eradicate the mosquito including a national holiday during which citizens were being asked to clear out standing water and other potential breeding areas used by the mosquitoes. The intense and spontaneous movement of solidarity from the civil society was another noteworthy dimension, not only cleaning but also voluntary donating blood to strengthen the stock of the Central Hospital Dr. Agostinho Neto, in Praia.
We used the data for human population related to Cape Verde [54] . Due to low surveillance and the fact of being the first ever of Dengue outbreak in the country, it was not possible to collect detailed data from the mosquito. However, the authorities speak explicitly about mosquitoes coming from Brazil [2] . Also the information that comes from the Ministry of Health in the capital of Cape Verde, Praia, confirms that the insects responsible for Dengue came most probably from Brazil, transported by means of air transport that perform frequent connections between Cape Verde and Brazil, as reported by the Radio of Cape Verde. With respect to Aedes aegypti, we have thus considered data from Brazil [132, 143] .
The simulations were carried out using the following values: N h = 480000, B = 1, β mh = 0.375, β hm = 0.375, µ h = 1/(71 × 365), η h = 1/3, µ m = 1/11, ϕ = 6, µ A = 1/4, η A = 0.08, η m = 1/11, ν h = 1/4, , m = 6, k = 3. The initial conditions for the problem were:
Considering nonexistence of control, i.e. c = 0, the basic reproduction number for this outbreak in Cape Verde is approximately R 0 = 2.396, which is in agreement to other studies of Dengue in other countries [100] . The control c affects the basic reproduction number, and our aim is to find a control that puts R 0 less than one. This value was obtained by Mathematica, calculating the inequality with the parameters values above. The computational investigations were carried out using c = 0.157, which means that the insecticide is continuously applied during a twelve week period.
The software used was Scilab [21] . It is an open source, cross-platform numerical computational package and a high-level, numerically oriented programming language. For our problem we used the routine ode to solve the set of differential equations. By default, ode uses the lsoda solver of package ODEPACK. It automatically selects between the nonstiff predictor-corrector Adams method and the stiff Backward Differentiation Formula (BDF) method. It uses the nonstiff method initially and dynamically monitors data in order to decide which method to use. The graphics were also obtained using this software, with the command plot (see code in [110] ). Figures 5.2 and 5.3 show the curves related to human population, with and without control, respectively. The number of infected people, even with small control, is much less than without any insecticide campaign. Figures 5.4 and 5.5 show the difference of the mosquito population with control and without control. When the control is applied, the number of infected mosquitoes is close to zero. Note that the intention is not to completely eradicate the mosquitoes but instead the number of infected mosquitoes. 
It has been algebraically proved that if a constant minimum level of a control is applied (c = 0.157), it is possible to maintain the basic reproduction number below unity, guaranteeing the BRDFE. This value is corroborated in another numerical study [113] .
The values of infected humans obtained by the model are higher when compared to what really happened in Cape Verde. Despite the measures taken by the government and Health authorities were not accounted, they have had a considerable impact in the progress of the disease.
Until here, it was considered a constant control. Using a theoretical approach [112] , we intend to find the best function c(t), using OC theory. Instead of finding a constant control it will be possible to study other types of control, such as piecewise constant or even continuous but not constant. Additionally we could consider another strategy, a more practical one: due to logistics and health reasons, it may be more convenient to apply insecticide periodically and at some specific hours at night.
In this section we investigate the best way to apply the control in order to effectively reduce the number of infected humans and mosquitoes, using pulse control.
In literature it has been proven that a DFE is locally asymptotically stable, whenever a certain epidemiological threshold, known as the basic reproduction number , is less than one. In the previous section, it was proven that if a constant minimum level of insecticide is applied (c = 0.157), it is possible to maintain the basic reproduction number below unity, guaranteeing the DFE. In this section, other kinds of piecewise controls that maintain the basic reproduction number less than one and could give a better solution to implement by health authorities are investigated.
To solve the system (5.1)-(5.2), in a first step, several strategies of control application were used: three different frequencies (weekly, bi-weekly and monthly) for control application, a constant control (c = 0.157) and no control (c = 0). The first three different frequencies mean that during one day (per week, bi-week and month), the whole (100%) capacity of insecticide (c = 1) is used during all day. Besides, also was used a constant control strategy (c = 0.157) that consists in the application of 15.7% capacity of insecticide 24 hours per day all the time (84 days). In this work, the amount of insecticide is an adimensional value and must be considered in relative terms.
The numerical tests were carried out using Scilab [120] , with the same ODE system and parameters of the previous section (code available in [110] ). Figures 5.6 and 5.7 show the results of these strategies regarding infected mosquitoes and individuals. Without control, the number of infected mosquitoes and individuals increase expressively. It is also possible to see that the weekly pulse control had the closest results to the continuous control.
Therefore, realizing the influence of the insecticide control, further tests were carried out to find the optimum periodicity of administration which, from gathered results, must rest between six and seven days. The second phase of numerical tests, Figures 5.8 and 5 .9, considered four situations: 6 days, 7 days, 10 days and continuously c = 0.157. To guarantee the DFE, the curves must remain below the one corresponding to c = 0.157.
The amount of insecticide, and when to apply it, are important factors for outbreak control. Table 5 .1 reports the total amount of insecticide used in each version during the 84 days.
The numerical tests conclude that the best strategy for the infected reduction is every 6/7 days application. The value spent in insecticide is similar to the continuous strategy, but is much easier to implement In this work, several piecewise strategies were studied to find the best way of applying insecticide, but always having in mind a periodic application of the product. But if the best strategy is not a periodic one? In the next section the optimal control solution for this problem will be studied.
The aim of this section is the study of optimal strategies for applying insecticide, taking into account different perspectives: thinking only on the insecticide cost, focusing on the cost of infected humans or even combining both perspectives.
To take this approach, and after several numerical experiences, we considered that it was better to normalize the ODE system (5.1)-(5.2). The reason for this transformation is due to the bad scaling of the variables: some vary from 0 to 480000, and others from 0 to 1, which can influence the software performance.
Consider the following transformations:
The ODE system (5.1)-(5.2) is transformed into:
with the initial conditions Let us consider the objective functional J considering the costs of infected humans and costs with insecticide:
where γ D and γ S are positive constants representing the costs weights of infected individuals and spraying campaigns, respectively.
Using OC theory, let λ i (t), with i = 1, . . . , 8, be the co-state variables. The Hamiltonian for the present OC problem is given by
By the Pontryagin Maximum Principle [106] , the optimal control c * should be the one that minimizes, at each instant t, the Hamiltonian given by (5.10), that is, H (x * (t), λ * (t), c * (t)) = min c∈[0,1] H (x * (t), λ * (t), c). In this way, the optimal control is given by
It is also necessary to consider the adjoint system λ i (t) = − ∂H ∂x i , i.e.,
As the OC problem just has initial conditions it is necessary to find the transversality conditions, that correspond to a terminal condition in the co-state equation,
Replacing the optimal control c * in the state system (5.7) and in the adjoint system (5.11) is possible to solve the differential system taking into account the initial and transversality conditions.
In order to solve this OC problem three approaches were tested. The first one is a direct method,
DOTcvp toolbox [67] for Matlab. It uses the differential system (5.7), the initial conditions (5.8) and it is necessary to transform the problem into a Mayer form (see Section 2.3). To solve the discretized OC problem, the Ipopt software is chosen as an option inside the DOTcvp. The functional is divided into time intervals, in this case ten intervals, with an initial value problem tolerance of 10 −7 and a NLP tolerance of 10 −5 . The optimal control function given by this toolbox is piecewise constant. Another direct method is OC-ODE [57] . It uses the differential system (5.7), the initial conditions (5.8) and it includes procedures for numerical adjoint estimation and sensitivity analysis and the feasibility tolerance considered was 10 −10 .
The last method used, an indirect one, it is coded in Matlab environment. It involves the backwardforward method (see Section 2.2) and the resolution of the ODE systems are made by the ode45 routine. The state differential system (5.7) is solved forward with initial conditions (5.8), while the adjoint system (5.11) is solved backwards using the terminal conditions (5.12). The absolute and relative tolerances were fixed in 10 −4 . The three codes are available at [110] .
All the parameters are assumed equal to the previous section. For this first simulation, the values for the weights of the costs are γ D = 0.5 and γ S = 0.5. Figure 5 .10 shows that, despite having distinct philosophies of resolution, the curves obtained by the three solvers are similar, which reinforces the confidence in the result. The optimal functional is 0.0470, 0.0498 and 0.0445, for DOTcvp, OC-ODE and Backward-Forward, respectively. The last solver achieved better value for the functional. This is expected, once its resolution contains more information about the OC problem because the adjoint system is supplied. The study considers three situations: A, B, and C. Situation A, that was previously presented, regards both perspectives in the functional (human infected and insecticide application). Situation B concerns only to infected humans whereas case C only considers insecticide campaigns. Table 5 .2 resumes the three situations.
Functional Values for weights Case A Both perspectives Since the three solvers presented similar solutions, only one of them was chosen to solve these cases, that was DOTcvp.
Figures 5.12, 5.13 and 5.14 show the results for optimal control c, infected humans and total costs, respectively.
In the medical perspective (case B), when only the costs related with ill people (absenteeism, drugs, ...) are considered, the number of infected is the lowest; however a huge quantity of insecticide is used, because it is considered cheaper. On the other hand, when people just think on the economical perspective (case C), the treatment for people is neglected. Optimal control is low, but the number of infected humans is high. The total cost is higher when both perspectives are considered.
Epidemiological modelling has largely focused on identifying the mechanisms responsible for epidemics but has taken little account on economic constraints in analyzing control strategies. Economic models have given insight into optimal control under constraints imposed by limited resources, but they frequently ignore the spatial and temporal dynamics of the disease. Nowadays the combination of epidemiological and economic factors is essential. The bioeconomic approach requires an equilibrium between economic and epidemiological parameters in order to give an efficient disease control and reflecting the nature of the epidemic. For this, the study goes on implementing both perspectives, but taking into account distinct weights of the parameters associated to the variables i h and c. Table 5 .3 resumes these approaches.
Values for weights Table 5 .3: Different weights for the functional In case D, it is studied a situation where the lack of insecticide in a country could be a reality and as a consequence its market value is high. This could happen due to an unprecedent outbreak where the authorities were not prepared or even due to financial reasons by the fact that the government does not have financial viability for this kind of measure. In case E, once again the human perspective gains strength and, as human life and quality of life is an expensive good, it was considered that it is more expensive to treat humans then apply insecticide.
The analysis from the Figures from 5.15 to 5.17 is consistent with what we expected in reality. In case D, as insecticide is expensive, the function for optimal control is lower than the other perspective. As consequence, the number of infected people is higher.
In case E, where the human factor is preponderant, the number of infected humans is low but expenses with insecticide are higher. Curiously the total cost, in case D and E, are of the same order of magnitude, with a slightly higher cost for case E. The total cost is reported in Table 5 .4. When both perspectives are considered, the total cost is higher than the single perspective.
For a last analysis, a mathematical perspective was carried out: what values should the γ D and γ S have, in order to minimize the functional? Let us call this perspective as Case F. In this perspective, we not only want to minimize the control c, but also the parameters γ D and γ S , enforcing the equality constraint γ D + γ S = 1. show the comparison of this case with the first one. As expected, the total cost and the number of infected humans are the lowest ones. Giving freedom to the parameters it is possible to see that the optimal control function is not periodic (as studied in the previous section), but gives a practical solution in applying insecticide.
In this chapter a model based on two populations, humans and mosquitoes, with insecticide control is presented. It has shown that as time goes by, depending on several parameters, the outbreak could disappear (leading to a DFE) or could transform the disease into an endemic one (leading to an EE). This Assuming that the parameters are fixed, the only variable that can influence this threshold is the control variable c, it has shown that with a steady insecticide campaign it is possible to reduce the number of infected humans and mosquitoes, and can prevent an outbreak that could transform an epidemiological episode to an endemic disease.
For a steady campaign has proven that c = 0.157 it is enough to maintain the R 0 below unit. However, this type of control is difficult to implement. A pulse insecticide campaign was studied to circumvent this difficulty. It has proven that applying insecticide every 6/7 days, is a better strategy to implement by health authorities and has the same efficacy level and financial costs.
Finally, the OC theory was used to find the best optimal control function for the insecticide. The optimal function varies, giving a different answer depending on the main goal to reach, thinking in economical or human centered perspective.
As future work it is important to study different kinds of controls. The accelerated increase in mosquito resistance to several chemical insecticides and the damage caused by these to the environment, has resulted in the search for new control alternatives. Among the alternatives available, the use of Bacillus thuringiensis israelensis (Bti) has been adopted by several countries [89] . Laboratory testing shows that Bti has a high larvicide property and its mechanism of action is based on the production of an endotoxin protein that, when ingested by the larvae, causes death. To ensure the minimization of the outbreaks, educational programmes that are customized for different levels of health care and that reflect local capacity should be supported and implemented widely. People should be instructed to minimize the number of potential places for the mosquito breeding. Educational campaigns can be included as an extra control parameter in the model. This chapter was based on work available in the peer reviewed journal [118] and the peer reviewed conference proceedings [113, 117] .
A new model with six mutually-exclusive compartments related to Dengue disease is presented. In this model there are three vector control tools: insecticides (larvicide and adulticide) and mechanical control. The human data for the model is again related to Cape Verde. Due to the rapid development of the outbreak on the islands, only a few control measures were taken, but not quantified. In this chapter, some of these measures are simulated and their consequences are analyzed.
In Chapter 5, a model with eight compartments and a single control was analyzed. However, after discussion with some researchers in this area, many of them suggested removing the compartment exposed for three main reasons: first, it is difficult to collect data for this compartment, since the disease at this stage does not show symptoms; second, the curve obtained is similar to the infected compartment with only an advance of time, not bringing novelty to the model, but possible difficulties to the numeric resolution; and finally, as the main goal is to study the effects of several controls centered in the infected humans, this compartment plays a secondary role. Thus, it was decided to remove the exposed compartments, in humans and mosquitoes, adjusting the other parameters to this new model and including three controls.
Taking into account the model presented in [43, 44] and the considerations of [111] - [113] , a new model more adapted to the Dengue reality is proposed. The notation used in the mathematical model includes three epidemiological states for humans:
It is assumed that the total human population (N h ) is constant and N h = S h (t) + I h (t) + R h (t) at any time t. The population is homogeneous, which means that every individual of a compartment is homogeneously mixed with the other individuals. Immigration and emigration are not considered.
There are three other state variables, related to the female mosquitoes:
Due to the short lifespan of mosquitoes, there is no resistant phase. It is assumed homogeneity between host and vector populations, which means that each vector has an equal probability to bite any host. Humans and mosquitoes are assumed to be born susceptible.
To analyze the effect of campaigns in the disease fight, three controls are considered:
proportion of mechanical control, 0 < α ≤ 1.
Larval control targets the immature mosquitoes living in water before they become biting adults. A soil bacterium, Bacillus thuringiensis israelensis (Bti), is applied from the ground or by air to larval habitats. This bacterium is used because when properly applied, it has virtually no effect on non-target organisms.
The control of adult mosquitoes is necessary when mosquito populations cannot be treated in their larval stage. It is the most effective way to eliminate adult female mosquitoes that are infected with human pathogens. Depending on the size of the area to be treated, either trucks for ground adulticide treatments or aircraft for aerial adulticide treatments can be used.
The purpose of mechanical control is to reduce the number of larval habitat areas available to mosquitoes. The mosquitoes are most easily controlled by treating, cleaning and/or emptying containers that hold water, since the eggs of the specie are laid in water-holding containers.
The aim is to simulate different realities in order to find the best policy to decrease the number of infected humans. A temporal mathematical model is introduced, with mutually-exclusive compartments, to study the outbreak occurred on Cape Verde islands in 2009 and improving the model described in [111] .
The model uses the following parameters:
average daily biting (per day) β mh :
transmission probability from I m (per bite) β hm :
transmission probability from I h (per bite) 1/µ h : average lifespan of humans (in days) 1/η h : mean viremic period (in days) 1/µ m : average lifespan of adult mosquitoes (in days) ϕ :
number of eggs at each deposit per capita (per day) 1/µ A : natural mortality of larvae (per day) η A : maturation rate from larvae to adult (per day) m :
female mosquitoes per human k :
number of larvae per human
The Dengue epidemic is modelled by the following nonlinear time-varying state equations:
with the initial conditions 
Due to biological reasons, only nonnegative solutions of the differential system are acceptable. More precisely, it is necessary to study the solution properties of the system (6.1)-(6.2) in the closed set
It can be verified that Ω is a positively invariant set with respect to (6.1)-(6.2). The proof of this statement is similar as in [118] . The system (6.1)-(6.2) has at most three biologically meaningful equilibrium points (cf. Theorem 6.1.1).
is said to be an equilibrium point for system (6.1)-(6.2) if it satisfies the following relations:
An equilibrium point E is biologically meaningful if and only if E ∈ Ω. The biologically meaningful equilibrium points are said to be disease free or endemic depending on I h and I m : if there is no disease for both populations of humans and mosquitoes (I h = I m = 0), then the equilibrium point is said to be a Disease Free Equilibrium (DFE); otherwise, if I h = 0 or I m = 0 (in other words, if I h > 0 or I m > 0), then the equilibrium point is called endemic. Theorem 6.1.1. System (6.1)-(6.2) admits at most three biologically meaningful equilibrium points; at most two DFE points and at most one endemic equilibrium point. More precisely, let
,
.
If M ≤ 0, then there is only one biologically meaningful equilibrium point, E * 1 , which is a DFE point. If M > 0 with ξ ≥ χ, then there are two biologically meaningful equilibrium points, E * 1 and E * 2 , both DFE points. If M > 0 with ξ < χ, then there are three biologically meaningful equilibrium points, E * 1 , E * 2 , and E * 3 , where E * 1 and E * 2 are DFEs and E * 3 endemic.
Proof. System (6.3) has four solutions easily obtained with a Computer Algebra System like Maple: E * 1 , E * 2 , E * 3 and E * 4 . The equilibrium point E * 1 is always a DFE because it always belongs to Ω with I h = I m = 0. In contrast, E * 4 is never biologically realistic because it has always some negative coordinates. The other two equilibrium points, E * 2 and E * 3 , are biologically realistic only for certain values of the parameters. The equilibrium E * 2 is biologically realistic if and only if M ≥ 0, in which case it is a DFE. For the condition M ≤ 0, the third equilibrium E * 3 is not biologically realistic. If M > 0, then three situations can occur with respect to E * 3 : if ξ = χ, then E * 3 degenerates into E * 2 , which means that E * 3 is the DFE E * 2 ; if ξ > χ, then E * 3 is not biologically realistic; otherwise, one has E * 3 ∈ Ω with I h = 0 and I m = 0, which means that E * 3 is an endemic equilibrium point.
By algebraic manipulation, M > 0 is equivalent to condition
is related to the basic offspring number for mosquitos. Thus, if M ≤ 0, then the mosquito population will collapse and the only equilibrium for the whole system is the trivial DFE E * 1 . If M > 0, then the mosquito population is sustainable. From a biological standpoint, the equilibrium E * 2 is more plausible, because the mosquito is in its habitat, but without the disease.
An important measure of transmissibility of the disease is now introduced: the basic reproduction number. It provides an invasion criterion for the initial spread of the virus in a susceptible population. For this case the following result holds. Theorem 6.1.2. The basic reproduction number R 0 associated to the differential system (6.1)-(6.2) is
Proof. In agreement with [42] , just the epidemiological compartments that have new infections, I h and I m , are considered. The two differential equations related to these two compartments can be rewritten
where F is the rate of production of new infections and V is the transition rates between states:
The Jacobian derivatives are
The quantity J F (x) J −1 V(x) gives the total production of new infections over the course of an infection. The largest eigenvalue gives the fastest growth of the infected population, which means that R 0 is the spectral radius of the matrix J F (x) J −1 V(x) in a DFE point. Maple was used to obtain
The basic reproduction number R 0 in (6.5) is obtained, replacing S h DF E and S m DF E in (6.6) by those of the DFE E * 2 .
The model has two different populations (host and vector) and the expected basic reproduction number should reflect the infection transmitted from host to vector and vice-versa. Accordingly, R 0 can be seen as R 0 = (R hm × R mh ) If R 0 < 1, then, on average, an infected individual produces less than one new infected individual over the course of its infectious period, and the disease cannot grow. Conversely, if R 0 > 1, then each individual infects more than one person, and the disease can invade the population. Proof. The only solution of (6.3) with I h > 0 or I m > 0, the only endemic equilibrium, is E * 3 . That occurs, in agreement with Theorem 6.1.1, in the case M > 0 and χ > ξ. The condition χ > ξ is equivalent, by Theorem 6.1.2, to R 0 > 1.
Using the methods in [42, 83] , it is possible to prove that if R 0 ≤ 1, then the DFE is globally asymptotically stable in Ω, and thus the vector-borne disease always dies out; if R 0 > 1, then the unique endemic equilibrium is globally asymptotically stable in Ω, so that the disease, if initially present, will persist at the unique endemic equilibrium level.
Assuming that parameters are fixed, the threshold R 0 is influenceable by the control values. Figure 6 .2 gives this relationship. It is possible to realize that the control c m is the one that most influences the basic reproduction number to stay below unit. Besides, the control in the aquatic phase alone is not enough to maintain R 0 below unit: an application close to 100% is required.
The simulations were carried out using the following numerical values: N h = 480000, B = 0.8, β mh = 0.375, β hm = 0.375, µ h = 1/(71 × 365), η h = 1/3, µ m = 1/10, ϕ = 6, µ A = 1/4, η A = 0.08, m = 3, k = 3. The initial conditions for the problem were: S h0 = N h − 10, I h0 = 10, R h0 = 0, A m0 = kN h , S m0 = mN h , I m0 = 0. With these values, one has M > 0. As in the previous chapter the values related to humans describe the reality of Cape Verde [32] and the information about mosquitoes is based on Brazil [30, 47] .
All computational calculus consider one year for time interval. Although the final time was t f = 365 days, the figures show graphics in suitable windows, in order to provide a better analysis. All the simulations and graphics were done in Matlab. To solve the differential equation system, the ode45 routine was used. This function implements a Runge-Kutta method with a variable time step for efficient computation (see [110] for more details about the code). The number of infected humans from the model (6.1)-(6.2) is higher when compared with what really happened in Cape Verde. As far as it was possible to investigate in the local news, the government of Cape Verde has done their best to banish the mosquito, with media campaigns appealing to people to remove or cover all recipients that could serve to breed the mosquito, and to use insecticide in critical areas. However, it was not possible to quantify those efforts in precise terms.
Next, follows a set of simulations using different controls. In each figure, only one control is used, continuously, which means that the others are not applied. The aim of this simulation is to see the importance of the control and what repercussions has on the model. Figures 6.5 and 6.6 concern the adulticide control, Figures 6.7 and 6.8 the larvicide control, and Figures 6.9 and 6.10 the mechanical control. Using a small quantity of each control, the number of infected people falls dramatically. In some cases, in spite of all graphs displaying five simulations, some curves are so close to zero that it is difficult to distinguish them. Figures 6.5 and 6.6 show that excellent results for the human population are obtained by covering only 25% of the country with insecticide for adult mosquitoes. The simulations were done considering that the aedes aegypti does not become resistant to the insecticide and that it is financially possible to apply insecticide during all time. Figures 6.7-6.10 are related to the applied controls in the aquatic phase of the mosquito. In these graphics the controls were studied separately, but one is closely related to the other. The application of these controls are not sufficient to decrease the infected human to zero, but the removal In the next section, using OC strategy we will find the best solution for the controls.
Epidemiological models may give some basic guidelines for public health practitioners, comparing the effectiveness of different potential management strategies.
In reality, a range of constraints and trade-offs may substantially influence the choice of practical strategy, and therefore their inclusion in any modelling analysis may be important. Frequently, epidemiological models need to be coupled to economic considerations, such that control strategies can be judged through holistic cost-benefit analysis. Control of livestock disease is a scenario when cost-benefit analysis can play a vital role in choosing between cheap, weak controls that lead to a prolonged epidemic, or expensive but more effective controls that lead to a shorter outbreak.
Normalizing the previous ODE system (6.1)-(6.2), we obtain:
with the initial conditions
The cost functional considered was
where γ D , γ S , γ L and γ E are weights related to costs of the disease, adulticide, larvicide and mechanical control, respectively. In a first approach to this problem, it is assumed that all weights are the same, which means γ D = γ S = γ L = γ E = 0.25 (Case A).
The OC problem was solved using two different packages: DOTcvp [67] and Muscod-II [79] . The mathematical formulation of the SIR+ASI problem, for the both packages, is available on [110] . The simulation behavior is similar, and we decided to show only the DOTcvp results. The optimal functions for the controls are given in Figure 6 .13.
The adulticide was the control that more influenced the decreasing of that ratio and, as consequence, the decreasing of the number of infected people and mosquitoes, matching with the results obtained for the basic reproduction number in Section 6.1. Therefore, the adulticide was the most used. We believe that the other controls do not assume an important role in the epidemic episode, because all the events happened at a short period of time, which means that adulticide has more impact. However the mosquito control in the aquatic phase can not be neglected. In situations of longer epidemic episodes or even in an endemic situation, the larval control represents an important tool. Figure 6 .14 presents the number of infected humans. Comparing the optimal control case with a situation with no control, the number of infected people decreased considerably. Besides, in the situation where OC is used, the peak of infected people is minor, which facilitates the work in health centers, because they can provide a better medical monitoring.
A second analysis was made, taking into account different weights on the functional. Table 6 .1 resumes the weights chosen for perspectives: not only economic issues (cost of insecticides and educational campaigns), but also human issues are considered. In case A, all costs were equal. In case B is given more impact on the infected people, considering that the treatment and absenteeism to work is very prejudicial to Figure 6 .14: Comparison of infected individuals under an optimal control situation and no controls the country, when compared with the cost of insecticides and educational campaigns. In case C, the costs with killing mosquitoes and educational campaigns have more impact in the economy. Higher total costs were obtained when human life had more weight than controls measures as can be checked in Table 6 Table 6 .1: Different weights for the functional and respective values Figure 6 .15 shows the number of infected human in each bioeconomic perspective. We can realize that Case A and Case C are similar. It can be explained by the low weight given to the cost of treatment (cases A and C) when compared with the heavy weight given in case B. Figure 6 .16 presents the behavior of the controls for the A, B and C cases. Again, as adulticide is the one that has more influence on the model, this is the control that most varies when the weights are changed. A third analysis was made to the model: it changed the functional in order to study the effects of each control when separately considered. Therefore, the new functional also considers bioeconomic perspectives, but only includes two variables: the costs with infected humans (with γ D = 0.5) and the costs with only one control (with γ i = 0.5, i ∈ {S, L, E}). Thus, in Figure 6 .17 are presented the proportion of adulticide (a) and infected humans (b), when the functional is 
Dengue disease breeds, even in the absence of fatal forms, significant economic and social costs: absenteeism, debilitation and medication. To observe and to act at the epidemics onset could save lives and resources to governments. Moreover, the under-reporting of Dengue cases is probably the most important barrier in obtaining an accurate assessment. A compartmental epidemiological model for Dengue disease composed by a set of differential equations was presented. Simulations based on clean-up campaigns to remove the vector breeding sites, and also simulations on the application of insecticides (larvicide and adulticide), were made. It was shown that even with a low, although continuous, index of control over time, the results are surprisingly positive. The adulticide was the most effective control, since with a low percentage of insecticide, the basic reproduction number is kept below unit and the number infected humans was smaller. However, to bet only in adulticide is a risky decision. In some countries, such as Mexico and Brazil, the prolonged use of adulticides has been increasing the mosquito tolerance capacity to the product or even they become completely resistant. In countries where Dengue is a permanent threat, governments must act with differentiated tools. It will be interesting to analyze these controls in an endemic region and with several outbreaks. We believe that the results will be quite different. Aedes aegypti eradication is not considered to be feasible and, from the environmental point of view, not desirable. The aim is to reduce the mosquito density and, simultaneously, amount the level of immunity on the human population. The increase of population herd immunity can be reached by two ways: increasing the resistant people to the disease implying an increasing of infected individuals or with a vaccination campaign. No commercially available clinical cure or vaccine is currently available for Dengue, but efforts are underway to develop one [10, 68] . The potential of prevention of Dengue by immunization seems to be technically feasible and progress is being made in the development of vaccines that may protect against all four Dengue viruses. In the next chapter, a model of this disease with a vaccine simulation as a new strategy to fight the disease will be analyzed.
This chapter was based on work accepted in the peer reviewed journal [115] and the peer reviewed conference proceedings [116] .
As the development of a Dengue vaccine is ongoing, it is simulated an hypothetical vaccine as an extra protection to the population. In a first phase, the vaccination process is studied as a new compartment in the model, and some different types of vaccines are simulated: pediatric and random mass vaccines, with distinct levels of efficacy and durability. In a second step, the vaccination is seen as a control variable in the epidemiological process. In both cases, epidemic and endemic scenarios are included in order to analyze distinct outbreak realities.
In 1760, the Swiss mathematician Daniel Bernoulli published a study on the impact of immunization with cowpox upon the expectation of life of the immunized population [119] . The process of protecting individuals from infection by immunization has become a routine, with historical success in reducing both mortality and morbidity.
The impact of vaccination may be regarded not only as an individual protective measure, but also as a collective one. While direct individual protection is the major focus of a mass vaccination program, the effects on population also contribute indirectly to other individual protection through herd immunity, providing protection for unprotected individuals [48] (see scheme in Figure 7 .1). This means that when we have a large neighborhood of vaccinated people, a susceptible individual has a lower probability in coming into contact with the infection, being more difficult for diseases to spread, which decreases the relief of health facilities and can break the chain of infection.
Vector control remains the only available strategy against Dengue. Despite integrated vector control with community participation, along with active disease surveillance and insecticides, there are only a [23] . Besides, the levels of resistance of Aedes aegypti to insecticides has increased, which implies shorter intervals between treatments, and only few insecticide products are available in the market due to the high costs for development and registration and low returns [76] .
Dengue vaccines have been under development since the 1940s, but due to the limited appreciation of global disease burden and the potential markets for Dengue vaccines, industry interest languished throughout the 20th century. However, in recent years, the development of Dengue vaccines has accelerated dramatically with the increase in Dengue infections, as well as the prevalence of all four circulating serotypes. Faster development of a vaccine became a serious concern [94] .
Economic analysis are now conducted periodically to guide public support for vaccine development in both industrialized and developing countries, including a previous cost-effectiveness study of Dengue [28, 123, 126] . The authors compared the cost of the disease burden with the possibility of making a vaccination campaign; they suggest that there is a potential economic benefit associated with promising Dengue interventions, such as Dengue vaccines and vector control innovations, when compared to the cost associated to the disease treatments.
Constructing a successful vaccine for Dengue has been challenging: the knowledge of disease pathogenesis is insufficient and in addition the vaccine must protect simultaneously against all serotypes in order to not increase the level of DHF [128] .
Nevertheless, several promising approaches are being investigated in both academic and industrial laboratories. Vaccine candidates include live attenuated vaccines obtained via cell passages or by recombinant DNA technology (such as those being developed by the US National Institutes of Allergy and Infectious Diseases, InViragen, Walter Reed Army Institute of Research/GlaxoSmithKline, and Sanofi Pasteur) and subunit vaccines (such as those developed by Merck/Hawaii Biotech) [60, 138] . Recent studies indicate that, by the progress in clinical development of Sanofi Pasteur's live attenuated tetravalent chimeric vaccine, a vaccine could be licensed as early as 2014 [87] . The team is carrying out an efficacy study on a vaccine covering four serotypes on 4000 children aged four to eleven years old in Muang district, Thailand.
At this time, the features of Dengue vaccine are mostly unknown. So, in this chapter we opt to present a set of simulations with different types of vaccines and we have explored also the vaccination process under two different perspectives. The first one uses a new compartment in the model and several kinds of vaccination are considered. A second perspective is studied using the vaccination process as a disease control in the mathematical formulation. In this case the theory of OC is applied. Both methods assume a continuous strategy vaccination.
In this section, a new compartment V is added to the previous SIR model related to the human population. This new compartment represents the new group of human population that is vaccinated, in order to distinguish the resistance obtained through vaccination and the one achieved by disease recovery.
Two forms of random vaccination are possible: the most common for human disease is pediatric vaccination to reduce the prevalence of an endemic disease; the alternative is random vaccination of the entire population in an outbreak situation. In both types, the vaccination can be considered perfect conferring 100% protection for all life or else imperfect. This last case can be due to the difficulty of producing an effective vaccine, the heterogeneity of the population or even the life span of the vaccine.
For many potentially human infections, such as measles, mumps, rubella, whooping cough, polio, there has been much focus on vaccinating newborns or very young infants. Dengue can be a serious candidate for this type of vaccination.
In the SV IR model, a continuous vaccination strategy is considered, where a proportion of the newborn p (where 0 ≤ p ≤ 1), was by default vaccinated. This model also assumes that the permanent immunity acquired through vaccination is the same as the natural immunity obtained from infected individuals eliminating the disease naturally. The population remains constant, i.e., N h = S h + V h + I h + R h . The new model for human population is represented in Figure 7 
We are assuming that it is a perfect vaccine, which means that it confers life-long protection. As a first step, it is necessary to determine the basic reproduction number without vaccination (p = 0).
Theorem 7.2.1. The basic reproduction number, R 0 , associated to the differential system (7.1) without vaccination is given by
Proof. The proof of this theorem is similar to the one in the previous chapter (see proof of Theorem 6.1.2). Just consider
In this chapter, we make all the simulations in two scenarios: an epidemic and an endemic situation (programming codes available in [110] ). For these, the following parameter values of the differential system and initial conditions were used (Tables 7.1 and 7.2).
There were two main differences between the epidemic episode and an endemic situation. Firsty, in the endemic situation there was a slight decrease in the average daily biting B and transmission probabilities β mh and β hm , that could be explained by the fact that the mosquito could have more difficulties to find a naive individual. The second difference is concerned with the strong increase of the initial human population that is resistant to the disease. This may be explained by the fact that the disease, in an endemic situation, already creates an immune resistance to the infection, i.e., the population already has herd immunity. With these values we obtain approximately R 0 = 2.46 and R 0 = 1.29 for epidemic and endemic scenarios, respectively.
During an outbreak, the disease transmission assumes different behaviors, according to the distinct scenarios, as can been seen in Figure 7 .3. In one year, the peak in an epidemic situation could reach more than 80000 cases. In contrast, instead in the endemic situation the curve of infected individuals has a more smooth behavior and reaches a peak less than 3000 cases. Figure 7 .4 relates to the mosquito population. In the endemic scenario, because a substantial part of the human population is resistant to the disease, the infected mosquitoes bite a considerable percentage of resistant host, and as consequence, the disease is not transmitted. Suppose that at time t = 0 a proportion p of newborns are vaccinated with a perfect vaccine that causes no side effects. Since this proportion, p, is now immune, R 0 is reduced, creating a new basic reproduction number.
Definition 19 (Basic reproduction number with pediatric vaccination). The basic reproduction number with pediatric vaccination, R p 0 , associated to the differential system (7.1) is given by
where R 0 is defined in (7.2).
Observe that R p 0 ≤ R 0 . Equality is only achieved when p = 0, i.e., when there is no vaccination. The constraint R p 0 < 1 implicitly defines a critical vaccination portion p > p c that must be achieved for eradication, where
Since vaccination entails costs, to choose the smallest coverage that achieves eradication is the best option. This way, the entire population does not need to be vaccinated in order to eradicate the disease. This phenomenon is called herd immunity.
Vaccinating at the critical level, p c , does not instantly lead to disease eradication. The immunity level within the population requires time to build up and at the critical level it may take a few generations before the required herd immunity is achieved. Thus, from a public health perspective, p c acts as a lower bound on what should be achieved, with higher levels of vaccination leading to a more rapid elimination of the disease. Figure 7 .5 shows the simulations related to the proportion of newborns vaccinated (p = 0, 0.25, 0.50, 0.75, 1) in both scenarios. Notice that at time t = 0 no person was vaccinated. In the epidemic situation, as the outbreak reached a peak at the beginning of the year, the proportion of newborns vaccinated at that time is minimum and cannot influence the curve of infected individuals, giving the optical illusion of a single curve. On the other hand, in the endemic case, as the outbreak occurs later, the vaccination campaign starts to produce effects, decreasing the total number of sick humans. This last graphic illustrates that a vaccination campaign centered in newborns is a bet for the future of a country, but does not produce instantly results to fight the disease. To produce immediate results, it is necessary to use random mass vaccination, which means that it is necessary to vaccine a significant part of the population.
A mass vaccination program may be initiated whenever there is an increase of the risk of an epidemic. In such situations, there is a competition between the exponential increase of the epidemic and the logistical constraints upon mass vaccination. For most human diseases it is possible, and more efficient, to not vaccinate those individuals who have recovered from the disease because they are already protected. Another situation could be the introduction of a new vaccine in a population that lives an endemic situation.
Let us consider the control technique of constant vaccination of susceptibles. In this scheme a fraction 0 ≤ ψ ≤ 1 of the entire susceptible population, not just newborns, is being continuously vaccinated. It is assumed that the permanent immunity acquired by vaccination is the same as natural immunity obtained from infected individuals in recovery. The epidemiological scheme is presented in Figure 7 The mathematical formulation for human population (the differential equations related to the mosquito remain equal to the previous subsection) is given by:
For this model, we define a new basic reproduction number.
Definition 20 (Basic reproduction number with random mass vaccination [145] ). The basic reproduction number with random mass vaccination, R ψ 0 , associated to the differential system (7.3) is
where R 0 is defined in (7.2) .
Comparing this model with the constant vaccination of newborns model, it is apparent that instead of constantly vaccinating a portion of newborns, a part of the entire susceptible population is now being continuously vaccinated. Since the natural birth rate µ h is usually small, the fraction pµ h of newborns being continuously vaccinated will be small, whereas in this model, will be also a larger group of susceptibles can be continuously vaccinated, ψS h . Due to this, we expect that this model should require a smaller proportion of ψ to achieve eradication.
Notice that R ψ 0 ≤ R 0 . Equality is only achieved in the limit ψ = 0, that is, when there is no vaccination. The constraint R ψ 0 < 1 implicitly defines a critical vaccination portion ψ > ψ c that must be achieved for eradication, where ψ c = (R 0 − 1) µ h . Observe that in spite of the calculations done in the period of 365 days, the figures only show suitable windows, in order to provide a better analysis. In both scenarios, even with a small coverage of population, vaccination dramatically decreases the number of infected. The epidemic scenario has changed from less than 80000 cases (with no vaccination, Figure 7 .3) to less than 1200 cases vaccinating only 5% of population. In the endemic scenario, the decrease is more accentuated.
Until here, we have considered a perfect vaccine, which means that every vaccinated individual remains resistant to the disease. However, a majority of the available vaccines for the human population does not produce 100% success in the disease battle. Usually, the vaccines are imperfect, which means that a minor percentage of cases, in spite of vaccination, are infected.
Most of the theory on the disease evolution is based on the assumption that the host population is homogeneous. Individual hosts, however, may differ and they may constitute very different habitats. In particular, some habitats may provide more resources or be more vulnerable to virus exploitation [56] . The use of models with imperfect vaccines can describe better this type of human heterogeneity. Another explanation for the use of imperfect vaccines is that until now we had considered models that assumed that as soon as individuals begin the vaccination process, they become immediately immune to the disease. However, the time it takes for individuals to obtain immunity by completing a vaccination process cannot be ignored, because meanwhile an individual can be infected.
In this section a continuous vaccination strategy is considered, where a fraction ψ of the susceptible class was vaccinated. The vaccination may reduce but not completely eliminate susceptibility to infection. For this reason we consider a factor σ as the infection rate of vaccinated members. When σ = 0 the vaccine is perfectly effective and when σ = 1 the vaccine has no effect at all. The value 1 − σ can be understood as the efficacy level of the vaccine.
The new model for the human population is represented in Figure 7 .8. Figure 7 .8: Epidemiological SV IR model for human population with an imperfect vaccine Therefore, the differential system is as follows:
(7.5)
For this system of differential equations, we have a new basic reproduction number. [85] ). The basic reproduction number with an imperfect vaccine, R σ 0 , associated to the differential system (7.5) is
where R 0 is defined in (7.2) .
Notice that R ψ 0 ≤ R σ 0 and when the vaccine is perfect (σ = 0), R σ 0 degenerates into R ψ 0 . In other words, a high efficacy vaccine leads to a lower vaccination coverage to eradicate the disease. However, it is noted in [85] that it is much more difficult to increase the efficacy level of the vaccine when compared to controlling the vaccination rate ψ. Figure 7 .7a, in the epidemic scenario with a perfect vaccine, the number of human infected has reached to a maximum peak of 1200 cases per day, in the worst scenario (ψ = 0.05).
Using an imperfect vaccine, with a level of efficacy of 80% (Figure 7 .9a), with the same values for ψ, the maximum peak increases until 9000 cases.
We conclude that production of a vaccine with a high level of efficacy has a preponderant role in the reduction of the disease spread. Figures 7.9c and 7 .9d reinforce the previous sentence. Assuming that 85% of the population is vaccinated, the numbers of infected cases decreases sharply with the increasing of the effectiveness level of the vaccine.
According to [34] , an acceptable level of efficacy is at least 80% against all four serotypes, and 3 to 5 years for the length of protection. These are commonly considered across countries, as the minimum acceptable levels.
In the next subsection we study another type of imperfect vaccine: a vaccine that confers a limited life-long protection.
Until the 1990s, this was an universal assumption of mathematical models of vaccination: there is no waning of vaccine-induced immunity. This assumption was routinely made because, for most of the major vaccines against childhood infectious diseases, it is approximately correct [119] .
Suppose that the immunity, obtained by the vaccination process, is temporary. Assume that immunity has the waning rate θ. Then the model for humans is given by Definition 22 (Basic reproduction number with an imperfect vaccine). The basic reproduction number with an imperfect vaccine, R θ 0 , associated to the differential system (7.6) is
where R ψ 0 is defined in (7.4 ).
According to [125] , the basic reproduction numbers are the same, R θ 0 and R ψ 0 , because the disease will still spread at the same rate with or without temporary immunity. However, we should expect that the convergence rate will be different between the random mass vaccination and random mass vaccination with waning immunity, since the disease will be eradicated faster in the constant treatment model without waning immunity compared to the other with waning immunity. Depending on the vaccine that will be available on the market, it will be possible to choose or even combine features. In the next section we will define the vaccination process as a control system.
In this section we consider a SIR model for humans and an ASI model for mosquitoes. The parameters remain the same as in the previous chapter.
The vaccination is seen as a control variable to reduce or even eradicate the disease. Let u be the control variable related to the proportion of susceptible humans that are vaccinated. A random mass vaccination with waning immunity is selected. In this way, a parameter θ associated to the control u represents the waning immunity process. Figure 7 .12 shows the epidemiological scheme for the human population. Figure 7 .12: Epidemiological SIR model for the human population using the vaccine as a control The model is described by an initial value problem with a system of six differential equations:
The main aim is to study the optimal vaccination strategy, considering both the costs of treatment of infected individuals and the costs of vaccination. The objective is to
where γ D and γ V are positive constants representing the weights of the costs of treatment of infected and vaccination, respectively. Using OC theory is possible to solve the problem.
Let us consider the following set of admissible control functions: ∆ = {u(·) ∈ (L ∞ (0, t f ))|0 ≤ u(t) ≤ 1, ∀t ∈ [0, t f ]}. Table (7. 2), admits a unique optimal solution (S * h (·), I * h (·), R * h (·), A * m (·), S * m (·), I * m (·)) associated with an optimal control u * (·) on [0, t f ], with a fixed final time t f . Moreover, there exists adjoint functions, λ * i (·), i = 1...6 such that
with the transversality conditionsλ i (t f ) = 0, i=1, . . . 6. Furthermore, u * = min 1, max 0,
Proof. The existence of optimal solutions (S * h (·), I * h (·), R * h (·), A * m (·), S * m (·), I * m (·)) associated to the optimal control u * (·) comes from the convexity of the integrand of the cost function (7.8) with respect to the control u and the Lipschitz property of the state system with respect to state variables (S h , I h , R h , A m , S m , I m ) (for more details see [24] ). According to the Pontryagin Maximum Principle [106] , if u * (·) ∈ ∆ is optimal for the problem considered, then there exists a nontrivial absolutely continuous mapping λ : [0, t f ] → R, λ(t) = (λ 1 (t), λ 2 (t), λ 3 (t), λ 4 (t), λ 5 (t), λ 6 (t)), called the adjoint vector, such that (7.11) and holds almost everywhere on [0, t f ]. Moreover, the transversality conditions λ i (t f ) = 0, i = 1, . . . 6, hold. System (7.9) is derived from (7.12), and the optimal control (7.10) comes from the minimality condition (7.13).
The simulations were carried out using the values of the previous section. The system was normalized, using the same strategy as in Chapter 5. It was considered that the waning immunity was at a rate of Method Epidemic scenario Endemic scenario Direct (DOTcvp) 0.07505791 0.00189056 Indirect (backward-forward) 0.06070556 0.00080618 Table 7 .3: Optimal values of the cost functional (7.8) θ = 0.05. The OC problem was solved using two methods: direct [9, 133] and indirect [81] . The direct method uses the cost functional (7.8) and the state system (7.7) and was solved by DOTcvp [67] . The indirect method used is an iterative method with a Runge-Kutta scheme, solved through ode45 of MatLab. Figure 7 .13 shows the optimal control obtained by both methods. Notice that DOTcvp only gives the optimal control as a constant piecewise function. Table 7 .3 shows the costs obtained by the two methods in both scenarios. The indirect method gives a lower cost. This method uses more mathematical theory about the problem, such as the adjoint system (7.11) and optimal control expression (7.10). Therefore it makes sense that the indirect method produces a better solution. Using the optimal solution as reference, some tests were performed, regarding infected individuals and costs, when no control (u ≡ 0) or upper control (u ≡ 1) is applied. Table 7 .4 shows the results for DOTcvp in the three situations. In both scenarios, using the optimal strategy of vaccination produces better costs with the disease, when compared to not doing anything. Once there is no control, the number of infected humans is higher and produces a more expensive cost functional. Figure 7 .14 shows the number of infected humans when different controls are considered. It is possible to see that using the upper control, which means that everyone is vaccinated, implies that just a few individuals were infected, allowing eradication of the disease. Although the optimal control, in the sense of objective (7.8), allows the occurrence of an outbreak, the number of infected individuals is much lower when compared with a situation where no one is vaccinated.
We conclude that a vaccination campaign in the susceptible population, and assuming a considerable efficacy level of the vaccine, can quickly decrease the number of infected people. 
The worldwide expansion of the Dengue fever is a growing health problem. Dengue vaccine is an urgent challenge that needs to be overcome. This may be commercially available within a few years, when the researchers find a formula that protect against all four Dengue viruses. A vaccination program is seen as an important measure used in infectious disease control and immunization and eradication programs.
In the first part of the chapter, different types of vaccine, as well as their features and some of coverage thresholds, were introduced. The main idea was to study some types of vaccines in order to cover most of the future vaccine features. The main goal of a vaccination program is to reduce the prevalence of an infectious disease and ultimately to eradicate it. It was shown that eradication success depends on the type of vaccine as well as on the vaccination coverage. The imperfect vaccines may not completely prevent infection but could reduce the probability of being infected, thereby reducing the disease burden. In this study all the simulations were done using epidemic and endemic scenarios to illustrate distinct realities.
A second analysis was made, using an OC approach. The vaccine behaved as a new disease control variable and, when available, can be a promising strategy to fight the disease.
Dengue is an infectious tropical disease difficult to prevent and manage. Researchers agree that the development of a vaccine for Dengue is a question of high priority. In the present study we have shown how a vaccine results in saving lives and at the same time in a reduction of the budget related with the disease. As future work we intend to study the interaction of a Dengue vaccine with other kinds of control already investigated in the literature, such as insecticide and educational campaigns [113, 118] .
This chapter was based on work accepted in the peer reviewed proceedings [114] .
Mathematical models can be a powerful tool to understand epidemiological phenomena. These models can be used to compare, plan, implement and evaluate several programs related to detection, prevention and control of infectious diseases. Indeed, one of the most important issues in epidemiology is to improve control strategies with the final goal to reduce or even eradicate the disease. In this thesis were constructed and analyzed some models for the spreading of diseases, particularly for Dengue Fever.
The links between health and sustainable development are illustrated for this disease. Any attempt in predicting and preventing the disasters caused by the disease will imply a global strategy that takes into account environmental conditions, levels of poverty and illiteracy and, eventually, degree of coverage by vaccination programs [35] .
Although vector control strategies were already available before the Second World War, Dengue pandemic was underestimated. It became a global public health problem in the past 60 years and a major concern for WHO.
The main contributions of this thesis can be classified into three main categories, namely, model formulation, mathematical and computational analysis, and contributions to public health/intervention design.
Deterministic models for assessing the combined impact of several control measures in Dengue disease were considered.
In Chapter 4 an old OC problem for Dengue was revisited. It was shown that new and robust tools bring refreshing solutions to the problem. Some analysis with discretization schemes were carried out in order to understand the best way to implement the direct methods in future approaches. For this problem, is was better to use robust solvers than to implement higher order discretization methods, due to the increasing of problem's dimension.
In Chapter 5, a SEIR+ASEI model was studied. Threshold criteria was established ensuring the disease eradication and hence convergence to the so-called disease free solution. Using real data from Cape Verde, the application of insecticide in the country during the outbreak was simulated and its repercussions were analyzed. The study of this outbreak was performed in several phases: firstly, using a previously calculated constant control for the whole period, with the aim of having a basic reproduction number below unity; then, several periodic strategies were studied in order to find the best logistics approaches to implement that keep R 0 less than one; finally, an OC approach to compare with the previous suboptimal approaches was used. For this final phase, several perspectives of the problem were analyzed, including bioeconomic, medical and economic approaches. Depending of the main target to achieve, the results for the control and infected individuals vary.
In Chapter 6, a SIR+ASI model for Dengue was presented. The lost of two differential equations, was compensated by the introduction of two more controls: in addition to adulticide, larvicide and mechanical control were introduced. Similarly to Chapter 5, a threshold criteria for the eradication of the disease was established. The influence of the controls in this threshold was analyzed. Then, varying the controls, separately and simultaneously, an analysis of the importance/consequence of each control in the development of the disease was made. Finally, an OC approach using different weights for the variables in the functional was studied, in order to establish the best optimal curve for each control and the respective effects in the development/erradication of the disease.
In Chapter 7, some simulations with different types of vaccines were made. As a Dengue vaccine is not yet available, distinct hypothetical models to introduce the vaccine were studied. In Section 7.2, was considered a new compartment V h producing a new model SVIR+ASI. This research comprised SVIR models with perfect vaccines -constant vaccination of newborns and constant vaccination of susceptibles -and imperfect vaccines -using a level of efficacy below 100% and also with waning immunity. In Section 7.3, the vaccination process was studied as a control of the epidemic model. For this, as in the previous chapter, an OC approach was presented. All the simulations were done in epidemic and endemic scenarios, in order to understand what type of repercussions could bring each kind of vaccines.
In this work, there was a concern in producing a mathematical analysis for all new models presented. Epidemiological concepts, such as the basic reproduction number and equilibrium points, were calculated. The equilibrium points were classified and their local stability analyzed. The OC theory was used in order to provide the best strategies for each model, involving direct and indirect approaches.
Throughout the thesis, a set of software packages were used, showing the importance of the these tools in the development of some mathematical fields. To solve ODE systems, codes in Matlab and Scilab were implemented. To calculate the equilibrium points and the basic reproduction number, Mathematica and Maple were used. For the OC approach, several packages were also selected: from programmed codes in AMPL and run in NEOS Server (Ipopt, Snopt, Knitro, Muscod-II) to Fortran codes in the Linux environment (OC-ODE) or even programs coded in Matlab (DOTcvp and indirect methods). The software choice was varying during the research process, due to availability, robustness and solving speed. The main codes developed in this thesis are available at: https://sites.google.com/site/hsofiarodrigues/home/ phd-codes-1 [110] . Using direct and indirect methods, the solutions obtained were similar, reinforcing the confidence in the results.
The study provides some important epidemiological insights into the impact of vector control measures. Dengue burden decreases with the increasing of vector control measures (adulticide, larvicide and mechanical control). Furthermore, the adulticide should be the first/main measure to apply when an outbreak occurs, whereas the other two measures should be considered as a long time prevention.
The last control measure to be studied was the vaccination. It was shown that the vaccine, when available, could bring advantages not only in the reduction of infected individuals, but also decreasing the disease costs.
A PhD thesis is always an unfinished process, but some day it is necessary to stop. Therefore, some topics were not explored and can be understood as future directions of the work.
A first suggestion is to use heterogeneity for both populations, dividing each one in more compartments [84] . Human population is not immunologically homogeneous, presenting groups with distinct level of risk, related with age, sex or the presence of a disease or immunosuppressive drugs which would be the case with transplant patients or cancer suffers. For example, the transmission probability in young children is higher and generally with more severe symptoms when compared to adult transmission. Moreover, the mosquito population does not have the same behavior during the whole year, depending specially on weather conditions. Temperature and humidity are key variables on vector population dynamics. It will be interesting to add some seasonality factors into the last models [45, 101] .
Another open question is the introduction of immigration and tourism issues. Along the work it was considered a constant population, but the addition of new individuals could induce new outbreaks.
Other aspect is related to the disease development in the presence of several serotypes. While in Cape Verde only one serotype was found, the interaction of several serotypes in Asia is already a reality. This will induce changes in the model, not only increasing of the number of variables but also causing a more expressive number of DHF cases [105] .
Currently, Portuguese researchers are developing a new repellent for mosquito that could be considered as a new control for the disease. This product does not kill the mosquito, but prevents the bite by deviating it from the target and consequently decreasing the chain of the disease transmission.
With the viability of the vaccine, it will be possible to fit a better model according to the vaccine used. One of the possibilities is using pulse vaccination, where children at a certain age cohorts are periodically immunized [41] . The theoretical challenge of pulse vaccination is the a priori determination of the pulse interval for specific values of R 0 , the proportion of vaccinated p and the population birth rate µ.
Most of the analysis and models presented in this thesis can be adapted for other vector-borne diseases -such as malaria, yellow fever, West Nile virus, chikungunya, Japanese encephalitis -by just fitting some variables/parameters and some initial assumptions intrinsic to the disease [18, 46, 72] . When attempting to model epidemics and control for public health applications, there is the compelling urge to make models as sophisticated as possible, including many details about the host and the vector. Although this strategy may be useful when such details are known or exist suitable data, it may lead to a false sense of accuracy when reliable information is not available. Another approach is to keep the model simple and, instead of using the conventional differential calculus, to apply fractional calculus to fit to the disease reality [107] or to use the general theory of time scales [15] .
"If people do not believe that mathematics is simple, it is only because they do not realize how complicated life is." -John Louis von Neumann, 1947