ON OPTIMAL SINGULAR CONTROLS FOR A GENERAL SIR-MODEL WITH VACCINATION AND TREATMENT

. A general SIR-model with vaccination and treatment is considered as a multi-input optimal control problem over a ﬁxed time horizon. Existence and local optimality of singular controls is investigated. It is shown that the optimal vaccination schedule can be singular, but that treatment schedules are not.


1.
Introduction. In spite of the improvement in sanitation, developments of antibiotics and vaccines, infectious diseases still contribute significantly to deaths worldwide. While the earlier recognized diseases like cholera or the plague still sometimes create problems in underdeveloped countries erupting occasionally in epidemics, in the developed countries new diseases are emerging like AIDS (1981) or hepatitis C and E (1989)(1990). Additionally, some of the diseases that we generally seem to have under control like tuberculosis, pneumonia, or gonorrhea develop antibiotic-resistant strains. Malaria or yellow fever continue to be a major problem in regions with climate changes. New threads constantly appear like the recent bird flu (SARS) epidemic in Asia or the very dangerous Ebola virus in Africa. Overall, infectious diseases continue to be one of the most important health problems worldwide.
Modeling of epidemiological phenomena has a very long history with the first model for smallpox formulated by Daniel Bernoulli in 1760. From the early twentieth century, in response to epidemics of various infectious diseases, a large number of models has been constructed and analyzed starting a field known as epidemiology modeling; for example, see [7,8]. Although many of these models are disease specific, like models for measles, smallpox, malaria, gonorrhea or HIV/AIDS (e.g., [9,19]), there are also common characteristics of all models which allows to analyze a general framework without a specific disease in mind [8].
Mathematical models have become an important tool in describing the dynamics of the spread of an infectious disease and the effects that vaccination and treatment have on its dynamics. These models can be particularly useful in comparing the effects of various prevention, therapy and control programs. Since a variety of these 982 URSZULA LEDZEWICZ AND HEINZ SCHÄTTLER programs are available, it is a natural objective to design optimal programs in terms of some pre-assumed criterion. This brings the application of the tools of optimal control to these problems. Optimal control has a long history of being applied to problems in biomedicine, particularly, to models for cancer chemotherapy. This research, that started in the seventies and eighties, e.g., [5,16,17], has continued uninterrupted to the present day [11,14,18] and has taken on new forms as mathematical models for novel therapies such as tumor anti-angiogenesis are formulated and analyzed, for example, see [12,13]. But until recently, little attention has been given to models in epidemiology. Gaff and Schaefer [6] consider several variations of standard epidemiologic models with an objective that takes into account vaccination and treatment. They compute optimal controls numerically and investigate the sensitivity of optimal solutions to parameter values. Behncke also analyzes SIR epidemiological models numerically and finds controls that exert maximum effort from the beginning [1]. In this paper, for a SIR-model with vaccination and treatment, we investigate the optimality of singular controls theoretically. These are controls that correspond to time-varying administration schedules and while we show that they cannot be optimal for treatment protocols, depending on the growth model taken, they become an option for optimal vaccination schedules.
2. Epidemiologic Models. The standard models in epidemiology include five compartments typically denoted by M , S, E, I and R corresponding to five classes of individuals [6]. The most typical epidemiology model, however, is of type SIR and only includes three of these classes. The class S represents the susceptible individuals, I stands for the class of infected ones and R denotes the recovered class of those who went through infection and emerge with permanent or temporary infection-acquired immunity. This model can be extended by adding a class M that represents infants with passive immunity. But generally maternal antibodies are cleared from the body relatively quickly and the infant enters the class S. Another class E (located between the classes S and I) that is taken into account in some models represents exposed individuals who are in the latent stage, i.e., have been infected, but are not yet infectious. In this paper, in order to keep the dimension of the model smaller, as a first step we only consider an SIR type model. Furthermore, depending on the type of disease considered, the immunity in the class R may not be permanent and the class R should be followed by the class S of individuals who regain their susceptibility when temporary immunity ends. Waning immunity clearly is important, but in the simplified model considered here we assume that immunity is permanent.
Let S(t) represent the number of susceptible individuals, I(t) the number of infective ones and R(t) the number of recovered ones, all at time t. We also denote the total number of individuals by N , N = S + I + R, and, like in [6] assume that all new births enter the susceptible class S. Naturally, it will depend on the disease specifics to what extent this holds. Different from classical formulations of SIR-models in the literature, however, we will not formulate the dynamics in terms of the variables S, I, and R, but we prefer to replace R by the total number of individuals, N . We believe this somewhat simplifies the structure of the equations and better brings out the theoretical connections between the variables. Therefore we consider the following dynamics similar to [6]: Equation (1) models the overall growth of the population. We only assume that (F): F : [0, K) → R + is a continuously differentiable, strictly increasing function defined on [0, K) with F (0) = 0 and K ≤ ∞ denoting a fixed carrying capacity for the population.
All standard growth models have these properties; e.g., with µ being a positive growth rate and K the carrying capacity, exponential growth given by F (N ) = µN The function F describes the normal growth of the population without any epidemic present and the term δI represents additional decrease through higher deaths of infectious individuals.
It is assumed that all new births will be susceptible individuals and thus the growth F (N ) enters equation (2). The main term in this dynamics, β IS N , represents the loss of the number of susceptible individuals that are being infected by individuals from class I with the parameter β standing for the average number of adequate contacts (i.e., contacts sufficient for transmission) of a person per unit of time. The last term, Su, represents the effect of vaccination. The variable u is a control that represents the rate at which susceptible individuals are vaccinated. It takes values in a compact interval, 0 ≤ u ≤ u max , and it is assumed that vaccination removes the fraction Su of individuals from the class S and makes them resistant.
In the I-dynamics, equation (3), β IS N thus represents the inflow of newly infected individuals from class S. The term γI models the outflow of infected individuals that recover without intervention while, in accordance with (1), δI represents deaths from the infection; accordingly γ and δ denote the recovery and the death rates, respectively. The additional outflow Iv is related to the cure of infected individuals due to treatment and v represents the rate at which infectious individuals are treated at each time period, the second control in the model with values in the interval Thus there are two possible mechanisms available as controls: immunization of the susceptible individuals and treatment of the infected ones. These actions are modeled by the two controls u and v that for mathematical reasons are taken as Lebesgue-measurable functions. The action of both controls enriches the class R of the recovered individuals by removing them from the class of susceptible and infected ones, respectively. The class R is defined as R = N − I − S. For the model to be realistic, we need to make sure that all the variables including R remain positive. The initial numbers of individuals in each of the populations are positive numbers denoted by Note that if there are no infected individuals initially, I 0 = 0, I remains identically zero. The model thus does not represent the onset of infection, but only its course.
But this is not possible in finite time τ < ∞: Let P = P (t) be the solution to the initial value probleṁ and it follows from standard comparison results (see, e.g., [10],) that the solution P (t) will always lie below N (t). Since F ∈ C 1 and F (0) = 0, we can write F in the form F (N ) = N φ(N ) with some continuous function φ. In a given small neighborhood W of P = 0 the right-hand side of the differential equation is therefore continuously differentiable and linearly bounded. Hence solutions are unique, P (t) ≡ 0 is an equilibrium solution and the solution P (t) can only approachP (t) ≡ 0 as t → ∞. In particular, it is bounded away from 0 on any finite interval. Hence N (t) > P (t) will be positive for all finite times.
3. Formulation as an Optimal Control Problem. Our goal is to solve the following problem: given initial population sizes of all three classes, S 0 , I 0 and R 0 , find the best strategy in terms of combined efforts of vaccination and treatment that would minimize the number of people who die from the infection while at the same time also minimizing the cost of vaccination and treatment of the population. Naturally, there are various ways of expressing such a goal mathematically. In this paper, for a fixed terminal time T , we consider the following objective: The first term in the objective, aI(t), represents the number of people who are infected at time t and is taken as a measure for the deaths associated with the outbreak. The terms, bu(t) and cv(t) represent the cost of vaccination and treatment, respectively, and are assumed to be proportional to the vaccination and treatment rates. This is a simplified first approach and does not reflect, for example, that the cost of efforts needed of vaccinating the population will increase drastically the closer we get to the saturation value K. Real life data show that the cost of vaccinating the first 80% of the population is relatively small compared with the cost for the remaining 20% (see [6]). But our model will be realistic if such high rates of vaccination are not reached. We shall apply methods of geometric optimal control theory to analyze the relations between optimal vaccination and treatment schedules. These techniques become more transparent if the problem is formulated as a Mayer-type optimal control problem; that is, one where we only minimize a penalty term at the terminal point. Such a structure can easily be achieved at the cost of one more dimension if the objective is adjoined as an extra variable. Defininġ we therefore consider the following optimal control problem: Introducing the state x = (Z, N, S, I) T , the dynamics of the system is a multiinput control-affine system of the forṁ with drift vector field f given by and control vector fields g 1 and g 2 given by We call an admissible control pair (u, v) with corresponding solution x a controlled trajectory of the system.

