Optimal Control of the COVID-19 Pandemic with Non-pharmaceutical Interventions.

The COVID-19 pandemic has forced societies across the world to resort to social distancing to slow the spread of the SARS-CoV-2 virus. Due to the economic impacts of social distancing, there is growing desire to relax these measures. To characterize a range of possible strategies for control and to understand their consequences, we performed an optimal control analysis of a mathematical model of SARS-CoV-2 transmission. Given that the pandemic is already underway and controls have already been initiated, we calibrated our model to data from the USA and focused our analysis on optimal controls from May 2020 through December 2021. We found that a major factor that differentiates strategies that prioritize lives saved versus reduced time under control is how quickly control is relaxed once social distancing restrictions expire in May 2020. Strategies that maintain control at a high level until at least summer 2020 allow for tapering of control thereafter and minimal deaths, whereas strategies that relax control in the short term lead to fewer options for control later and a higher likelihood of exceeding hospital capacity. Our results also highlight that the potential scope for controlling COVID-19 until a vaccine is available depends on epidemiological parameters about which there is still considerable uncertainty, including the basic reproduction number and the effectiveness of social distancing. In light of those uncertainties, our results do not constitute a quantitative forecast and instead provide a qualitative portrayal of possible outcomes from alternative approaches to control.

Introduction

The first known case of COVID-19 in the USA arrived from China on January 15, 2020 (Holshue et al. 2020). In the weeks that followed, deficiencies in testing allowed the virus that causes this disease to spread largely undetected (Perkins et al. 2020), with the number of reported cases growing to 429,319 by April 8, 2020 [39]. Social distancing—such as closing schools, working from home, and sheltering in place—has been adopted widely and appears capable of impacting transmission (Kraemer et al. 2020; Cowling et al. 2020). These practices have tremendous economic and social consequences, however, which has precipitated growing desire to relax them [17].

A danger in relaxing the use of social distancing is that the virus could resurge in its absence (Prem et al. 2020). A safe, effective, and widely used vaccine would obviate the need for social distancing, but there are many challenges associated with developing a vaccine against COVID-19 (Lurie et al. 2020) and one is not likely to be available until at least spring 2021 (Amanat and Krammer 2020). Strategies for successfully controlling COVID-19 until then will depend on a suite of non-pharmaceutical interventions (Aledort et al. 2007), including some degree of social distancing but also diagnostic testing, contact tracing, and case isolation (Bedford et al. 2020).

In models of pathogen transmission based on mass-action assumptions, non-pharmaceutical interventions all act in a similar way by reducing the transmission coefficient. This parameter, often denoted , reflects the product of the rate of contact between susceptible and infectious people and the probability of transmission upon contact (Keeling and Rohani 2007). Due to the extent to which they reduce contact, full-scale lockdowns appear to be the most effective strategy for reducing transmission available at present (Flaxman et al. 2020). Until alternative non-pharmaceutical interventions can be implemented that are similarly effective but less disruptive economically, it is urgently important to determine the minimal extent to which contact must be reduced to achieve public health objectives.

Optimal control theory offers a way to understand how to apply one or more time-varying control measures to a nonlinear, dynamical system in such a way that a given objective is optimized (Lenhart and Workman 2007). These techniques have been widely applied to a variety of pathogen transmission systems before (e.g., Blayneh et al. 2010; Miller Neilan et al. 2010; Choi and Jung 2014; Agusto and Khan 2018), including in the context of a pandemic of a respiratory pathogen (e.g., Lin et al. 2010; Tchuenche et al. 2011; Shim 2013). The latter studies indicate that the level of non-pharmaceutical interventions required is dependent on model parameters and that application of non-pharmaceutical interventions can be required at a high level and for a long duration in the absence of a vaccine.

