OPTIMAL CONTROL APPLIED TO VACCINATION AND TREATMENT STRATEGIES FOR VARIOUS EPIDEMIOLOGICAL MODELS

. Mathematical models provide a powerful tool for investigating the dynamics and control of infectious diseases, but quantifying the underlying epidemic structure can be challenging especially for new and under-studied diseases. Variations of standard SIR, SIRS, and SEIR epidemiological models are considered to determine the sensitivity of these models to various parameter values that may not be fully known when the models are used to investigate emerging diseases. Optimal control theory is applied to suggest the most ef-fective mitigation strategy to minimize the number of individuals who become infected in the course of an infection while eﬃciently balancing vaccination and treatment applied to the models with various cost scenarios. The optimal control simulations suggest that regardless of the particular epidemiological structure and of the comparative cost of mitigation strategies, vaccination, if available, would be a crucial piece of any intervention plan.


1.
Introduction. Accidental and intentional introduction of infectious diseases to previously naïve geographic regions has brought more focus and attention to the development of response plans to such scenarios. All levels of government and public health officials are searching for answers to identify the best strategies for intervention prior to what may be an inevitable event. Given our global connectedness, diseases that were previously isolated to various parts of the world have spread even across oceans [4].
In developing response plans to disease outbreaks, decision makers are often faced with trade-offs in choosing among various treatment options, quarantine, and/or vaccines. For example, cholera, human papilloma virus, and influenza all have approved medical treatment options as well as available vaccines. The challenge is to find the optimal response balancing treatment and vaccination that will minimize incidence and disease-related mortality while being mindful of the costs of each mitigation strategy.
Mathematical models have been used to explore the dynamics of diseases since Bernoulli used a simple model to estimate smallpox mortality and argue for the advantages of inoculation [3]. Most mathematical models for infectious diseases use a similar underlying methodology based on the work formalized by Anderson and May [1]. The first step is to represent the epidemiology of the disease being studied by dividing the population into subpopulations, called compartments, that represent the various stages of disease progression. For example, individuals are identified as "susceptible" (S) to a disease if they don't currently have the disease nor any immunity to the disease, e.g., they have not been vaccinated. Individuals are "exposed" (E) if they have been infected with the disease pathogen but are not able to infect others, and they are "infectious" (I) if they are infected and infectious with the disease pathogen. Finally, they are "removed" (R) if they have cleared the infection and have immunity to recurrence for at least some period of time, have been successfully vaccinated against the disease, or are otherwise isolated from the population so that they cannot spread the disease. Using these compartments and possibly additional ones as needed, the epidemiology of a disease is represented as a series of subpopulations connected by the flow from one compartment to another that is dictated by the disease dynamics. The rates of flow between compartments are estimated from experiments and field data.
Hethcote provides several notable surveys of basic epidemiological models [13,16], and he and other authors have considered in detail the effects of variations on these models. In addition, there have been several articles considering optimal control applied to specific diseases, for example [18,19,22,31]. Hethcote suggests that the most needed contributions to the literature would contrast the effects of altering the basic components of the various models [13].
Brauer considers several less traditional variations, including compartments for asymptomatic, quarantined, and isolated individuals [5]. In considering the effect of treatment in these models, he finds that asymptotically the models' qualitative behaviors agree with those of more simple models. However, he notes that it would be of interest to consider the trade-off between vaccination and treatment and how the more complicated structures might influence those results.
Behncke applies optimal control to a number of general epidemiological models with general terms allowing many variations of those models [2]. The article deals mostly with the costs and timely application of vaccination to various epidemiological structures, and provides clear theoretical results in those cases.
Several key components of such models are those that consider mitigation strategies such as treatment and vaccination. Optimal control theory provides a valuable tool to begin to assess the trade-offs between vaccination and treatment strategies. Optimal control is a mathematical technique derived from the calculus of variations. There are a number of different methods for calculating the optimal control for a specific mathematical model. For example, Pontryagin's maximum principle allows the calculation of the optimal control for an ordinary differential equation model system with a given constraint. Variations of Pontryagin's maximum principle have been derived for other types of models including partial differential equations and difference equations [20,23]. These techniques are powerful when applied to disease models and can provide great insights into the best pathway to reduce disease burden. For example, with a given mathematical model for a disease, one can calculate the best vaccination schedule balancing the cost of the vaccine and the cost of the disease burden [15].
In using optimal control to suggest an intervention strategy for a specific disease, one of the first steps is to identify the epidemic structure of a disease. During an emerging infectious disease, it is often unclear what the exact epidemic structure should be, as can be seen from the variety of models created for the 2002-2003 SARS epidemic [24,29,33]. Even for existing diseases, there is often a lack of agreement on whether a given disease has, for example, a significant incubation period prior to the individual being infectious as is the case for smallpox [8,21]. These types of questions are vital to the understanding of the disease dynamics in general. The process of developing a mathematical model forces quantitative thinking that often highlights gaps in data and current knowledge that need to be addressed.
The question we study in this paper is whether this underlying epidemic structure significantly impacts the predicted optimal control strategy for administering vaccination and/or treatment. This is an important question given that often at the beginning of a newly emerged or previously under-studied disease, the underlying epidemic structure is unknown. If the optimal control strategy remains qualitatively the same regardless of the underlying epidemic structure, then less time needs to be spent investigating that structure before initial response policies can be made. Additional research will be needed subsequently to understand the disease more fully and assess and update the response policy as needed. We present results from the exploration of SIR, SIRS, and SEIR models including vaccination and treatment options under varying rates for incidence and disease related death. Below we present the SIR model analysis in detail including the proof of uniqueness and existence of the optimal control solutions followed by summaries for the results of the similar analyses of the SIRS and SEIR models. Finally, we present the results of numerical simulations for each model under various parameter values.

