Constrained optimal control applied to vaccination for influenza.

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.

Introduction

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.

Mathematical models

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 and the initial conditions Fig. 1 shows a flow diagram for model (2.1).

Fig. 1

Flow chart for the SEIAR model.

The model classifies individuals into five key compartments of susceptible (), exposed (), symptomatic infective (), asymptomatic infective () and removed (). The number of contact events sufficient for transmitting an infection is by mass action incidence, where is the total population size. A fraction 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 . The compartment represents the latent stage when and the initial asymptomatic and mildly infectious stage when . Infective members leave the compartment at the rate with a fraction recovering from the disease, whereas the rest dying of infection. On average infective individuals have their contact rate reduced by a factor of . Asymptomatic members have their infectivity reduced by a factor of , with , and progress to the removed compartment at the rate . The time dependent control function 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 is the maximum daily vaccination and 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

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,

Therefore, (2.3) is transformed into Note that the inequality constraint is equivalent to due to the monotonically increasing property of .

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.

There exists an optimal solution to the problem   (2.5) .

Penalty method

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 and are constants that can be chosen to balance the relative costs of the infected individuals and vaccination. and are penalty parameters which change the relative severity of violation of the constraints. and denote the Heaviside step functions, i.e.

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. with and the initial conditions . with the transversality conditions .

Structured model

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. with and the initial conditions where is the frequency of contact between an individual in the th subgroup and individuals of the th 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 such that

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 and represent the weight constants of infective individuals and controls. and are penalty parameters which balance objective functional and feasibility. and denote the Heaviside step functions, that is

Now, the problem (2.13) is transformed into

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 with the transversality conditions .

The optimal control functions satisfy

Numerical simulations

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 , which associated with transmissibility level, ranges from 1.4 to 2.4. The weight constants in the penalized objective functional (2.14) are , , and .

Table 1

Parameters in numerical simulations.

Parameter	Description	Value
ϵ	Infectivity reduction factor for the exposed	0
q	Contact reduction by isolation	0.5
δ	Infectivity reduction factor for the asymptomatic	0.5
p	Fraction of developing symptoms	0.667
κ	Transition rate for the exposed	0.7143/day
f	Complement to fatality rate (one minus fatality rate)	0.999
α	Recovery rate for the (symptomatic) infective	0.1667/day
η	Recovery rate for the asymptomatic	0.1667/day
τ	Antiviral treatment rate	0/day
ψ	Efficacy of vaccination	70%

For simplicity and straightforward analysis of causality, we divide the population into two subgroups with the ratio of sizes , 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 , , , and , for . 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.

Population structure

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 . 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 between th and th subgroups. Fig. 2, Fig. 3 show optimal vaccinations under the different contact rates and compare the percent reductions in infectives relative to the baseline scenario with other strategies not incorporating heterogeneous transmission dynamics.

Fig. 2

Time schedule of vaccine allocations, the vaccinated proportions of subgroups, the incidence curves and the reductions in infectives are plotted under optimal vaccine distribution (A), early possible vaccination (B), and uniform distribution (C) for comparison.

Fig. 3

Optimal vaccine allocations, the vaccinated proportions of subgroups, the incidence curves and the reductions in infectives are displayed when (top), (middle) and (bottom).

Time dependent optimal vaccine distribution, early possible vaccination, and uniform distribution over the time interval are displayed in Fig. 2A1–C1 when , and . The cumulative proportion vaccinated per month in each subgroup using different vaccination strategies are depicted in 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 with fixed and (Fig. 3). The results for the contact of and are similar to case of , 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 , 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 and are fixed and 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.

Fig. 4

The vaccinated proportions of subgroups are compared using different ratio of population sizes 5:5 (top) and 2:8 (bottom). The contact frequency and are fixed and is varied by using 1, 3, 5 and 10, respectively.

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).

Fig. 5

The ratio of cumulative proportions of vaccinated in the higher contact group relative to the lower contact group is displayed in A. The threshold value of at which the cumulative proportions of vaccinated in two groups become equal is plotted in B. in the ratio of population sizes varies from 1 to 5 and varies from 1 to 10, fixing and .

Transmission levels and vaccination coverage

We study the effects of transmission levels on optimal vaccination strategies using different values for basic reproduction number . In Fig. 6 , optimal vaccination controls and their impacts on the incidence of infected are compared under three different values of , 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 results in earlier epidemic peaks due to rapid spread of the disease. It is also observed that high requires optimal control that manages to reduce the infectives in the group with big population since large quickly depletes the susceptible population. Incidence curves of infected proportions are shown in Fig. 6A2–C2 and higher values of 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).

Fig. 6

Optimal vaccine allocations, the vaccinated proportions of subgroups, the incidence curves and the reductions in infectives are displayed under different values of : 1.4 (A), 1.9 (B), and 2.4 (C).

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 polices on the dynamics of disease under different vaccination coverage and maximum daily vaccination capacity. In both simulations, the ratio of subgroup population sizes is 3:7, contact frequencies are , and .

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 . Time schedules of vaccinated proportion under three different vaccination coverage are monotonically decreasing in time with different slopes in Fig. 7 A1–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)).

Fig. 7

Time schedules of vaccinated proportion, the incidence curves are evaluated under vaccination coverage of 50%, 25%, and 10%. Optimal vaccination is compared with early possible and uniform distribution in terms of the reductions in infectives.

Fig. 8

Cumulative vaccinated proportion in each subgroup (A) and reduction in number of infectives (B) are shown for three different vaccination coverage of 50%, 25% and 10%.

The optimal value of maximum daily vaccination is related to total vaccination coverage and the transmission level. This relation is illustrated in Fig. 9, Fig. 10 under three different values of . 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 value of maximum daily vaccination as a function of vaccination coverage and in terms of reduction of infected individuals. For instance, in Fig. 10 , 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 .

Fig. 9

Cumulative number of vaccinated is shown as maximum daily vaccination varies from 0.25% to 1.5% under three different values of .

Fig. 10

The reduction of infectives is plotted as a function of maximum daily vaccination and total vaccination coverage under three different values of .

Delay of vaccination

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. 11

Optimal vaccine allocations, the corresponding incidence curves and the reductions in infectives are displayed under the time delay of vaccination by 0 (A), 30 (B), 60 (C) and 90 days (D).

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 .

Fig. 12

Total number of vaccinated in each group using three different vaccination coverages of 50%, 25%, and 10% are presented under the time delay of vaccination by 0, 30, 60 and 90 days (left). The corresponding percent reductions in infected individuals are shown (right).

Fig. 13

Optimal vaccination coverage as a function of delay time is plotted.

Conclusion

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 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.