Here, we apply optimal control theory to determine optimal strategies for the implementation of non-pharmaceutical interventions to control COVID-19, with a focus on the USA. An optimal strategy in this sense involves weighing the relative costs of control and COVID-related mortality, and finding an approach to control that minimizes that combined cost. Other optimal control analyses of COVID-19 are beginning to emerge (Djidjou-Demasse 2020; Mallela 2020; Patterson-Lomba 2020; Piguillem and Shi 2020; Shah et al. 2020), although these have been less focused on any particular geographic setting. One contextual feature of the USA that our analysis considers is the level of control that was implemented early in the pandemic. We account for this, because our objective is less about understanding what level of control would have been optimal earlier on and more about understanding what level of control would be optimal later on. We assume that a vaccine is necessary as a long-term solution to COVID, so we oriented our analysis around the time frame between May 2020 and a hypothetical vaccine introduction beginning in March 2021. In addition, we emphasize the sensitivity of our results to epidemiological parameters that are not well characterized but appear influential in determining what range of impacts of non-pharmaceutical interventions on COVID-19 are possible. All code and data used in our analysis are available at http://github.com/TAlexPerkins/covid19optimalcontrol.

Methods

Model

We modeled SARS-CoV-2 transmission according to a system of ordinary differential equations,with variables and parameters defined in Tables 1 and 2.

Table 1

State variables in the model

Symbol	Definition
S	Susceptible individuals
E	Exposed but not yet infectious
A	Asymptomatic infections
I	Symptomatic infections
H	Hospitalized infections
V	Vaccinated individuals who have not been infected
u	Control with non-pharmaceutical interventions

Table 2

Model parameters

Symbol	Definition	Min	Mid	Max	Source
α	Relative infectiousness of asymptomatic infections	0.046	0.602	0.981	Perkins et al. (2020)
β	Transmission coefficient	0.514	0.600	0.686	Park et al. (2020)
γ	Rate of progression through infectious stage	0.125	0.31	0.5	Chinazzi et al. (2020)
δ	Baseline per capita death rate	–	3.18×10^-5	–	Martin et al. (2019)
Δ_-	Minimum probability of death following hospitalization	0.0865	0.104	0.1213	Centers for Disese Control and Prevention (2020)
Δ_+	Maximum probability of death following hospitalization	0.233	0.28	0.328	Onder et al. (2020)
ϵ	Per exposure protection from vaccination	0.5	0.8	0.9	Assumed
η	Rate of progression through hospitalization	0.0485	0.075	0.172	http://virological.org/t/epidemiological-data-from-the-ncov-2019-outbreak-early-descriptions-from-publicly-available-data/337
ι	Importation rate of infected people	1.07×10^-7	3.24×10^-7	9.7×10^-7	Calibrated
κ	Probability of hospitalization among symptomatic infections	0.207	0.26	0.314	Centers for Disese Control and Prevention (2020)
μ	Birth rate	–	3.18×10^-5	–	Martin et al. (2019)
ν	Rate at which vaccines are administered	0.001516	0.00197	0.002467	https://www.cdc.gov/h1n1flu/vaccination/
ρ	Rate of progression to infectiousness following infection	0.166	0.2	0.333	Chinazzi et al. (2020)
σ	Proportion of infections that result in symptoms	0.5	0.82	0.9	Mizumoto et al. (2020)
τ_ν	Day on which vaccination begins	365	456	548	Amanat and Krammer (2020)
ω	Proportion of deaths observed	0.5	–	0.8	Assumed
h	Steepness of increase in death probability when H_max is exceeded	1403	701	350	Assumed
H_max	Maximum hospital capacity	–	2.74×10^-3	–	https://www.medrxiv.org/content/early/2020/03/30/2020.03.27.20043752
u_max	Maximum value of the control	0.5	0.7	0.9	Jarvis et al. (2020)

Parameterization of , , , , , , , , and is explained in further detail in Sect. 2.4. The timescale of all parameters related to time is daily

In the model, all individuals are initially susceptible, S, with infections introduced at a constant rate, , from outside the population. Individuals then transition to the exposed class, E, where they reside for an average of days. A proportion experience a symptomatic infection, I. The remainder experience an asymptomatic infection, A, and have a fraction, , of the infectiousness of symptomatic infections. Individuals reside in the I and A classes for an average of days. All asymptomatic infections and a proportion of symptomatic infections then recover and become fully immune to subsequent infection. The remaining proportion of symptomatic infections transition to the hospitalized class, H, from which they exit through either recovery or death after an average of days.

Rather than track deaths due to COVID-19, D(t), as a state variable in Eq. (1), we assume that they follow directly from hospitalizations, H(t), according toThis results in the probability of death moving beyond a minimum of toward a maximum of as H(t) exceeds , as illustrated in Fig. 1. The motivation for this choice is to account for the possibility that patients could experience increased mortality when the demand for certain resources, such as intensive care unit beds or ventilators, exceeds their availability. One other optimal control analysis of COVID-19 has incorporated a similar phenomenon (Piguillem and Shi 2020), albeit with a different functional form.