2.
Equations for a standard SIR model. We begin our explorations with a standard SIR model for a population with underlying logistic growth as defined by the following state equations: subject to the boundary conditions As described in the previous section, the variables S, I, and R represent the susceptible, infectious, and removed classes, respectively. We assume logistic growth of the total population N = S + I + R with carrying capacity K, and we additionally assume that all new births enter the susceptible class. Note that many models treat the population size as constant, and in the short term the models provide essentially the same results. However, when considering disease-related deaths [13] the logistic model should provide more accurate predictions. It should be noted that the direct applicability of a logistic growth model to a human population remains a subject of debate [27], as human birth and death rates are far more complicated than simple resource allocation. Nonetheless, the logistic model, which provides mathematical stability, will produce reasonable results for our generic model population over a time frame that is significantly shorter than a human lifespan as is used in this optimal control study.

HOLLY GAFF AND ELSA SCHAEFER
The parameter β approximates the average number of contacts with infectious individuals needed to make one person ill in each unit of time. Note that I N approximates the fraction of the population that is infectious. If we assume that each person has one contact in each time period (measured in days for our analysis), and thus the term SI N approximates the number of susceptible/infectious contacts in each time period, and β SI N , the "standard incidence," approximates the number of susceptibles who become infectious in each time period [13]. Liu, Hethcote, and others have published multiple studies examining the need for and effect of other measures of incidence [12,14,17,25,26]. In particular, if I is raised to a power, interesting dynamics can be observed. As these researchers write, it is difficult to assess in a particular case which incidence model is most appropriate.
We assume that some proportion δ of the infectious population dies due to the disease in each time period. The control function ν(t) measures the rate at which susceptibles are vaccinated in each time period, and the control function τ (t) measures the rate at which infectious individuals are treated in each time period. Note that in equations (1) -(3) the control ν moves individuals from the S class to the R class, and the control τ moves individuals from the I class to the R class. Since quarantine effectively removes infectious individuals from the population, the control τ can provide a rudimentary model for quarantine effects in addition to modeling medical treatments of the disease. This observation may be useful for determining mitigation strategies for diseases for which no treatment is known. Note that an alternative model for quarantine introduced by Hethcote et al. adjusts the incidence term to account for a smaller mixing population, which can produce periodic disease dynamics for some parameter sets [12]. In addition, as mentioned earlier, Brauer's recent model also distinguishes between quarantine and isolation, and also adjusts other key parameters for these populations [5].
A summary of parameter meanings and values used in our numerical computations appears in Table 1. Note that for these simulations the total initial population is at carrying capacity, and the models' logistic population dynamics allow the the population to recover its size as disease-related deaths are permitted. Also, note that initial conditions have been chosen to allow our optimal control pairs to be well-defined in equation (14) that will be derived in Section 4 following.
Given initial population sizes S 0 , I 0 , and R 0 , we seek the best mitigation strategy for the outbreak modeled in equations (1) -(3) by optimally defining bounded, Lebesque integrable control functions ν(t) and τ (t). Our goal is to minimize the number of people who become infected, and thus also the number of people who die due to the infection, while also minimizing the effort of vaccinating and treating the population. Thus, we seek to minimize the objective functional where m ≥ 1. The constants B 1 , B 2 , and B 3 have a dual role. On one hand, they are needed to balance the units in the integrand because the number of infected individuals will be measured in the hundreds in this paper, while τ and ν are treatment rates and will necessarily lie between 0 and 1. On the other hand, in numerical runs we further vary the constants to place stronger importance on the minimization of infected individuals or on treatment and/or vaccination efforts. The terms B 1 I(t) and B 3 τ 2 (t) are rather standard terms in objective functionals with similar goals (see, for example, [18,19], or, for a contrasting objective functional  [22,30]). Treatment is viewed as a nonlinear function since implementation of any public health intervention does not have a linear cost, but rather, there are increasing costs with reaching higher fractions of the population. The quadratic is the most simple nonlinear function and so is used for treatment. The remaining term in our objective functional seeks to increase the expense of vaccination when most of the population has either been vaccinated or has immunity from a prior infection. When m is chosen to be an even integer such as 10, the term B 2 [ R K ] m ν 2 (t) reflects the reality that the cost of vaccinating the first approximately 80 percent of a population is relatively small as compared with the cost of vaccinating the remaining 20 percent. For example, this could reflect the costs of lab tests that could used to determine a person's immune status prior to vaccination. Indeed, numerical investigations suggested unrealistic vaccination schemes when the vaccination term was allowed to have the more standard quadratic form m = 2. An alternate way to address the cost differential caused by the density of succeptibles in the population is suggested by Behncke, who introduces a cost function and otherwise maintains a linear functional resulting in a bang-bang control for the vaccination strategy [2].
We assume there are practical limitations on the maximum rate at which individuals who may be vaccinated or treated in a given time period and we define positive constants ν max and τ max accordingly. We define the set of admissible controls to be We seek an optimal control pair (ν * , τ * ) such that As is listed in Table 1, we assume a maximum vaccination rate of 0.1 and a maximum treatment rate of 0.6. The maximum vaccination rate is based on the 1947 smallpox outbreak in New York City. At the height of this outbreak, a half million to a million of the population of 6 million were vaccinated daily [32]. The maximum treatment rate was chosen more arbitrarily, but qualitative results were not sensitive to a change of this parameter.
3. Existence of an optimal control pair. We begin by examining sufficient conditions for the existence of a solution to the optimal control problem. A high death rate from disease could theoretically cause the population to vanish, and the existence theorem must assume that the external death rate δ is smaller than the population intrinsic growth, or turnover, rate µ.
Proof. We refer to the conditions in Theorem III.4.1 and its corresponding Corollary in Fleming and Rishel [10]. The requirements there on the set of admissible controls Ω and on the set of end conditions are clearly met here. The following nontrivial requirements from Fleming and Rishel's theorem are listed and later verified below: A. The set of all solutions to system (1) -(4) with corresponding control functions in Ω (as given in equation (6)) is nonempty. B. The state system can be written as a linear function of the control variables with coefficients dependent on time and the state variables. C. The integrand L in equation (5) is convex on Ω and additionally satisfies In order to establish condition A, we refer to Theorem 3.1 by Picard-Lindelöf in Coddington and Levinson [7]. If the solutions to the state equations are a priori bounded and if the state equations are continuous and Lipschitz in the state variables, then there is a unique solution corresponding to every admissible control pair in Ω.
Thus, we begin establishing bounds on N = S + I + R, and, by extension, S, I, and R. Note that N satisfies the modified logistic equation and assume that N (0) = N 0 . If N > K, then N is decreasing. Thus, N is bounded above by max(N 0 , K). Note that the quantities S, I, and R decrease only proportional to their present sizes, respectively, and thus none of S, I, and R can be negative. Therefore, the upper bound for N is also an upper bound for S, I, and R. Note that equations (1) and (2) have division by N , and thus we need a strictly positive lower bound for N . To produce this bound, we note that if the solution to the differential equation is bounded below then surely the solution to (8) is bounded below by the same bound. Since dÑ dt is proportional toÑ we know that in finite timeÑ remains strictly positive. Thus, for a given end-time T solutions to (8) will have a positive lower bound.
With the bounds established above, it follows that the state system is continuous and bounded. It is equally direct to show the boundedness of the partial derivatives with respect to the state variables in the state system, which establishes that the system is Lipschitz with respect to the state variables (see [6], page 248). This completes the proof that condition A holds.
Condition B is verified by observing the linear dependence of the state equations on controls ν and τ . Finally, to verify condition C we note that since the integrand L of the objective functional is quadratic in the controls, L is clearly convex in the controls. To prove the bound on the L we note that by the definition of Ω we have B 2 ν 2 ≤ B 2 , and so B 2 ν 2 − B 2 ≤ 0. Thus, 4. Characterization of optimal controls. We apply Pontryagin's Maximum Principle [28] and convert the optimization problem described in the previous section to the problem of finding the point-wise minimum relative to ν and τ of the Hamiltonian.
Theorem 4.1. Given the optimal control pair (ν * , τ * ) and corresponding solutions to the state system (1) -(4) S * , I * , R * that minimize the objective functional (5) there exist adjoint variables λ 1 , λ 2 , and λ 3 satisfying with transversality conditions Furthermore, as long as the optimal removed class R * is nonzero, we may characterize the optimal pair by the continuous functions We note that the restriction on R * can easily be resolved by assuming a nonzero initial value for R 0 for mathematical convenience. Physically this requirement is not significant.
Proof. The result follows from a direct application of a version of Pontryagin's Maximum Principle for bounded controls (see [20,23,28]). We form the Hamiltonian H: As dictated by the Maximum Principle, the adjoint equations are given by the equations dλ1 dt = − ∂H ∂S , dλ2 dt = − ∂H ∂I , dλ3 dt = − ∂H ∂R , and must satisfy transversality conditions λ i (T ) = 0 for values i = 1, 2, 3. Finally, the optimality conditions dictate that ∂H ∂ν = ∂H ∂τ = 0 for the optimal pair (ν * , τ * ) on the interior of the control set, and this condition is simplified in equations (14) with attention to the bounds on the pair as given in the definition of Ω in equation (6). Note that as a result of the transversality condition, the optimal vaccination and treatment will be zero at the end time. Observe that by the characterization of the controls given in (14) and the nonzero assumption for R, it follows that the controls are continuous in time.
The optimality system is defined as is the compilation of the state equations (1) -(3), the initial conditions (4), the adjoint equations (10) - (12), and the transversality conditions (13), with the optimality equations (14) substituted into the state and adjoint equations. Uniqueness of the optimality system can be shown for a small time interval.
Theorem 4.2. For T sufficiently small the optimality system is unique.
The proof of the the theorem is rather lengthy due to the complexity of the optimality system. The approach is explained in [9,18], and the proof itself is given in the appendix.

