key: cord-0687324-3j4mqf5x authors: Kim, Jungeun; Kwon, Hee-Dae; Lee, Jeehyun title: Constrained optimal control applied to vaccination for influenza date: 2016-01-21 journal: Comput Math Appl DOI: 10.1016/j.camwa.2015.12.044 sha: 89e5e252b3311d452f628ba23f5bc2f65f275af5 doc_id: 687324 cord_uid: 3j4mqf5x The efficient time schedule and prioritization of vaccine supplies are important in mitigating impact of an influenza pandemic. In practice, there are restrictions associated with limited vaccination coverage and the maximum daily vaccine administration. We extend previous work on optimal control for influenza to reflect these realistic restrictions using mixed constraints on state and control variables. An optimal control problem is formulated with the aim of minimizing the number of infected individuals while considering intervention costs. Time-dependent vaccination is computed and analysed using a model incorporating heterogeneity in population structure under different settings of transmissibility levels, vaccine coverages, and time delays. An epidemic outbreak causes a crisis in public health accompanying social fear as well as a direct loss due to the disease. In developing a response plan to minimize socio-economic losses in disease outbreak, it is important to predict the transmission dynamics of an infectious disease, to compare the effects of different control strategies, and to design the best strategy. The World Health Organization(WHO) issued warning of pandemic influenza and encouraged to prepare countermeasures in 1999, but it did not draw much attention until the emergence of Severe Acute Respiratory Syndrome(SARS) in 2002-2003. Since the first SARS case was reported from Canton, China on December 2002, more than 31 countries worldwide reported 7956 confirmed cases including 666 fatal cases to the WHO by the end of 21 May, 2003. The SARS aroused great interest in the use of mathematical models to predict the course of the epidemic and to evaluate various management strategies [1] . This interest triggered studies for investigating the dynamics and control of infectious diseases and has been reinforced by the treatment of a pandemic including the influenza A(H1N1) 2009 [2] . Vaccination is the principal control measure for reducing the spread of many infectious diseases. The optimization of vaccination policies before their implementation is essential to better allocate resources and to minimize disease burdens. Recent advances have been made in optimizing vaccine distribution policies in specific settings [3] [4] [5] . Several studies have applied optimal control techniques to determine immunization strategies [6, 7] . However, it is still complicated and controversial to address the best possible strategies because the dynamics of epidemics and the optimal use of vaccines depend on various factors including the structure of the population, vaccine availability, and transmissibility levels. In this context, we take into account of the size and heterogeneous dynamics of groups, different transmission levels, the total and daily maximum amount of vaccines available, and time delays in vaccine implementation. In this paper, we seek optimal time-dependent vaccination strategies within the vaccine availability. We begin by introducing a deterministic, compartmental model of influenza transmission incorporating the structure of the population with different dynamics. We then present the formulation of the optimal control problem to minimize the incidence while satisfying the constraints of the total and daily maximum amount of vaccines available. We use the penalty method to approximate this constrained optimization problem and derive an optimality system that characterizes the optimal control. Finally, we present the results of numerical simulations under various settings and conclude with a summary. The influenza model, SEIAR, is an extension of the standard SEIR model incorporating asymptomatic compartment [8] . In the SEIAR model, infected individuals in the exposed stage can either develop symptoms and move to an infective stage or develop no symptoms and move to an asymptomatic infective stage. This baseline model is modified to include control measures of vaccination and antiviral treatment. The nonlinear system of ODEs describing the influenza dynamics is given by with Λ(t) = ϵE(t) + (1 − q)I(t) + δA(t) and the initial conditions The model classifies individuals into five key compartments of susceptible (S), exposed (E), symptomatic infective (I), asymptomatic infective (A) and removed (R). The number of contact events sufficient for transmitting an infection is βN by mass action incidence, where N is the total population size. A fraction p of individuals in the exposed stage proceeds to the infective stage at the rate κ and the remainder goes to the asymptomatic infective stage also at the rate κ. Exposed individuals are assumed to reduce infectivity by a factor of ϵ, with 0 ≤ ϵ ≤ 1. The compartment E represents the latent stage when ϵ = 0 and the initial asymptomatic and mildly infectious stage when ϵ > 0. Infective members leave the compartment at the rate α with a fraction f recovering from the disease, whereas the rest dying of infection. On average infective individuals have their contact rate reduced by a factor of q. Asymptomatic members have their infectivity reduced by a factor of δ, with 0 ≤ δ ≤ 1, and progress to the removed compartment at the rate η. The time dependent control function ν(t) measures the rate at which susceptible individuals are vaccinated with vaccine efficacy ψ and the infective individuals are treated at the rate τ during the epidemic period. In developing response plans for disease outbreaks, one seeks strategies that can minimize the incidence and/or disease related mortality while considering the cost of intervention strategies. The goal is to minimize the number of people who become infected at a minimal efforts of vaccination. Thus, the objective functional is given by In general, vaccination coverage and the maximum daily vaccine administration are limited during an epidemic. In our optimal control formulation, realistic restrictions associated with vaccination are incorporated using state variable inequality constraints. This can be stated where ν max is the maximum daily vaccination and ν total is vaccine coverage. Now, we are faced with the problem to minimize the number of infected individuals using a limited total vaccination, or stated mathematically minimize J(ν) subject to (2.1) and (2.2). (2. 3) The necessary conditions of the optimal solutions of problem (2.3) cannot be directly derived from Pontryagin's Maximum Principle due to the constraint on vaccination coverage. We can convert this problem to a more convenient form by introducing an auxiliary state variable Then, it follows, Note that the inequality constraint z(T ) ≤ ν total is equivalent to z(t) ≤ ν total , ∀t due to the monotonically increasing property of z(t). The existence of an optimal solution can be established by verifying conditions of the Filippov-Cesari existence theorem [9, 10] . The detailed results are shown in Appendix. Theorem 2.1. There exists an optimal solution to the problem (2.5). One popular approach to handle the constrained optimal control problems is to convert them into unconstrained optimization problems. This can be accomplished by using a penalty function method, which attempts to attain feasibility and to optimize the objective functional simultaneously by optimizing a weighted combination (refer to [11] for details). A suitable penalty function of the constrained minimization problem (2.5) is where P and Q are constants that can be chosen to balance the relative costs of the infected individuals and vaccination. In summary, we wish to approximate the solution of constrained optimization problem (2.5) by solving Pontryagin's Maximum Principle is used to derive the optimality system which provides necessary conditions of the optimal solutions of (2.7). The optimality system consists of the state system (2.8) with initial conditions, the adjoint system (2.9) with the transversality conditions, and the optimality conditions (2.10). The details are given in Appendix. (2.8) with Λ(t) = ϵE(t) + (1 − q)I(t) + δA(t) and the initial conditions (2.10) The dynamics of epidemics and the effectiveness of an intervention strategy depend on the structure of the population. For example, the potential benefit of prioritizing school-age children for vaccination has been discussed because this group is disproportionately responsible for influenza transmission [12, 13] . Thus, the SEIAR model of an epidemic is extended to capture characteristics regarding age and severity of disease. To incorporate these features, we divide the population into subpopulations with different transmission dynamics. where a ij is the frequency of contact between an individual in the jth subgroup and individuals of the ith subgroup. The objective functional to be minimized is given by and restrictions associated with vaccination coverage and the maximum daily vaccine administration are incorporated using state variable inequality constraints: Therefore, we seek the optimal controls ν = [ν 1 (t), . . . , ν m (t)] T such that minimize F (ν) subject to (2.11) and (2.12). (2.13) To convert this problem to a more convenient form from which optimal solutions can be approximated, we introduce an auxiliary state variable and a penalty functional where P i and Q i represent the weight constants of infective individuals and controls. µ 1 and µ 2 are penalty parameters which balance objective functional and feasibility. H 1 and H 2 denote the Heaviside step functions, that is The derivation of optimality system is similar to the SEIAR model. In summary, we find the optimal controls by solving the state equations with initial conditions and the adjoint equations The optimal control functions satisfy We present simulation results by solving the constrained optimal control problems using penalty method as described in the previous chapter. Applying the Pontryagin's Maximum Principle, optimality system is derived from which optimal solutions may be obtained. This system is a two-point boundary value problem, where the state system with the initial conditions (2.8) and the adjoint system with the terminal conditions (2.9) are coupled. Among many practical approaches to solve the coupled optimality system, we use a gradient-based algorithm to uncouple the system [5] . In the numerical simulations, we compare the optimal intervention strategies under different settings of the population structures, transmissibility levels, vaccine coverages, and delays in vaccine production. The values of model parameters were obtained by the best available local data and literature reviews of previous studies [8, 14] . Table 1 provides a summary of the definitions and values of the parameters. The default values for parameters in addition to those listed in Table 1 are vaccination coverage, maximum daily vaccination, and the basic reproduction number. In general, vaccination coverage for a pandemic influenza is far lower than full, for example, the US or Canada held 30%-40% coverage for 2009 H1N1 pandemic influenza [15, 16] . The daily rate of vaccine administration has been estimated to be below 2% of the total population [16, 17] . Therefore, we varied vaccination coverage in the range of 10%-50% and maximum daily administration of less than 2% of the total population. The basic reproduction number R 0 , which associated with transmissibility level, ranges from 1.4 to 2.4. The weight constants in the penalized objective functional (2.14) are P i = 1, Q i = 1, µ 1 = 1 and µ 2 = 1. For simplicity and straightforward analysis of causality, we divide the population into two subgroups with the ratio of sizes r 1 : r 2 , and start simulations in early stage of epidemic by introducing one infective to the whole population of susceptible. In other words, we take initial conditions S 1 (0) = r 1 r 1 +r 2 × 5 × 10 7 , S 2 (0) = r 2 r 1 +r 2 × 5 × 10 7 , I i (0) = 1, and All the baseline values introduced here are used throughout the paper unless otherwise specified. Results from simulations are presented in terms of number of vaccine doses, proportion vaccinated, number of infectives, and reduction in infectives. The proportion vaccinated is the number of vaccinated individuals in each subgroup divided by the corresponding population size of each subgroup. The reduction in infectives, the number of infectives reduced as a result of control relative to the number of infectives in the absence of vaccination, assesses the effectiveness of control strategies. The optimal control is compared with early possible vaccination and uniform distribution. In order to investigate whether each factor of population structure plays an important role in determining efficient vaccination strategy, we perform numerical experiments under different settings of contact rates and subgroup sizes. We assume a 50% vaccination coverage, a 1.5% maximum daily administration of the total population and R 0 = 1.9. First, we consider two subgroups of the fixed ratio 3:7 of sizes where the former subgroup has a higher contact than the latter. The mixing rates are varied by using different frequencies of contacts a ij between jth and ith subgroups. Fig. 2A2-C2 . Under the optimal vaccination, control efforts are gradually decreasing in time. Vaccination is highly concentrated in the early stage of epidemic for early possible strategy and the vaccinated proportion stays constant over the time for uniform distribution. A larger fraction is vaccinated in the higher contact group compared to the lower contact group under the optimal strategy. In contrast, the number of individuals vaccinated in each subgroup is proportional to the size of subgroup for both early possible and uniform strategy, or the proportion vaccinated in each subgroup is kept the same. This is because neither of them reflects the heterogeneous mixing rates to determine the vaccination schedule. Fig. 2A3 -C3 demonstrate that the implementation of control measures yields substantial reductions in the peak size of infections while generating shifted, longer pandemic durations. The reduction of cumulative number of infectives relative to the case without intervention are approximately 68.2%, 59%, 45.5% under optimal, early possible, and uniform distribution, respectively (Fig. 2A4-C4 ). This result proposes that control measures may be optimally targeted to minimize the number of cases for a given control effort. The impact of contact rates on optimal vaccination is evaluated by varying a 11 with fixed a 21 = 1 and a 22 = 1 (Fig. 3) . The results for the contact of a 11 = 7 and a 11 = 10 are similar to case of a 11 = 5, in general. Our findings suggest that more intensive vaccination must be implemented in the early stage on the targeted group as the gap of contact rates between two subgroups grows (Fig. 3A2-C2) . As a result, for a 11 = 10, infectives are significantly reduced by 80.1% compared to early possible vaccination by 51.3% and uniform distribution by 31.8% (Fig. 3C4) . We may conclude that differences in efficiency between optimally targeted control and others in which heterogeneous transmission dynamics is not considered, is growing as the contact frequency differences increase. Then we explore the impact of the population sizes of subgroups on optimal targeting strategy. The proportions of vaccinated in each group using the population ratios of 5:5 and 2:8 are illustrated in Fig. 4 where a 11 = 10 and a 21 = 1 are fixed and a 22 is varied by using 1, 3, 5 and 10. It is obvious that the proportion of vaccinated individuals for each subgroup corresponds to contact frequencies when the subgroups are equal in population sizes ( Fig. 4A1-A4 ). However, a switch in the main targeted group occurs as the gap of contact rates between two subgroups diminishes ( Fig. 4B1-B4 ). Similar behaviour is observed with different ratios of the subgroup population sizes. Optimal targeting strategy is investigated in terms of cumulative proportion of vaccinated in each group under different ratios of populations sizes. The ratio of the vaccinated proportion in the higher contact group relative to the lower contact group as a function of population sizes and contact frequencies is illustrated in Fig. 5(A) . For all cases of subgroup population sizes, the main targeted vaccinated group is the higher contact group when the gap of contacts between two groups is large enough. As the gap decreases, the targeted group is switched to the one with more population and the bigger the difference in sizes of two groups, the earlier the turning point comes. This threshold behaviour is shown in Fig. 5(B) . We study the effects of transmission levels on optimal vaccination strategies using different values for basic reproduction number R 0 . In Fig. 6 , optimal vaccination controls and their impacts on the incidence of infected are compared under three different values of R 0 , 1.4, 1.9, and 2.4. Fig. 6A1 -C1 display the time series of vaccinated proportion for each subgroup. Vaccination must be implemented in a shorter period of time, for higher level of transmission because large R 0 results in earlier epidemic peaks due to rapid spread of the disease. It is also observed that high R 0 requires optimal control that manages to reduce the infectives in the group with big population since large R 0 quickly depletes the susceptible population. Incidence curves of infected proportions are shown in Fig. 6A2-C2 and higher values of R 0 generate earlier epidemic peaks with higher and narrower in shape. For all cases, optimal vaccination yields more reduction in infected than early possible and uniform distribution (Fig. 6A3-C3) . In developing response plans for disease outbreaks, vaccination coverage and maximum daily vaccination are parameter values with restriction, in practice. It is of great interest to seek the best strategies to better allocate limited resources and minimize disease burdens. In the optimal control problem, vaccination coverage and maximum daily vaccination are expressed as mixed constraints on control and state to reflect realistic restrictions. We study the effects of optimal vaccination To evaluate the impact of constraints on the availability of vaccine supplies, we vary vaccination coverage from 50% to 10% of the population, fixing maximum daily vaccination as 1.5% when R 0 = 1.9. Time schedules of vaccinated proportion under three different vaccination coverage are monotonically decreasing in time with different slopes in Fig. 7A1-C1 . Cumulative vaccinated proportion in each subgroup is shown in Fig. 8(A) , which demonstrates that targeting strategy remains similar with changes in vaccination coverage. As vaccination coverage decreases, number of vaccinated decreases, resulting in an increase in the overall number of infected (Fig. 7A2-C2) . Optimal control is more efficient than other strategies in terms of reduction under all vaccination coverage (Fig. 7A3-C3) . A significant difference is observed in the reduction of number of infected under different vaccination coverage and this result underscores the importance to secure sufficient vaccine supplies (Fig. 8(B) ). The optimal value of maximum daily vaccination is related to total vaccination coverage and the transmission level. This relation is illustrated in Figs. 9 and 10 under three different values of R 0 . Fig. 9 displays cumulative number of vaccinated as maximum daily vaccination varies from 0.25% to 1.5% with 50% vaccination coverage. The total number of vaccinated increases as daily maximum increases to 1% and 1% is enough to use 50% coverage in total. We also analyse the threshold For instance, in Fig. 10 R 0 = 1.9, the reduction of infectives increases as daily maximum increases to 1% for vaccination coverage of 50%. However, further increase in maximum daily vaccination does not result in more reduction of infected individuals. Our results indicate that maximum daily vaccination affects the optimal control or the dynamics of influenza only in a certain range which depends on vaccination coverage and R 0 . In the previous simulations, we found that sufficient vaccination should be implemented in the early phases to effectively control an epidemic. In practice, however, the timing of vaccine delivery or availability may be much later than the onset of the pandemic. For the pandemic A(H1N1) 2009 outbreak, the influenza began spreading in April 2009 and vaccination started in October 2009. Therefore, we incorporated the time delay by varying the start of vaccination, 30, 60, and 90 days after the pandemic onset. We investigate the influence of the timing of the vaccination on optimal strategy and vaccination coverage. Fig. 11 displays the vaccinated proportions of subgroups at high and low contact, corresponding incidence curves of infected, and the reduction in the cumulative number of infected individuals relative to the ones in the absence of interventions. As the time delay increases, more intensive vaccination must be implemented in a shorter period of time as shown in Fig. 11A1 -D1. The group with higher contact remains as the main target of optimal vaccine allocation for most delay cases. However, a longer period of time delay leads to the switch in the vaccinated proportion (Fig. 11D1 ). This is because peak occurs earlier in higher contact group than lower contact group and vaccination has little effect after the peak. Delay of vaccination yields significant increases in the number of infected individuals and decreases in the percent reductions ( Fig. 11A2-D2 A3-D3) . While general performance of optimal strategy is superior to others, the difference becomes insignificant for a longer period of time delay (Fig. 11A3-D3) . Fig. 12 presents the cumulative number of vaccinated of each group under three different vaccination coverages and the impact of optimal vaccinations in terms of reductions of infected. As vaccination coverage increases, the number of infected individuals decreases, resulting in an increase in the reductions of infected. However, with long time delay and large vaccination coverage, the whole coverage is not distributed due to the lack of implementation period. For instance, if vaccination starts 60 days after the beginning of the transmission, only 39% coverage is allocated when 50% vaccination coverage is allowed. Not surprisingly, 50% coverage and 25% coverage lead to the same results in cumulative number of vaccinated and the reductions of infectives with 90 days of delay. The optimal vaccination coverage as a function of time delay is presented in Fig. 13 . Vaccination is among the most important control measures for reducing the spread of many infectious diseases. Thus, it is great interest to develop an efficient time schedule and prioritization of limited vaccine supplies. This study uses a mathematical model of the transmission dynamics of pandemic influenza and employs techniques from control theory to derive optimal intervention strategies. Time-dependent vaccination is computed and analysed using a model incorporating heterogeneity in population structure under different settings of transmissibility levels, vaccine coverages, and time delays in the implementation of vaccination. Vaccination coverage and the maximum daily vaccination rate are critical factors in developing response plans for disease outbreaks. These are parameters with limited values and associate with weight constants and a control upper bound in optimal control formulation of previous studies. In our mathematical framework, mixed constrained optimization problem is considered to reflect realistic restrictions. The impact of heterogeneous contact rates on optimal allocation is evaluated and compared with other strategies not incorporating population structure. Our analysis confirms that more intensive control measures must be implemented in the early stage on the group with higher contact, in general. However, the targeted group is switched to the one with more population, as the gap of contact rates between two groups diminishes, which suggests that prioritization strategy should be tailored. Higher values of R 0 requires rapid implementation of optimal control policies with more concern on the group with big population. While the reduction of infected significantly decreases as the vaccination coverage decreases, prioritization schemes on subgroups remain similar with changes in total number of vaccine doses available. Finally, the benefit of vaccination is significantly reduced as the time of the start of vaccination from pandemic onset is delayed. Consider the following optimal control problem: and The existence of a solution to the optimal control problem can be obtained by verifying conditions of the Filippov-Cesari existence theorem [9] . The boundedness of solutions to the system (A.1) for the finite time interval is needed to establish these conditions. Note that the quantities S, E, I, and A decrease only in proportion to their present sizes, respectively, and thus, all variables remain nonnegative if the initial values are nonnegative. Then so are R and z, since the changes in these variables are nonnegative. To establish the upper bounds for the solutions, we consider an equation for the total population size N. N satisfies N ′ = −(1 − f )αI from the system (A.1) and N is bounded above by N(0). Because S, E, I, A, and R are all nonnegative, the upper bound for N is also the upper bound for S, E, I, A, and R. Boundedness of z follows from the boundedness of ν and S. Let x(t) ≡ (x 1 (t), . . . , x n (t)) ∈ R n be a state vector and u(t) ≡ (u 1 (t), . . . , u r (t)) ∈ R r be a control vector. Consider the following optimal control problem: subject tȯ the terminal conditions and the constraints Assume that the functions F : R n × R r × R → R, f : R n × R r × R → R n and g : R n × R r × R → R s are continuously differentiable with respect to all their arguments. We call (x(t), u(t)) an admissible pair if u(t) is any piecewise continuous control and x(t) is a continuously differentiable function such that (A.4)-(A.6) are satisfied. Theorem A.1 (Filippov-Cesari's Existence Theorem). Suppose that there exists an admissible pair (x(t), u(t)) and further that 1. U is closed. 2. N(x, t) = {ỹ ≡ (y, y n+1 ) : y = f (x, u, t), y n+1 ≥ F (x, u, t), g(x, u, t) ≥ 0, u ∈ U} is convex for all (x, t) ∈ R n × [t 0 , t 1 ]. 3. There exists a number θ > 0 such that ∥x(t)∥ < θ for all admissible pairs (x(t), u(t)), and all t ∈ [t 0 , t 1 ]. There exists an open ball B(0, γ ) ⊂ R r which contains the set Ω(x, t) = {u ∈ U : g(x, u, t) ≥ 0} for all x ∈ B(0, θ ). Then there exists an optimal pair (x * (t), u * (t)) to the problem (A.3)-(A.6) with u * (t) measurable. For a proof of theorem 5.1, see Cesari [9] . Then, we may differentiate the Lagrangian L with respect to ν to obtain ∂L ∂ν = 2Q ν(t) + 2µ 1 (ν(t)S(t) − ν max )S(t)H 1 (ν(t)S(t) − ν max ) − (λ 1 (t) − λ 5 (t) − λ 6 (t))S(t) − ω 1 (t) + ω 2 (t) = 0. Solving for the optimal control, we obtain ν * (t) = 2µ 1 ν max S(t)H 1 + (λ 1 (t) − λ 5 (t) − λ 6 (t))S(t) + ω 1 (t) − ω 2 (t) To determine an explicit expression for the optimal control without the penalty multipliers ω 1 and ω 2 , we consider the following three cases: 1. On the set {t|0 < ν * (t) < 1}, we have ω 1 (t) = ω 2 (t) = 0. Hence the optimal control is ν * (t) = 2µ 1 ν max S(t)H 1 + (λ 1 (t) − λ 5 (t) − λ 6 (t))S(t) 2  Q + µ 1 S 2 (t)H 1  . 2. On the set {t|ν * (t) = 1}, we have ω 1 (t) = 0 and ω 2 (t) ≥ 0. Hence 1 = ν * (t) = 2µ 1 ν max S(t)H 1 + (λ 1 (t) − λ 5 (t) − λ 6 (t))S(t) − ω 2 (t) which implies that ν * (t) = 2µ 1 ν max S(t)H 1 + (λ 1 (t) − λ 5 (t) − λ 6 (t))S(t) 3. On the set {t|ν * (t) = 0}, we have ω 2 (t) = 0 and ω 1 (t) ≥ 0. Hence ν * (t) = 2µ 1 ν max S(t)H 1 + (λ 1 (t) − λ 5 (t) − λ 6 (t))S(t) + ω 1 (t) which implies that ν * (t) = 2µ 1 ν max S(t)H 1 + (λ 1 (t) − λ 5 (t) − λ 6 (t))S(t) Combining these three cases, the optimal control ν * is characterized as ν * (t) = min  1, max  0, 2µ 1 ν max S(t)H 1 + (λ 1 (t) − λ 5 (t) − λ 6 (t))S(t) For the structured model, we obtain the optimality system with the similar arguments beginning by defining the Lagrangian as follows: where ω ij (t) ≥ 0 are the penalty multipliers satisfying ω i1 (t)ν i (t) = ω i2 (t)(1 − ν i (t)) = 0 at ν i = ν * i . Here ν * i is the optimal control pair yet to be found. Transmission dynamics and control of severe acute respiratory syndrome Pandemic potential of a strain of influenza A (H1N1): early findings Use of an inactivated vaccine in mitigating pandemic influenza A(H1N1) spread: a modelling study to assess the impact of vaccination timing and prioritisation strategies Optimal Control Theory with Economic Applications Perspectives in Flow Control and Optimization Modeling optimal age-specific vaccination strategies against pandemic influenza Optimal control of an influenza model with seasonal forcing and age-dependent transmission rates A model for influenza with vaccination and antiviral treatment Optimization Theory and Applications: Problems with Ordinary Differential Eqttations The Mathematical Theory of Optimal Processes Linear and Nonlinear Programming Who should get influenza vaccine when not all can? Prevention Effectiveness: A Guide to Decision Analysis and Economic Evaluation Containing pandemic influenza with antiviral agents Public Health Agency of Canada Peterborough County-city health unit pandemic influenza plan, Annex A: Mass vaccination plan The work of Hee-Dae Kwon was supported by the NRF Grant funded by the Korean government (NRF-2014R1A1A2 056498). The work of Jeehyun Lee was supported by NRF grant 2015 R1A5A1009350. Proof of Theorem 2.1. We verify nontrivial requirements listed in the Filippov-Cesari's existence theorem. In order to verify these conditions, we write x(t) = (S(t), E(t), I(t), A(t), R(t), z(t)) T , u(t) = ν(t), F (x(t), u(t), t) = Px 3 (t) + Qu 2 (t),It is clear that F , f and g are of class C 1 and f is bounded. Due to this property, there exists a solution for (2.8) for a zero constant control, which guarantees an admissible pair (x(t), u(t)). Conditions 1 and 4 hold vacuously as the control setWe note that U is convex, f is expressed as a linear function of the control variable with coefficients dependent on the state variables and time in (A.7), and F and g are convex on U. Then N(x, t) is convex, which is Condition 2. To show this, letThus,ρ ∈ N, which shows that N is convex. Condition 3 follows from the boundedness of solutions to the system (2.8) for the finite time interval. The Lagrangian, the Hamiltonian augmented with penalty terms for the control constraints is, given by, where ω i (t) ≥ 0 are the penalty multipliers satisfyingHere ν * is the optimal control pair yet to be found. Setting to zero the first variations with respect to state variables, S, E, I, A, R and z, yield the adjoint equations= βδS(t) (λ 1 (t) − λ 2 (t)) + η (λ 4 (t) − λ 5 (t)) λ ′ 5 (t) = 0 λ ′ 6 (t) = −2µ 2 (z(t) − ν total ) H 2 (z(t) − ν total ), with the transversality conditions λ 1 (T ) = λ 2 (T ) = λ 3 (T ) = λ 4 (T ) = λ 5 (T ) = λ 6 (T ) = 0.