Fig. 1

Relationship between the number of hospitalizations and the probability of death from COVID-19 among hospitalized patients. The parameters and represent lower and upper bounds on the probability of death, and represents the hospital capacity above which the probability of death exceeds . Hospitalizations are quantified as a proportion of the overall population

Individuals who are recovered and immune are not followed explicitly, as the model’s assumption of density-dependent transmission only requires specification of susceptible and infectious classes in the transmission term. Due to our assumption that rates of birth, , and death due to reasons other than COVID-19, , are equal, changes in population size over the course of the epidemic are modest, making the distinction between density- and frequency-dependent transmissions negligible.

The primary form of control in the model is achieved through the variable u, which represents a proportional reduction in the transmission coefficient, . As is standard for mass-action models of directly transmitted pathogens, reflects the product of the rate at which susceptible and infectious individuals come into contact and the probability of transmission given that a contact has occurred (Keeling and Rohani 2007). Thus, a wide range of non-pharmaceutical interventions could result in changes in u, including school closures, work from home policies, and shelter in place mandates, as well as more targeted approaches, such as isolation based on self-awareness of symptoms or contact tracing. In addition to control through u, the V class represents individuals who have been vaccinated, with individuals entering V from S at rate beginning on day . We assume that vaccination may not provide complete protection, resulting in vaccinated individuals becoming infected at a fraction of the rate at which fully susceptible individuals become infected.

Basic Reproduction Number,

At its core, the behavior of this model is similar to that of an SEIR model with demography. Because analyses of the transient and asymptotic properties of this class of models are plentiful in textbooks and elsewhere (e.g., Keeling and Rohani 2007), we omit such an analysis here. We do, however, derive a formula to express the basic reproduction number, , as a function of model parameters. This relationship plays a role in how we parameterize the model.

We use the next-generation method (van den Driessche and Watmough 2008) to obtain a formula describing as a function of model parameters. This method depends on matrices and , whose elements are defined as the rates at which secondary infections increase the ith compartment and the rates at which disease progression, death, and recovery decrease the ith compartment, respectively. For “disease compartments” E, A, I, and H, these matrices are defined asThese matrices are then used to define two others,where x represents disease compartments and y non-disease compartments. The inverse of V iswhich, along with F, defines the matrix . Specifyingwe obtainas the maximum eigenvalue of K.

Optimal Control Problem

The optimal control problem is to minimize the objective functionalsubject to the constraints of the state dynamics described in Eq. (1) and the initial conditions that and all other state variables equal zero at time . The parameter c weights the extent to which the control, u(t), is prioritized for minimization relative to deaths, D(t). Squared terms for D(t) and u(t) are chosen both for mathematical convenience in the case of u(t) and to more heavily penalize solutions for u(t) that permit relatively high values of D(t) or u(t) at any given time.

To find the optimal control, , that minimizes J(u), we follow standard results from optimal control theory applied to systems of ordinary differential equations (Lenhart and Workman 2007). These techniques make use of Pontryagin’s Maximum Principle to determine the pointwise minimum of the Hamiltonian of the system, , using adjoint variables, , that correspond to each of the state variables.

To solve for the adjoint variables, we first define the Hamiltonian,We then define differential equations describing the behavior of each adjoint variable as the negative of the partial derivative of with respect to the state variable corresponding to each adjoint variable. This yieldswhereEquation (9) can be solved backward in time with transversality conditions at time equal to zero for each adjoint variable.

To find the pointwise optimal control, , we find the value of u that minimizes , which yieldsThe optimal control is subject to a lower bound of 0 and an upper bound of , with those values used for whenever the right-hand side of Eq. (10) yields values outside those bounds.

