Optimal Intervention Strategies for a SEIR Control Model of Ebola Epidemics

A SEIR control model describing the Ebola epidemic in a population of a constant size is considered over a given time interval. It contains two intervention control functions reflecting efforts to protect susceptible individuals from infected and exposed individuals. For this model, the problem of minimizing the weighted sum of total fractions of infected and exposed individuals and total costs of intervention control constraints at a given time interval is stated. For the analysis of the corresponding optimal controls, the Pontryagin maximum principle is used. According to it, these controls are bang-bang, and are determined using the same switching function. A linear non-autonomous system of differential equations, to which this function satisfies together with its corresponding auxiliary functions, is found. In order to estimate the number of zeroes of the switching function, the matrix of the linear non-autonomous system is transformed to an upper triangular form on the entire time interval and the generalized Rolle’s theorem is applied to the converted system of differential equations. It is found that the optimal controls of the original problem have at most two switchings. This fact allows the reduction of the original complex optimal control problem to the solution of a much simpler problem of conditional minimization of a function of two variables. Results of the numerical solution to this problem and their detailed analysis are provided. Mathematics 2015, 3 962


Introduction
Ebola is a lethal virus for humans.It was previously confined to Central Africa, but recently was also identified in West Africa ( [1]).As of October 8, 2014, the World Health Organization (WHO) reported 4656 cases of Ebola virus deaths, with most cases occurring in Liberia ( [2]).The extremely rapid increase of the disease, high mortality rate, and the nonexistence of a vaccine make this virus a major problem for public health.Ebola is transmitted through direct contact with blood, bodily secretions and tissues of infected humans and primates ( [3,4]).Therefore, the practical control interventions of public health measures are major factors for stopping Ebola transmission.To understand the dynamics of the spread of Ebola infection in the affected countries, it is crucial to have models of this process to simulate.
In this paper we focus our attention on a SEIR deterministic model and its use for the description of the Ebola epidemics ( [5][6][7]).In order to deal with the epidemics of Ebola, governments of affected countries decided to implement intervention measures: some applied social distancing, early detection of infectious individuals, quarantine procedures, campaigns for information and education, and acceleration of the burial process.In this context, we introduce into the SEIR model the control functions reflecting these intervention measures and consider the optimal control problem to study the effect of these intervention control strategies on the epidemic spread.Optimality is measured as the minimization of the weighted sum of total fractions of infected end exposed individuals and the total cost of intervention control constraints over a given time interval.
We model the effects from intervention control measures by reducing the rate at which the disease is contracted from an average individual during such measures (called shortly transmission rate).We justify this as follows: suppose during an outbreak one starts the measures orienting susceptible individuals to avoid contracting the virus (assuming some protective behavior, for example, washing hands, avoiding close environments, etc.).The effect of these measures will be that the probability of a susceptible individual contracting the virus will decrease.This reduction (increase) is bounded below (above) and costs of such measures are linear on the controls.Hence, the original problem renders itself to analytical treatment and we can establish the main properties about optimal strategies of intervention control measures.It simplifies the problem and allows its complete numerical study.Influence of health-promotion and educational campaigns on the spread of epidemic, which is described by various SIR models, was previously studied in [8,9].
The paper is organized as follows.In Section 2, we introduce a SEIR control model, describe its variables, control functions, parameters and study the properties of these variables.In Section 3 we state the optimal control problem and consider in detail the epidemiological significance of the introduced functional.Next, we discuss the existence of the optimal solutions in our problem.In Section 4, for the analysis of these solutions, we use the Pontryagin maximum principle to write the Hamiltonian of the optimal control problem and corresponding adjoint system and find the formulas which determine the appropriate optimal controls through solutions of this system.Based on these formulas we find the switching function that defines the type of the optimal controls.Finally, we obtain the system of differential equations to which this function satisfies with its corresponding auxiliary functions.Section 5 is devoted to studying the properties of the switching function, which lead to the conclusion regarding bang-bang type of the optimal controls and help in the estimation of the number of switchings.These results are the basis of the method of the solution to the optimal control problem, a description of which is given in Section 6.Here the Lipschitz constant is found for use in the numerical algorithm.Results of numerical calculations and their detailed analysis are presented in Section 7. Finally, Section 8 contains our conclusions.