5.
Variations of the SIR model. As Hethcote details in [13], there are hundreds of variations of the SIR model described in Section 2. Two classic variations result from adding an incubation period before a subject becomes infectious, the SEIR model, or from including waning immunity, the SIRS model. In this section, we provide the system and optimal control differential equations for several variations of the SIR model. We omit the derivations and theory that are similar to the results presented in the previous section for the SIR model. In the section on numerical results we contrast the models and isolate the parameters and terms that cause the most variation in optimal control mitigation strategies.
5.1. The SEIR model. Adding an incubation period to the disease creates the "exposed" class E of subjects who have be exposed to the disease and will ultimately become infectious after a waiting period of 1 ǫ . An obvious consequence of the addition of a waiting period is the increased window of time to vaccinate subjects before the first big wave of infectious individuals is released into the population. We will examine other qualitative model differences in Section 6. The state equations for this model are For comparison purposes, we seek to minimize the same objective functional given in equation (5), which we repeat here for convenience: Following Theorem 4.1, we consider the Hamiltonian equation given by and the optimal controls are found by solving equations (16) - (19) with initial conditions S 0 , E 0 , I 0 , and R 0 along with the adjoint equations

HOLLY GAFF AND ELSA SCHAEFER
and transversality conditions as given in equation (13) of Theorem 4.1. The optimal control pair satisfies τ * = min max 0, 5.2. The SIRS model. We now consider the consequences when the SIR model is modified to account for waning immunity. We assume in this model that after a time period 1 ω , subjects in the removed class return to the susceptible class. Again, one consequence of this variation is intuitively clear. The susceptible class will likely have a long term positive equilibrium extending the need for vaccination. We discuss this and other significant outcomes in the consideration of this model in Section 6.
The state equations for the SIRS model are The objective functional is as in equation (5) of Section 2. The Hamiltonian is formed as before, resulting in the adjoint equations The transversality conditions are as in equation (13) of Section 2, and the optimal control pair will satisfy the same condition given in equations (14) in Theorem 4.1.
6. Numerical results. The optimal control problem was solved using an iterative scheme derived by Hackbush [11] and described briefly here. Using an initial guess for the control variables, ν and τ , the state variables, S, E (if applicable), I, R, are solved forward in time, and then the adjoint variables, λ i , are solved backwards in time. If the new values of the state and adjoint variables differ from the previous values, the new values are used to update ν and τ and the process is repeated until the system converges. For the SIR model, the state equations are given by equations (1) -(3) with the initial conditions given by (4). The adjoint equations are given by equations (10) -(12) with the final time conditions given by (13). Finally, the control variables are given by equations (14). Similarly, the SEIR model has state equations (16) - (19), adjoint equations (20) - (23), and control variables (24) - (25). Finally, the SIRS model has state equations (27) - (29), adjoint equations (30) - (32), and control variables (13). For each run, the control variables are set to an initial guess of 0.0 for all time. In the numerical runs presented here we use the parameter assumptions listed in Table 1.
In sensitivity analyses, we found that the weights in the objective functional (5) had very little qualitative impact on the suggested vaccination scheme, though as we mentioned in Section 2 the inclusion of the term ( R K ) 10 in the functional was essential to the outcome. Likewise, the carrying capacity K can be varied widely provided that it is not extremely small with respect to the number of initially infectious individuals. The incubation period ǫ specified in the SEIR model does affect the number of infectious individuals in the population at any time -the longer the incubation period the more infectives -but the dependence is smooth and the qualitative behaviors the same. Likewise, varying the waning rate ω has the expected impact in the SIRS model -the higher the waning rate, the higher the number of newly susceptibles, and thus the higher the resulting maintenance vaccination level. On the other hand, the models are sensitive to changes in the incidence rate β and the death rate δ.
In our numerical runs, we compare optimal mitigation strategies as we vary incidence rates, death rates, and weights for the relative costs of vaccination and treatment. With a death rate of 0 we varied the incidence rates from 0.05 to 0.55, and for an incidence rate of 0.3 we varied the death rate from 0 to 0.1. Recalling that a factor of 1000 is needed to balance units for treatment and vaccination with the number of infected individuals in the objective functional, we shifted a total weight of 2000 between treatment and vaccination costs. To investigate the required effort for the predicted optimal control, we multiplied the controls by the corresponding states. For treatment, we multiplied the treatment rate by the number of infected individuals in the population at that time to give us an estimate for how much effort this represented. For vaccination, we multiplied the vaccination rate by the number of susceptible individuals in the population for a similar estimate of effort. Our graphs show by color the number of individuals who are infected, treated, or vaccinated in any time step for the various scenarios.
The advice suggested by the various runs is consistent, as is detailed below and illustrated in Figures 1 through 4. Regardless of the relative costs of vaccination and treatment and also regardless of the disease structure, it is essential for the effective control of an epidemic to vaccinate at the highest possible rate as early as possible to minimize the number of individuals who become infectious and to minimize the effort of controlling the epidemic. However, if the relative costs of treatment and vaccination allow treatment to coincide with vaccination, we see that the presence of treatment reduces the size of the infected population correspondingly. While it is not evident in the time scale pictured for this article, we further observe in the SIRS case the need for continued vaccination for small portions of the large susceptible population, as one would expect to see where there is waning immunity.
In Figures 1 and 2 we see the baseline recommendations for vaccination and treatment and the resulting numbers of infected when vaccination and treatment are equally weighted and incident rates and disease death rates vary. Note that the true maximum numbers who are vaccinated or treated at any one time are about  Figure 1. For this figure, vaccination and treatment are weighted equally, the incidence rate is β = 0.3, and the death rate δ varies. All plots show the solution of the model systems using the predicted optimal control solutions for each given value of δ. The graphs in the first column are the results of the SIR model, the second column the SIRS model, and the third column the SEIR model. The first row shows the size of the infected class versus time for each diseaserelated death rate for each model, the second row shows the number vaccinated, and the third shows the number treated. It is important to note that the color scale for the infected class is 0-600, while the scale for the vaccinated and treated is only 0-200. For example, we see that when δ = 0, which is the top-most bar in each plot, there is little variation between the SIR and SIRS model for the size of the infected class. However, as the disease-related death rate increases to 0.1, there is a noted increase in the size of the infected class for SIRS and SEIR models that is not reflected in the SIR model. The optimal number of vaccinations does not vary widely from model to model for varying δ. The timing of treatment does vary from model to model and across values for δ.
550 and 400, respectively, and our smaller scale in the graph is to visually aid with comparisons as time progresses. One can observe in these initial baseline cases that treatment is sometimes foregone to support the necessary immunization; however, there is a corresponding price in the number of individuals that become infected. One can also observe the continuing vaccination in the SIRS cases until the disease is eradicated.  Figure 2. For this figure, vaccination and treatment are weighted equally, the incidence rate β varies, and the death rate δ is zero. All plots show the solution of the model systems using the predicted optimal control solutions for each given value of β. The graphs in the first column are the results of the SIR model, the second column the SIRS model, and the third column the SEIR model. The first row shows the size of the infected class versus time for each disease-related death rate for each model, the second row shows the number vaccinated, and the third shows the number treated.
It is important to note that the color scale for the infected class is 0-600, while the scale for the vaccinated and treated is only 0-200. Similar to the results of Figure 1, the optimal solution shows maximum vaccination across all models and parameter values for at least the first half of the run. Treatment levels, however, vary by model and the incidence values, β. Figure 3 considers the case in which vaccination carries a weight of 200 and treatment carries a weight of 1800, so that treatment is 9 times as expensive as vaccination. This figure should be compared to Figure 2 in which treatment and vaccination are equally weighted. We see that the SIRS population benefits from the less expensive vaccination (observe the missing "bumps" in the second column of Figure 3) because some of the resources saved in the inexpensive vaccination effort can be applied to the necessary treatment. In contrast, the SEIR population suffers because there are not enough resources to treat at the higher rate used in the baseline case for higher incidence rates. The graph shown is for varying incidence, but the results are similar as the death rate varies. The results are also qualitatively  Figure 3. For this figure, treatment was nine times more expensive than vaccination (vaccination weight = 200, treatment weight = 1800), the incidence rate β varies, and the death rate δ is zero. All plots show the solution of the model systems using the predicted optimal control solutions for each given value of β. The graphs in the first column are the results of the SIR model, the second column the SIRS model, and the third column the SEIR model. The first row shows the size of the infected class versus time for each disease-related death rate for each model, the second row shows the number vaccinated, and the third shows the number treated.
It is important to note that the color scale for the infected class is 0-600, while the scale for the vaccinated and treated is only 0-200. Note the increased size of the infected class for the SEIR model as compared with that in Figure 2 as a function of increased cost for treatment.
the same when treatment is only 3 times as expensive or as much as 200 times more expensive. Finally, Figure 4 shows the results as we vary incidence when vaccination is 9 times as expensive as treatment. Again, this figure should be compared to Figure  2 in which treatment and vaccination have equal cost. Ironically, even though treatment is relatively cheap we have a big hole in the numbers treated for the SIRS case, reflecting the urgency of vaccination, and yet the lack of treatment produces a surge of infected individuals. As we would expect from the previous example, the SEIR audience fairs better because the less expensive treatment has a  Figure 4. For this figure, vaccination was nine times more expensive than treatment (vaccination weight = 200, treatment weight = 1800), the incidence rate β varies, and the death rate δ is zero. All plots show the solution of the model systems using the predicted optimal control solutions for each given value of β. The graphs in the first column are the results of the SIR model, the second column the SIRS model, and the third column the SEIR model. The first row shows the size of the infected class versus time for each disease-related death rate for each model, the second row shows the number vaccinated, and the third shows the number treated. It is important to note that the color scale for the infected class is 0-600, while the scale for the vaccinated and treated is only 0-200. Again, the SEIR infected class shows the most dramatic impact of the cost of treatment in comparison to Figure 2.
strong effect on the number who become infected when there is an exposed period for the disease. Finally, Figure 4 shows the results as we vary incidence when vaccination is 9 times as expensive as treatment. Again, this figure should be compared to Figure  2 in which treatment and vaccination have equal cost. Ironically, even though treatment is relatively cheap we have a big hole in the numbers treated for the SIRS case, reflecting the urgency of vaccination, and yet the lack of treatment produces a surge of infected individuals. As we would expect from the previous example, the SEIR audience fairs better because the less expensive treatment has a strong effect on the number who become infected when there is an exposed period for the disease.
Another view of the mitigation strategies suggested by the optimal control simulations can be seen in Figure 5. In this figure, we observe the number of individuals vaccinated and treated on top, and the rates of vaccination and treatment on the bottom. We see that it can be misleading to only record only the recommended rates because at the end of the time period the numbers of infected are so low that a high treatment rate corresponds to a small number of infected individuals being treated. We note that the strategy requires vaccinating at the maximum rate allowed, and when that translates to fewer succeptible individuals remaining and being vaccinated, then additional healthcare dollars can be devoted to treatment while maintaining the recommended vaccination rate. This time-series view of a given scenario helps to illustrate the optimal control recommendation of maximum vaccination as a primary strategy.   Figure 3 where the vaccination weight = 200, treatment weight = 1800, incidence rate =0.35 and the death rate as zero. The top time series shows the predicted numbers of individuals vaccinated (solid line) and treated (dotted line) in the first 25 time steps. In contrast the lower time series shows the predicted optimal control rates for vaccination (solid) and treatment (dotted). These plots show that a prediction of a constant vaccination or treatment rate does not imply a constant level of effort, but rather the number of individuals affected provides a better view of the true costs.
7. Conclusion. While underlying epidemic structure is crucial to answering many questions, it appears that for the range of models we studied there is only subtle variation in the optimal strategies. Our conclusion from our simulations is that vaccination, if available, is an essential tool in fighting an epidemic. However, if resources allow for the provision of treatment as well, this additional tool is a valuable resource in decreasing the number of individuals who are affected by an epidemic, particularly if a disease has an exposed period and/or a high incidence rate.
We recall that N = S + I + R, that S, I, and R are nonnegative, and that N > 0. Thus, 0 ≤ I(N −S) , IS N 2 , S(N −I) N 2 ≤ 1. We recall additionally the bounds 0 ≤ ν ≤ ν max and 0 ≤ τ ≤ τ max , and with this substitution we see that the adjoint system is bounded by linear systems with bounded coefficients. Thus, the sub-and supersolutions are uniformly bounded, establishing bounds for the adjoint system in finite time.
We make the change of variables for each of the six differential equations in our optimality system. For ease of notaiton, we generally omit the dependence on time in the following except in the case taht a specific time is intended. We consider the difference of the resulting equations for s ands, i andī, and so on, and we then simplify the resulting equations by integration with the use of appropriate integrating factors.