To find numerically, we use the forward–backward sweep method (Lenhart and Workman 2007), which involves first solving for the state variables forward in time, next solving for the adjoint variables backward in time, and then plugging the solutions for the relevant state and adjoint variables into Eq. (10), subject to bounds on u(t). As is often the case for optimal control problems (Lenhart and Workman 2007), we found that we needed to perform these steps iteratively and with a convex combination of controls across iterations to achieve convergence. Specifically, we repeated the forward–backward sweep process 50 times, after which we took newly proposed solutions of as the average of the 20 most recent solutions, until the algorithm was stopped after 2000 iterations. We assessed convergence with a statistic defined aswhere smaller values indicate better convergence. All numerical solutions were obtained with a Runge–Kutta 4 routine implemented with the ode function from the deSolve package (Soetaert et al. 2010) in R.

Model Parameterization

In Table 2, we specify low, intermediate, and high values of each parameter. For some parameters, estimates matching our parameter definitions were taken from other studies. For other parameters, additional steps were necessary to match our parameter definitions. We elaborate on the latter below.

Transmission coefficient, We based values of this parameter on assumptions about by solving for as a function of and other parameters in Eq. (6). Because estimates of for COVID-19 vary widely (Park et al. 2020), we chose values that span a range of estimates that may be applicable to the USA.

Background birth and death rates, and We parameterized consistent with a rate of 3,791,712 births in a population of 331 million in the USA in 2018 (Martin et al. 2019). To achieve a constant population size in the absence of COVID-19, we set equal to .

Probability of death among hospitalized cases, and h We assume that is equal to early estimates from the USA (2.6%) (Centers for Disese Control and Prevention 2020) and that is equal to estimates from Italy (7.2%) (Onder et al. 2020). Estimates of how quickly might be approached have not been made empirically, so the value of h was assumed. The intermediate value of corresponds to an increase in H of 50% beyond resulting in 50% of the maximum increase from to .

Progression through hospitalization, We used line-list data [1] to estimate a mean (13.2 days) and standard deviation (7.39 days) of the time between hospital admission and either discharge or death.

Timing of vaccine introduction, Experts have stated that a vaccine against COVID-19 could be available to the public by spring 2021 (Amanat and Krammer 2020). We used April 1, 2021, as our default value for this parameter.

Vaccination rate, In the 2009 H1N1 pandemic, 100 million doses of vaccine were administered between October 2009 and April 2010 [8,9]. Consistent with that, our default value of resulted in 0.197% of the population being vaccinated each day.

Hospital capacity, We based our definition of this parameter on hospital beds and did not include other aspects of hospital resources, such as ICU beds, ventilators, or hospital staffing. Specifically, we adopted an estimate of 312,090 beds available for COVID-19 patients in the USA by the Institute for Health Metrics and Evaluation [54].

Maximum effect of control, A model incorporating survey data and age-based contact patterns estimated that social distancing could reduce transmission of SARS-CoV-2 by 73% (Jarvis et al. 2020). Based on this, we used 0.7 as an intermediate value of this parameter and 0.5 and 0.9 as lower and upper values.

Model Calibration

To obtain realistic behavior of the model, we calibrated it to match the cumulative number of reported deaths in the USA within the first 100 days of 2020 (i.e., by April 9), which we obtained from the New York Times [39]. We focused our calibration on the parameter , due to the difficulty of empirically estimating the rate at which imported infections appear. To obtain a value of that resulted in the model matching the reported number of deaths, we simulated the model across 300 values of evenly spaced between and on a log scale, performed a linear interpolation of the simulated number of deaths across those values of , and found the value of that most closely matched reported deaths.

We calibrated the model under a total of 18 different parameter scenarios, crossing low, intermediate, and high values of with low, intermediate, and high values of and low and high values of . The latter represents the proportion of all deaths caused by COVID-19 that were reported. Because non-pharmaceutical interventions began going into effect in the USA within the timeframe of this calibration period, we calibrated the model subject to an assumed pattern of u(t) through the first 100 days of 2020. We chose a logistic functional form for this, with a minimum of 0 and a maximum of . Two parameters that control the midpoint and slope of the increase from 0 to were selected by an informal process of trial and error, with the goal of having the model’s predictions of deaths over time match the timing of reported deaths under all 18 parameter scenarios. We also used this process to select the date on which importations were initiated through .

Results

Under each of the 18 different scenarios we considered about , , and , our model was able to reproduce the numbers and timing of reported deaths in the USA reasonably well (Fig. 2). Different values of were required to do so under different values of those other three parameters, but all values of ranged (Table 3). Taking into account the population of the USA, this equates to a range of 33–331 imported infections per day. Higher values of resulted in lower values of , given that fewer importations were required to generate the reported numbers of deaths when there was more local transmission. For similar reasons, was lower when was lower. The value of had a negligible effect on .