Necessary Conditions for Optimality.
First-order necessary conditions for optimality of a controlled trajectory are given by the Pontryagin maximum principle [4,15]: For a row-vector λ = (λ 1 , λ 2 , λ 3 , λ 4 ) ∈ (R 4 ) * , we define the Hamiltonian H = H(λ, x, u, v) as the dot product, , ·, · , of the row vector λ with the column vector that defines the dynamics, that is Then, if (u * , v * ) is an optimal control defined over the interval [0, T ] with corresponding trajectory x * = (Z * , N * , S * , I * ) T , there exists an absolutely continuous co-vector, λ : [0, T ] → (R 4 ) * , such that the following conditions hold 1 : (a) λ satisfies the adjoint equations (written as a row vector and with Df and Dg i denoting the Jacobian matrices of the partial derivatives) with terminal condition (b) for almost every time t ∈ [0, T ] the optimal controls (u * (t), v * (t)) minimize the Hamiltonian along (λ(t), x * (t)) over the control set [0, u max ] × [0, v max ] and (c) the Hamiltonian is constant along the optimal solution. We call a pair (x, (u, v)) consisting of admissible controls (u, v) with corresponding trajectory x for which there exist multipliers λ such that the conditions of the Maximum Principle are satisfied an extremal (pair) and the triple (x, (u, v), λ) is an extremal lift (to the cotangent bundle).
Note that the dynamics does not depend on the auxiliary variable Z and thus by the adjoint equation (6) the multiplier λ 1 is constant; by the terminal condition (7), it is thus given by λ 1 (t) ≡ 1. In particular, the overall multiplier λ(t) is never zero. For almost any time t, the optimal controls (u * (t), v * (t)) minimize the Hamiltonian Since H is linear in the controls, this minimization problem splits into two separate one-dimensional problems that can easily be solved. Defining the so-called switching functions Φ 1 and Φ 2 as Φ 1 (t) = λ(t), g 1 (x * (t)) = b − λ 3 (t)S * (t) and Φ 2 (t) = λ(t), g 2 (x * (t)) = c − λ 4 (t)I * (t), it follows that the optimal controls satisfy The minimum condition alone does not determine the control at times when Φ i (t) = 0. If Φ i (τ ) = 0, butΦ i (τ ) = 0, then the control switches between the value 0 and its maximum value depending on the sign ofΦ i (τ ). Controls with this property are called bang-bang controls and we refer to the constant controls with values in the endpoints of the control intervals as bang controls. The other extreme occurs when a switching function vanishes over an open interval. In this case also all derivatives of Φ i (t) must vanish and this typically allows to compute such a control. Controls of this kind are called singular [2]. While the name (which has historical reasons) might give the impression that these controls are less important, quite the contrary is true. Singular controls (if they exist) tend to be either the best (minimizing) or the worst (maximizing) strategies and in either case they are essential in determining the structure of optimal controls. These typically then need to be synthesized from bang and singular controls through an analysis of the switching function. Thus singular controls generally play a major role in a synthesis of optimal controlled trajectories and in this paper we analyze their existence and local optimality for the problem [OC]. An essential tool in this analysis is the Lie bracket of vector fields which naturally arises in the formulas for the derivatives of the switching function. Given two differentiable vector fields f and g defined on a common open subset of R n , their Lie bracket can be defined as The Lie-bracket is anti-commutative, i.e., [f, g] = −[g, f ], and for arbitrary vector fields f , g and h it satisfies the Jacobi identity [3] [f, [g, h]] The following result provides an elegant and important framework for efficient computations of the derivatives of the switching functions. It is easily verified by a direct computation.