SEIR Model and Its Properties
At a given time interval [0, T ] we consider a SEIR model described by the following system of differential equations: (1) Such a model describes the spread of the epidemic in a population of the constant size N .Indeed, considering that the equality holds, we add together the equations of system (1).Taking into account relationship (2), we find the equality: System (1) was used in [5] to describe the Ebola outbreak in Zaire (1976), then, in [6], for studying the effects of public health measures during Ebola outbreaks in Congo (1995) and Uganda (2000).Finally, this system was applied in [7] for the study of recent Ebola outbreaks in Guinea, Liberia and Sierra Leone (2014).
In system (1) the quantity S(t) is the number of susceptible individuals.The quantity E(t) denotes the number of individuals exposed to the virus of infected but not yet infectious.Individuals that are infected with the disease and are suffering the symptoms of Ebola are classified as infectious individuals, the number of them is denoted by I(t).Similarly, the number of deceased or recovered individuals is denoted by R(t).In system (1) the value βN −1 I(t)S(t) is the number of individuals infected due to direct contact with an infected individual, and the value αN −1 E(t)S(t) is the number of individuals infected due to direct contact with an exposed individual.Here β is the transmission rate; α = qβ, where q ∈ [0, 1] is a weight factor taking into account that a susceptible individual has a higher chance of getting infected from an infectious individual than from an exposed individual.The value δE(t) is the individuals in the exposed stage, which show symptoms of the disease and pass on to the infectious stage; δ is the infectious rate.The value γI(t) is the individuals in the infectious stage, which die; γ is the death rate.We consider that death and recovery are taken to be the same, since there has not been a case in which a person who survived Ebola contracts the disease again.Now, we will make controllable SEIR model, described by system (1).For this, we introduce two control functions u(t), v(t) implying the efforts of preventing susceptible individuals from becoming infectious individuals as a result of contact with infectious and exposed ones.For our controls we have the following constraints: where u max = β, v max = α.We observe that the control functions u(t), v(t) regulate the goals and efforts of two similar interventions.In the case of control u(t), the equality u(t) = β corresponds to not having an intervention affecting the susceptible individuals and the equality u(t) = u min corresponds to the maximum effort that can be made.Similar conclusions are valid for the control v(t).Today there are no licensed treatment and approved vaccines for Ebola; there are examples of experimental vaccines as proof of the possibility of treatment by vaccinating infected individuals (see [10]).Therefore, controlled interventions are the major factor for stopping Ebola transmission.Thus, we have the following SEIR control model of the type: in which the equation for the function R(t) is excluded, and this function can be found from equality (3) by the formula: Moreover, in system (5) we introduce the new variables: with corresponding initial values: for which the following inequalities hold: s 0 , e 0 , i 0 > 0; s 0 + e 0 + i 0 < 1.
These variables are fractions of the quantities S(t), E(t), I(t) in a population of size N .Then, for the variables s(t), e(t), i(t) we obtain the following nonlinear control system: For this system the set of all admissible controls is formed by all possible Lebesgue measurable functions u(t), v(t), which for almost all t ∈ [0, T ] satisfy constraints (4).Now, we introduce a region: Λ = (s, e, i) ∈ 3 : s > 0, e > 0, i > 0, s + e + i < 1 .
Here a sign means transpose.
It is easy to establish the boundedness, positiveness, and continuation of solutions for system (6) by the following lemma.
Lemma 1.If the inclusion (s 0 , e 0 , i 0 ) ∈ Λ holds, then for any admissible controls u(t), v(t) the corresponding solutions s(t), e(t), i(t) for system (6) are defined on the entire interval [0, T ] and satisfy the inclusion: Proof of this fact is standard, and so we omit it.Relationship (7) implies that the region Λ is a positive invariant set for system (6).

Optimal Control Problem for a SEIR Model
For controlled system (6) on the set of all admissible controls (4) we consider the functional: which defines the weighted sum of the total fractions of exposed and infected individuals on the given time interval [0, T ].
Next, as in [9], we introduce the following costs of intervention control constraints: These costs are linear in the controls u(t), v(t) and supposed to be proportional to the fractions of the exposed and infected individuals.If one assumes that these fractions are proportional to the fraction of the number of regions where the disease occurs, to the population size, and therefore to the the fraction of the number of regions to be converted by our intervention controls to the population size as well.The higher the fractions of exposed and infected, the higher corresponding costs.Finally, for system (6) on the set of all admissible controls (4) we consider the optimal control problem of minimization of the sum of values ( 8), (9): The existence in problem (10) of the optimal controls u * (t), v * (t) and corresponding optimal solutions s * (t), e * (t), i * (t) for system (6) follows from Lemma 1 and Theorem 4 ( [11], Chapter 4).