Fig. 2

Reported (red) and simulated (black) numbers of daily deaths in the USA resulting from model calibration under 18 different parameter scenarios (black lines) (Color Figure Online)

Table 3

Calibrated estimates of the importation rate, , under 18 scenarios about the values of , , and

R_0	u_max	ω	ι
3	0.5	0.5	9.35×10^-7
3.5	0.5	0.5	3.16×10^-7
4	0.5	0.5	1.07×10^-7
3	0.7	0.5	9.69×10^-7
3.5	0.7	0.5	3.32×10^-7
4	0.7	0.5	1.14×10^-7
3	0.9	0.5	9.35×10^-7
3.5	0.9	0.5	3.16×10^-7
4	0.9	0.5	1.07×10^-7
3	0.5	0.8	9.69×10^-7
3.5	0.5	0.8	3.32×10^-7
4	0.5	0.8	1.14×10^-7
3	0.7	0.8	9.35×10^-7
3.5	0.7	0.8	3.16×10^-7
4	0.7	0.8	1.07×10^-7
3	0.9	0.8	9.69×10^-7
3.5	0.9	0.8	3.32×10^-7
4	0.9	0.8	1.13×10^-7

In addition to , we also calibrated the date on which importations began (February 1) and u(t) during the first 100 days of 2020. Even though the first imported cases in the USA began appearing before February 1 (Holshue et al. 2020), it is not surprising that our model better matched the data with a slightly later start date, given that the model did not allow for increases in importation over time, as likely occurred (Perkins et al. 2020). The calibration of u(t) was consistent with u(t) starting at zero, reaching half of on April 4, 2020, and increasing by a maximum of per day around that time. The timing of these changes is approximately two weeks later than changes in Google mobility data from the USA [22], which may be a consequence of our model moving some individuals from exposure to death too quickly given that residence time in each compartment is exponentially distributed.

In general, we do not interpret any of the calibrated parameter values as reliable estimates of empirical quantities. Such an interpretation would require analyses that more carefully account for data-generating processes and sources of uncertainty. Rather, our objective in this calibration exercise was to ensure that the model’s behavior is reasonably consistent with empirical data, which we feel was accomplished.

Convergence

For our main results, we obtained solutions to under a total of 126 parameter combinations, crossing the 18 parameter sets in Table 3 with seven values of c spanning to by factors of ten. In one representative example (Fig. 3, left), we observed that solutions of converged to a set of similar solutions after several hundred iterations. As expected, the objective functional decreased during those initial iterations and remained low thereafter (Fig. 3, right). To assess convergence by the statistic in Eq. (11), we selected the iterations with the ten lowest values of J(u) in the last 100 iterations. Of 126 parameter combinations, 78% had final solutions of with convergence statistics below , 94% below , and all below .

Fig. 3

Convergence of solutions of under parameters , , , and . Left: Colors indicate values of for each day in 2020 and 2021 across 2,000 iterations of the forward–backward sweep algorithm. Right: Across iterations, the value of the objective functional, J(u), decreased steadily until cycling for the remaining iterations (Color Figure Online)

Minimized Objective Functional Under Different Parameters

Comparing values of the two components of the objective functional, and , across different parameter values provides insight into the range of behavior of the model and its response to control. Because they are similar to the components of the objective functional but more easily interpretable, we describe effects of model parameters on (cumulative deaths) and (cumulative time under control). The parameter c, which controls the weight of the two terms in J(u), led to a transition from minimizing cumulative deaths to minimizing cumulative time spent under control as it varied from to (Fig. 4). The parameters and both influenced the extent to which cumulative deaths could be minimized. With a high value of , deaths could be kept relatively low across all values of explored, provided that . With a low value of and a high value of , cumulative deaths could only be reduced by around 10% as c varied across its entire range from to . The parameter had essentially no influence on the components of the objective functional (Fig. 4).

Fig. 4

Dependence of time under control (blue) and cumulative deaths (red) on c (x-axis), (columns), (rows), and (markers). Deaths are quantified as a proportion of the overall population (Color Figure Online)

Optimal Control Over Time