Proposition 2.
Let (x, (u, v)) be a controlled trajectory of the system and let λ be a solution to the corresponding adjoint equations. Given a continuously differentiable vector field h, define Ψ(t) = λ(t), h(x(t)) . Then the derivative of Ψ is given bẏ

5.
The Structure of Singular Controls. We investigate the existence and local optimality of singular controls for the system [OC]. By Proposition 2 the derivatives of the switching functions Φ 1 (t) = λ(t), g 1 (x(t)) and Φ 2 (t) = λ(t), g 2 (x(t)) are given byΦ By anti-commutativity of the Lie bracket [g i , g i ] ≡ 0 and a simple computation verifies that the control vector fields g 1 and g 2 commute, i.e., [g 1 , g 2 ] ≡ 0 as well. We thus have thaṫ Elementary calculations verify that We first analyze the control u, i.e., vaccination schedules. Applying Proposition 2 once more toΦ 1 , it follows thaẗ A direct calculation shows that g 2 and [f, g 1 ] commute as well, [g 2 , [f, g 1 ]] ≡ 0, and that

URSZULA LEDZEWICZ AND HEINZ SCHÄTTLER
The relationΦ implies that λ(t), [g 1 , [f, g 1 ]](x(t)) = −2λ 3 (t)F (N (t)) and Φ 1 (t) = b − λ 3 (t)S(t) ≡ 0 gives that λ 3 (t) must be positive along a singular arc. Hence we have that Singular controls of this type, i.e., for which λ(t), [g 1 , [f, g 1 ]](x(t)) does not vanish, are said to be of order 1 and it is a second-order necessary condition for minimality, the so-called Legendre-Clebsch condition, that this quantity be negative [2]. Thus for this model singular controls u are locally optimal. Furthermore, in this case, and taking into account that [g 2 , [f, g 1 ]] ≡ 0, we can compute the singular control as Evaluating the vector fields, this equation can be simplified. A direct, but somewhat lengthy computation shows that for some smooth functions ϕ and ψ we can write Since λ(t), [f, g 1 ](x(t) ≡ 0, it follows from (11) that u sin (t) = 1 2 δβI(t)S(t) N (t))F (N (t)) − ψ(x(t)).
Once more using (10), the first term can be simplified to 1 2 δ λ2 λ4 and we obtain the following result: Proposition 3. A singular control u is of order 1 and satisfies the Legendre-Clebsch condition for minimality. The singular control is given as a function depending both on the state and the multiplier in the form with ψ given by (12).
Proposition 4. The control v cannot be singular.
Proof. Suppose v is singular on an open interval J. The second derivative of the switching function Φ 2 is given bÿ ](x(t)) .
A singular control v is therefore of order higher than 1. Hencë and additional non-trivial conditions need to be satisfied that are obtained by differentiating (14) further. These conditions prevent the existence of singular arcs. Since λ 1 ≡ 1, the condition that λ vanishes against the vector fields g 2 , [f, g 2 ] and [f, [f, g 2 ]] determines the multiplier λ(t). Differentiating (14), we obtain Φ (3) Using the Jacobi identity once more, it follows that Hence the multiplier λ in addition needs to vanish against the vector field [f + g 1 u, [f, [f, g 2 ]]](x(t)). this leads to additional restrictions on either λ or on the second control u that would need to be singular. In either case, these conditions can not be satisfied and lead to contradictions.
6. Conclusion. In this paper we initiated the analysis of an optimal control problem for an SIR-model with vaccination and treatment. We analyzed the structure of singular controls, but it still remains to complete this analysis by determining the structure of possible concatenations with bang-bang controls in order to determine an optimal synthesis of controlled trajectories. Based on our computations so far, it is expected that the optimal treatment schedule will be bang-bang, most likely with just one switch from u max to u = 0 while we expect a singular regimen for the optimal vaccination strategy.