Pontryagin Maximum Principle
In order to analyze problem (10), we apply the Pontryagin maximum principle ( [12]).According to it, we write the Hamiltonian as where ψ 1 , ψ 2 , ψ 3 are the adjoint variables.We calculate the necessary derivatives: Then, by the Pontryagin maximum principle, for the controls u * (t), v * (t) and corresponding solutions s * (t), e * (t), i * (t) there exists a nontrivial solution ψ such that the controls u * (t), v * (t) maximize the Hamiltonian for almost all t ∈ [0, T ], and therefore satisfy the relationships: Here the function is the switching function, which defines the types of the optimal controls u * (t), v * (t) according to formulas (12), (13).
In order to investigate the behavior of the switching function L(t) we will introduce the following auxiliary functions: Now, using the equations of systems ( 6), (11), we obtain the differential equations for the functions L(t), G(t), P (t): Using the initial conditions of system (11), we find the initial conditions for the functions L(t), G(t), P (t): Combining the obtained differential equations and the found initial conditions, we have the system: System ( 14) can be simplified by introducing the following functions: Finally, we obtain the system of differential equations for the functions L(t), G(t), P (t): This system will be the subject of our investigation.

Properties of the Switching Function
Analyzing system (15), we establish the properties of the switching function L(t).The following lemma is valid.
Lemma 2. The switching function L(t) cannot be zero on any subinterval of the interval [0, T ].
Proof of Lemma 2. We assume the contrary.Let the function L(t) become zero everywhere on the interval ∆ ⊂ [0, T ].Then, it is obvious that L(t) = 0 almost everywhere on this interval.Therefore, from the first equation of system (15) we find that G(t) = 0 for all t ∈ ∆.Hence, the derivative Ġ(t) becomes zero almost everywhere on the interval ∆.From the second equation of this system we obtain that P (t) = 0 for t ∈ ∆.Therefore, we have the equality Ṗ (t) = 0 almost everywhere on the interval ∆.Then, the third equation of system ( 15) is also satisfied.Thus, we find the equalities: which are valid for all t ∈ ∆.System ( 15) is the linear non-autonomous homogeneous system of differential equations.Therefore, relationships ( 16) hold on the entire interval [0, T ], and, in particular, at t = T .This fact contradicts the initial conditions of system (15).Our assumption was wrong, and the switching function L(t) does not vanish on any subinterval of the interval [0, T ].
Lemma 2 and formulas ( 12), (13) show that the optimal controls u * (t), v * (t) are bang-bang functions taking values {u min ; u max }, {v min ; v max }, respectively.Moreover, these controls switch from maximum values to minimum values and vice versa at the same moments of switching.
The following important fact is contained in the lemma presented below.
Lemma 3. The switching function L(t) has at most two distinct zeros on the interval [0, T ].
Proof of Lemma 3. Justification of this statement consists of three stages.At the first stage, using special substitutions of variables we reduce the matrix of linear non-autonomous system (15) to the upper triangular form.At first, the following change of the variables is made: where the functions g 1 (t), g 2 (t) will be further defined.Using system (15), we write the system of differential equations for the functions r(t), q(t), µ(t) as We choose the functions g 1 (t), g 2 (t) in such a way that the expressions inside the braces of the last equation of system (17) become zero.Then, this system can be rewritten as and the corresponding differential equations for the functions g 1 (t), g 2 (t) have the following form: Now, a change of variables of the following form in system ( 18) is executed: where the function g 3 (t) will be further defined.Using system (18), we write the system of differential equations for the functions r(t), q(t), µ(t) as We choose the function g 3 (t) in such a way that the expression in the braces of the second equation of system (20) becomes zero.Then, this system has the form: and the corresponding differential equation for the function g 3 (t) is written as follows Thus, using relationships ( 19), ( 22), we have the non-autonomous quadratic system of differential equations for the functions g 1 (t), g 2 (t), g 3 (t): The problem of continuability of solutions of non-autonomous quadratic systems of differential equations was considered in [13,14].In these studies, as well as in the classic textbook on differential equations [15], it was noted that, when the initial conditions g j (0) = g 0 j , j = 1, 3 were added to system (23), an appropriate solution g(t) = (g 1 (t), g 2 (t), g 3 (t)) will be determined on the interval [0, t max ), which is the maximum possible interval for the existence of such solution, and either t max = +∞ or t max < +∞.In the latter case, it may be so that t max ≤ T .Therefore, at the second stage of our arguments we show the existence of such solution g(t) for system (23), that is defined on the entire interval [0, T ].For this, we write the system in a matrix form.Introduce the symmetric matrices: vector functions: and scalar functions: Then, system (23) can be written as follows Next, we assume the contrary, i.e., whatever a solution g(t) for system (24) we do not consider, it is defined on the interval [0, t 1 ), where t 1 ≤ T , and this interval is the maximum possible of the existence of such solution.Then, from Lemma ( [16], § 14, Chapter 4) for solution g(t) it follows the relationship: where || • || is the usual Euclidean norm of a vector.This relationship implies the existence of the number ρ > 0 and the value t 0 ∈ [0, t 1 ) for which the inequality ||g(t)|| ≥ ρ holds for all t ∈ [t 0 , t 1 ).Note that the values ρ and t 0 will be defined below.Now, on the interval [t 0 , t 1 ) we calculate the derivative of the function ||g(t)|| with respect to system (24).We have the equality: Using Lemma 1 and the inequalities from (4), we evaluate the upper boundary of the expressions inside the square brackets in formula (26).At first, we consider the third expression, and the corresponding chain of inequalities: where C = 2u 2 max + v 2 max .Now, we consider the second expression, and the following chain of inequalities: where B = 5u 2 max + 9v 2 max + 4γ 2 + 4δ 2 .Finally, we consider the first expression, and the corresponding inequality: For the first two terms of the right hand side of this relationship, we obtain the inequalities: Next, let us find the eigenvalues of matrix Q 3 .We have the formulas: Then, for the last term in the right hand side of relationship (27) the inequality is valid: Using inequalities ( 28), (29), we continue evaluating the right hand side in (27): where A = 7γ 2 4 + δ 2 .Substituting the obtained relationships for the expressions in the square brackets into formula (26), we finally have the differential inequality: Now, let us consider the quadratic equation We define the sign of its discriminant: It is easy to see that this discriminant is positive.Let us denote by K 0 the biggest root of equation ( 31).Then, the following formula holds: Now, for all vectors g ∈ 3 such that ||g|| ≥ ρ we introduce the function W (g) = ||g|| + K 0 , and then rewrite differential inequality (30) for this function.We have the inequality: After its transformations regarding relationship AK 2 0 − BK 0 + C = 0 we find the differential inequality: Finally, let us consider the auxiliary Cauchy problem: It is easy to see that the value h 0 satisfies the inequality Let us find the solution to Cauchy problem (33).We solve the corresponding Bernoulli equation and satisfy the initial condition.Then, the following formula can be obtained: In this equality we assume that the values ρ and t 0 are such that the expression inside the brackets is defined for all t ∈ [t 0 , t 1 ].We can do this, for example, by choosing for a given ρ, a value t 0 such that the difference (t 1 − t 0 ) is sufficiently small.By inequality (34) the expression inside the square brackets in (35) is negative, and hence h(t) is a finite positive function increasing on the interval [t 0 , t 1 ].Then, we have the inequality h(t) < h(t 1 ) for all t ∈ [t 0 , t 1 ).Therefore, from differential inequality (32), Cauchy problem (33), the Chaplygin's Theorem ( [17]), and the equality: the chain of inequalities follows: which contradicts relationship (25).Our assumption was wrong.Therefore, there exists the solution g(t) for system (24), and hence for system (23), which is defined on the entire interval [0, T ].Then, system ( 21) is also defined on this interval.We rewrite this system using the components of the solution Finally, at the third stage, we show that the function r(t) = L(t) has at most two distinct zeros on the interval [0, T ].Again, we assume the contradiction.Let the function r(t) has on the interval [0, T ] at least three distinct zeros τ 1 , τ 2 , τ 3 , where 0 ≤ τ 1 < τ 2 < τ 3 ≤ T .Applying the generalized Rolle's Theorem ( [18]) to the first differential equation of system (36), we conclude that the function q(t) has on the interval (0, T ) at least two distinct zeros η 1 , η 2 , where 0 < η 1 < η 2 < T .Applying now this theorem to the second differential equation from (36), we find that the function µ(t) has at least one zero ξ on the interval (0, T ).However, this function satisfies the linear homogeneous differential equation, and therefore µ(t) = 0 everywhere on the interval [0, T ].Considering then the second equation of system (36), we conclude that the function q(t) is also zero everywhere on the interval [0, T ].Therefore, using the first equation of this system, the similar conclusion is obtained for the function r(t), i.e., r(t) = 0 for all t ∈ [0, T ].This fact contradicts Lemma 2. Hence, our assumption was wrong, and the switching function L(t) has at most two distinct zeros on the interval [0, T ].
Let us note one important fact of the proof of Lemma 3. Due to various ways of evaluation of expressions inside square brackets in formula (26), the values A , B, C are not uniquely defined.They can be made as big in value or not.Analysis of the formula of K 0 , relationship (35), and the initial condition at Cauchy problem (33) shows that the values B, (B 2 − 4AC) should be the smallest possible.Therefore, the values A , C should be the biggest possible as well, and satisfy the inequality 4AC < B 2 .
Arguments similar to those presented in Lemma 3, were previously used in [19,20] to estimate the number of zeros of switching functions in optimal control problems for the SIR model of the epidemic spread.