We first consider the scenario where and , because those are the conditions under which control has the greatest potential to influence the pandemic. Under , remains at until late June 2020, after which it settles down to around 75% of until late 2021 (Fig. 5). This results in hospitalizations dropping from their peak in April 2020 and remaining very low through 2021. The susceptible population remains very high and only begins eroding once a vaccine is introduced.

Fig. 5

Optimal control under parameters with maximal ability to control the pandemic, with maximal weighting on minimization of deaths. Panels show the optimal control (bottom) and its impacts on the dynamics of hospitalized (middle) and susceptible (top) compartments with , , , and . Vaccination is introduced at the time indicated by the arrow, with the vaccinated population (orange) reducing the susceptible population. Through April 2020 (gray shading), the value of the control is fixed according to its calibrated trajectory. The dashed horizontal line indicates hospital capacity (Color Figure Online)

With a higher value of , drops to around 50% of in May 2020 (Fig. 5). As a result, hospitalizations rebound and exceed hospital capacity by around a third in June and July before falling again, after returns to . From August 2020 onward, a steady decline in allows hospitalizations to be maintained at moderate levels before rebounding and exceeding hospital capacity once again for several months in 2021.

With the highest value of , drops almost to zero at the beginning of May 2020 (Fig. 7). This results in a rapid increase in hospitalizations, which is followed by an increase in . By the end of June, the susceptible population has been depleted to the point that herd immunity begins to obviate the need for control and declines to zero by October 2020. During this large second wave in summer 2020, hospital capacity is exceeded by more than 20-fold. This results in cumulative deaths equaling 5% of the population, which is approximately one order of magnitude greater than when and two orders of magnitude greater than when .

Fig. 7

Optimal control under parameters with maximal ability to control the pandemic, with minimal weighting on minimization of deaths. Panels show the optimal control (bottom) and its impacts on the dynamics of hospitalized (middle) and susceptible (top) compartments with , , , and . Vaccination is introduced at the time indicated by the arrow, with the vaccinated population (orange) reducing the susceptible population. Through April 2020 (gray shading), the value of the control is fixed according to its calibrated trajectory. The dashed horizontal line indicates hospital capacity (Color Figure Online)

Under other conditions about the transmissibility of the virus and the potential for non-pharmaceutical interventions to reduce transmission, control is less capable of curtailing the epidemic. With and , hospitalizations peak in early summer 2020 at twice hospital capacity and take until December 2020 to fall below hospital capacity, even with the most aggressive value of c that we considered () (Fig. 8). With , , and , hospitalizations peak in summer 2020 at 30 times hospital capacity (Fig. 9). This results in extensive herd immunity and the relaxation of control in 2021, albeit at the cost of extensive deaths (Fig. 4). Optimal controls under a wider variety of parameter combinations can be explored interactively at http://covid19optimalcontrol.crc.nd.edu.

Fig. 8

Optimal control under parameters with maximal ability to control the pandemic, with maximal weighting on minimization of deaths. Panels show the optimal control (bottom) and its impacts on the dynamics of hospitalized (middle) and susceptible (top) compartments with , , , and . Vaccination is introduced at the time indicated by the arrow, with the vaccinated population (orange) reducing the susceptible population. Through April 2020 (gray shading), the value of the control is fixed according to its calibrated trajectory. The dashed horizontal line indicates hospital capacity (Color Figure Online)

Fig. 9

Optimal control under parameters with maximal ability to control the pandemic, with maximal weighting on minimization of deaths. Panels show the optimal control (bottom) and its impacts on the dynamics of hospitalized (middle) and susceptible (top) compartments with , , , and . Vaccination is introduced at the time indicated by the arrow, with the vaccinated population (orange) reducing the susceptible population. Through April 2020 (gray shading), the value of the control is fixed according to its calibrated trajectory. The dashed horizontal line indicates hospital capacity (Color Figure Online)

Across the scenarios that we considered, control generally either relaxes slowly once vaccination commences (Figs. 5, 6) or abruptly once herd immunity is attained (Fig. 7). Near the end of some scenarios, control drops markedly as the end of 2021 is approached (e.g., Fig. 8), due to the fact that deaths resulting from relaxation of control near the end of 2021 would not manifest appreciably until 2022. Such behavior should be interpreted as an artifact of the finite time horizon of our analysis. Were we to extend the analysis until the end of 2022, for example, we would expect behavior similar to that in mid-2021 to extend forward through much of 2022, and for the artifactual relaxation of control, we observe toward the end of 2021 to apply to the end of 2022 instead.

Fig. 6

Optimal control under parameters with maximal ability to control the pandemic, with intermediate weighting on minimization of deaths. Panels show the optimal control (bottom) and its impacts on the dynamics of hospitalized (middle) and susceptible (top) compartments with , , , and . Vaccination is introduced at the time indicated by the arrow, with the vaccinated population (orange) reducing the susceptible population. Through April 2020 (gray shading), the value of the control is fixed according to its calibrated trajectory. The dashed horizontal line indicates hospital capacity (Color Figure Online)

Optimal Control Following Different Starting Conditions

Our calibration procedure resulted in a single assumption about the trajectory of u(t) prior to April 30, 2020. Control during this period was fixed in analyses in Sect. 3.4, with the flexibility to define only allowed during the period from May 1, 2020, to December 31, 2021. Here, we explore how non-pharmaceutical interventions initiated one to three weeks earlier or later in spring 2020 would affect optimal controls in the period after. We focus this analysis on scenarios in which , which results in the optimal control problem seeking to minimize deaths as aggressively as any scenario that we explored. Consequently, we interpret changes in the optimal control in this section to reflect changes in constraints on what solutions of are possible, rather than changes in the balance between and in the optimization.

Across all combinations of , , and that we considered, cumulative deaths through 2020 and 2021 decrease when control begins earlier and increase when control begins later (Fig. 10, red). In the scenario in which a delay in the initiation of control has the smallest effect (, ), cumulative deaths increase by 10% with a three-week delay. In the scenario in which a delay in the initiation of control has the largest effect (, ), cumulative deaths increase 28-fold with a three-week delay.

Fig. 10

Dependence of time under control (blue) and cumulative deaths (red) on a shift in the timing of u(t) before April 30, 2020 (x-axis), (columns), (rows), and (markers). Deaths are quantified as a proportion of the overall population (Color Figure Online)

The overall amount of time spent under control throughout 2020 and 2021 increases when control begins earlier and decreases when control begins later (Fig. 10, blue). This is the case across all combinations of , , and that we considered. In part, this owes to less time spent under control through April 30, 2020, when u(t) is fixed and not subject to optimization. At the same time, delays in the initiation of control result in a higher prevalence of infection by the beginning of the optimization period, which results in higher levels of subsequent transmission, greater depletion of the susceptible population, and less need for control later in the period of optimization (compare Fig. 11 with Fig. 5).

Fig. 11

Optimal control (bottom) and its impacts on the dynamics of hospitalized (middle) and susceptible (top) compartments with , , , , and the initiation of control delayed by 21 days. Vaccination is introduced at the time indicated by the arrow, with the vaccinated population (orange) reducing the susceptible population. Through April 2020 (gray shading), the value of the control is fixed according to its calibrated trajectory. The dashed horizontal line indicates hospital capacity (Color Figure Online)

Discussion

Under our model, the ability of non-pharmaceutical interventions to minimize deaths depends to a large extent on the maximum effect that they could have on transmission, as captured by the parameter . If is high (0.9), our model predicts that maintaining transmission at a low enough level that hospital capacity will not be exceeded could be possible. If is low (0.5), our model predicts that there may be limited scope to curtail the pandemic, even if control is sustained for a long period of time. If is intermediate (0.7), our model predicts that the potential impact of control is sensitive to the value of . These results emphasize the importance of careful estimation of these parameters as the pandemic progresses (Flaxman et al. 2020). This includes accounting for geographic differences in both and (Gilbert et al. 2020; Hilton and Keeling 2020).

The balance between minimizing deaths versus days under control is determined by the parameter c in our model, which mirrors the weighting between these factors that government leaders will need to consider as they make decisions in coming months. In the event that non-pharmaceutical interventions are effective, our analysis shows that they would need to be sustained at a high level until at least sometime in summer 2020, at which time they could potentially be relaxed to a small degree but would still need to be maintained at a relatively high level thereafter. Scenarios that place greater weight on minimizing days under control provide insight into the possible consequences of relaxing non-pharmaceutical interventions prematurely. In a scenario in which the effect of non-pharmaceutical interventions is reduced moderately in May, resuming non-pharmaceutical interventions at high levels soon after becomes necessary to react to a second wave later in the summer. In a scenario in which non-pharmaceutical interventions are reduced more drastically in May, an extremely large second wave occurs in summer 2020 that exceeds hospital capacity many times over.