Solution to Optimal Control Problem
We describe the method of solving the optimal control problem (10).At first, we define the set: Next, for arbitrary point θ ∈ Ω we define the controls u θ (t), v θ (t) by formula: It is easy to see that the controls u θ (t), v θ (t), defined in this way, include the type of the optimal controls u * (t), v * (t) given by formula (37).Finally, substituting the controls u θ (t), v θ (t) into system (6), and integrating this system over interval [0, T ], we find the functions s θ (t), e θ (t), i θ (t), which correspond to these controls.Then, using the controls u θ (t), v θ (t) and the functions e θ (t), i θ (t), we find the value of functional J(u θ , v θ ) by formula (10), corresponding to these controls, or, in other words, to the point θ ∈ Ω.
Thus, the function of two variables is defined: and problem (10) is now reduced to the problem of constrained minimization: Methods for numerical solution of such problem are well-known (see [21]).Problem (40) is considerably simpler than optimal control problem (10) and can be solved numerically by using, for example, the overlapping method.Using this method assumes a knowledge of the Lipschitz constant of the function F (θ).The following lemma is devoted to finding of this constant.4-6 show the graphs of the corresponding optimal controls u * (t), v * (t).The control u * (t) is shown in red, and the control v * (t) is shown in blue.For these Ebola outbreaks, Table 2 shows the dependence of the values of the switchings θ * 1 , θ * 2 of the controls u * (t), v * (t), and the optimal value of the functional J * = N J(u * , v * ) on the values of the weighting coefficients σ and ν.
Let us discuss the results presented by Table 2 and obtained for Zaire (1976).We can see that the type of the optimal controls u * (t), v * (t) for 0.01 ≤ σ = ν ≤ 0.1 is very different from the type of these controls corresponding to 0.5 ≤ σ = ν ≤ 10.In the first case, if 0.01 ≤ σ = ν ≤ 0.1, the moments of switching θ * 1 , θ * 2 of the controls u * (t), v * (t) are almost the same, which means that almost no intervention is needed in order to minimize the functional from (10), and all preventive control measures would take less than one day.In the second case, if 0.5 ≤ σ = ν ≤ 10, the type of the optimal controls u * (t), v * (t) are qualitatively different, that is θ * 1 = 0, and the switching θ * 2 moves to the end of the interval [0, T ] as the values of the weighting coefficients σ, ν increase.In particular, if σ = ν = 0.5, active intervention control measures take place during the first 29-30 days, then they stop.If σ = ν = 1, then these measures continue at the maximum effort during the first 80-81 days.Finally, if σ = ν = 10, then intervention control measures continue at the maximum rate almost during the entire time interval (90 days).Note that similar behavior of the optimal controls u * (t), v * (t) can be observed for the Ebola outbreaks in Uganda (2000) and Liberia (2014).Such behavior of these controls is explained by the fact that the weighting coefficients σ and ν from formula (8) determine the relative importance of the functional J 1 (u, v) to J 2 (u, v).If the values of these weighting coefficients are small (0.01 ≤ σ = ν ≤ 0.1), then the costs of intervention control constraints are relatively more important than the total fractions of exposed and infected individuals, and vice versa.Finally, Figures 7-9 represent the surfaces describing the dependence of values of the function J(θ) = N J(u θ , v θ ) on the moments of switching θ = (θ 1 , θ 2 ) ∈ Ω.One of the points of such surfaces correspond to the optimal values J * .Analyzing these surfaces it is easy to see that the function J(θ) has not only the global minimum, the search of which is the subject of our investigation, but also local minima.This once again demonstrates that the Pontryagin maximum principle is only a necessary optimality condition.

Conclusions
In this paper, a problem of minimizing the weighted sum of total fractions of infected and exposed individuals and total costs of intervention control constraints over a given time interval for the SEIR control model describing the spread of an Ebola epidemic in a population of a constant size was considered.Analysis of the optimal solutions of this problem was made using the Pontryagin maximum principle.The established properties of the switching function allowed the reduction of the original optimal control problem to the solution of the problem of a constrained minimization of the function of two variables.Results of the numerical calculations and their analysis are presented.