Our conclusion that prolonged control is needed to avoid a resurgence that would greatly exceed healthcare capacity is in line with results from other modeling studies (Davies et al. 2020; Ngoonghala et al. 2020; Ferguson et al. 2020; Kissler et al. 2020; Tuite et al. 2020). Some differ though in that they consider scenarios in which control is implemented intermittently, with periods of relaxed control in between periods of maximal control (Ferguson et al. 2020; Kissler et al. 2020; Tuite et al. 2020). In the limited number of direct comparisons of intermittent and continuous strategies that have been performed for COVID-19 to date, continuous strategies appear to be capable of more effectively limiting transmission to low levels (Djidjou-Demasse 2020; Yap and Raja 2020). The sustainability of either strategy would benefit from the ability to transition away from heightened social distancing and more toward diagnostic testing, contact tracing, and case isolation (Tuite et al. 2020; Piguillem and Shi 2020). Although our analysis focuses on a single timing of vaccine introduction in April 2021, we assume that similar results would apply until a vaccine is available, whenever that may be.

Our results tend to agree with those from other optimal control analyses of COVID-19, although our analysis goes beyond those studies in some ways. In general, ours and other studies share the conclusion that heightened control early in the pandemic is important for achieving long-term success. Like our analysis, some anticipate that partial relaxation of controls may be possible over time (Djidjou-Demasse 2020; Shah et al. 2020), whereas others focus on strategies intended to have a more limited duration to begin with (Patterson-Lomba 2020; Morris et al. 2020; Piguillem and Shi 2020). Importantly though, our analysis goes further than others in exploring sensitivity of the optimal control to model parameters. Specifically, we show that preventing a large wave that overwhelms health systems may not even be possible under some parameter combinations (low , high ) and that prioritizing the minimization of deaths versus days under control leads to vastly different outcomes. We also constrain levels of control applied through April 2020, making the optimization more relevant to decision making thereafter. Shifting the timing of the initiation of control shows that constraints about past levels of control strongly affect future possibilities for the extent to which deaths can be minimized and the level of control required to achieve that, consistent with other results (Lai et al. 2020).

The goal of our analysis is to provide qualitative insights into the implications of alternative approaches to control, rather than to make quantitative predictions about future events. Some key limitations of using our model for the latter purpose include our omission of subnational variation in epidemic dynamics (Perkins et al. 2019; Team 2020), differentiation among alternative non-pharmaceutical interventions (Flaxman et al. 2020), and age differences in contact patterns (Jarvis et al. 2020), susceptibility (Davies et al. 2020), and risk of hospitalization (Centers for Disese Control and Prevention 2020). Whereas a previous optimal control analysis of pandemic influenza (Shim 2013) suggested that age-specific optimal controls were all relatively similar, recent work on COVID-19 (Richard et al. 2020; Gondim and Machado 2020) suggests that optimal controls should be higher for older age-groups due to their higher risk of severe disease and death. Inclusion of age structure is important for other reasons too, such as realistically capturing transmission dynamics (Britton et al. 2020) and accounting for age-specific interventions, such as school closures (Head et al. 2020). Additional limitations that affect our model’s suitability for making future predictions include its deterministic nature and the rudimentary calibration procedure that we performed, which was sufficient to provide a basis for qualitative analyses but that would need refinement for application of our model to inference or forecasting.

In conclusion, our analysis suggests that May 2020 was a critical juncture in the pandemic, when decisions about the continuation or relaxation of non-pharmaceutical interventions had major implications for the possibility of keeping transmission below levels that health systems could cope with. At the same time, our analysis highlights the role that constraints play in determining optimal levels of control, both in terms of constraints on epidemiological parameters and on levels of control prior to the time that a decision is made about future actions. At any point during the pandemic, reducing transmission in the near term would give decision makers greater flexibility in the range of decisions available to them in the long term, and gathering high-quality data could help reduce uncertainty about the consequences of those decisions.