key: cord-340851-lnv8yi9l
authors: Perkins, Alex; Espana, Guido
title: Optimal control of the COVID-19 pandemic with non-pharmaceutical interventions
date: 2020-04-23
journal: nan
DOI: 10.1101/2020.04.22.20076018
sha: 
doc_id: 340851
cord_uid: lnv8yi9l

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

The first known case of COVID-19 in the United States arrived from China on January 15, 2020 [22] . In the weeks that followed, deficiencies in testing allowed the virus that causes this disease to spread largely undetected [41] , with the number of reported cases growing to 429,319 by April 8, 2020 [36] . Social distancing-such as closing schools, working from home, and sheltering in place-has been adopted widely and appears capable of impacting transmission [26, 12] . These practices have tremendous economic and social consequences, however, which has precipitated growing desire to relax them [16] .

A danger in relaxing the use of social distancing is that the virus could resurge in its absence [44] . 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 [30] and one is not likely to be available until at least spring 2021 [4] . Strategies for successfully controlling COVID-19 until then will depend on a suite of non-pharmaceutical interventions (NPIs) [3] , including some degree of social distancing but also diagnostic testing, contact tracing, and case isolation [5] .

In models of pathogen transmission based on mass-action assumptions, NPIs 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 [24] . 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 [18] . Until alternative NPI-based strategies 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 [28] . These techniques have been widely applied to a variety of pathogen transmission systems before (e.g., [6, 33, 11, 2] ), including in the context of a pandemic of a respiratory pathogen (e.g., [29, 48, 46] ). The latter studies indicate that the level of NPI-based control required is dependent on model parameters, and that application of NPIs 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 NPIs to control COVID-19, with a focus on the US. 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 [14, 31, 40, 43, 45] , although these have been less focused on any particular geographic setting. One contextual feature of the US 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 going forward, conditional on what has happened already. In addition, we emphasize the sensitivity of our results to epidemiological parameters that are not yet well characterized but appear influential in determining what range of impacts of NPIs on COVID-19 are possible. All code and data used in our analysis is available at http://github.com/TAlexPerkins/covid19optimalcontrol.

We modeled SARS-CoV-2 transmission according to a system of ordinary differential equations,

with variables and parameters defined in Tables 1 and 2. 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 ρ −1 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 γ −1 days. All asymptomatic infections and a proportion 

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 Hmax represents the hospital capacity above which the probability of death exceeds ∆ − . Hospitalizations are quantified as a proportion of the overall population.

1 − κ 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 η −1 days.

Rather than track deaths due to COVID-19, D(t), as a state variable in eqn.

(1), we assume that they follow directly from hospitalizations, H(t), according to

This results in the probability of death moving beyond a minimum of ∆ − towards a maximum of ∆ + as H(t) exceeds H max , as illustrated in Figure 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 [43] , albeit with a different functional form. Individuals who are recovered and immune are not followed explicitly, as the model's assumption of density-dependent transmission only requires spec-. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04. 22.20076018 doi: medRxiv preprint ification 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 frequencydependent transmission 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 [24] . Thus, a wide range of NPIs 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 1 − of the rate at which fully susceptible individuals become infected.

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., [24] ), we omit such an analysis here. We do, however, derive a formula to express the basic reproduction number, R 0 , as a function of model parameters. This relationship plays a role in how we parameterize the model.

We use the next-generation method [15] to obtain a formula describing R 0 as a function of model parameters. This method depends on matrices F and V, whose elements are defined as the rates at which secondary infections increase the i th compartment and the rates at which disease progression, death, and recovery decrease the i th compartment, respectively. For "disease compartments" E, A, I, and H, these matrices are defined as

These matrices are then used to define two others, . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted April 23, 2020. 

where x represents disease compartments and y non-disease compartments. The inverse of V is

which, along with F , defines the matrix K = F V −1 . Specifying

we obtain

as the maximum eigenvalue of K.

The optimal control problem is to minimize the objective functional

subject to the constraints of the state dynamics described in eqn.

(1) and the initial conditions that S(t 0 ) = 1 and all other state variables equal zero at time t 0 . 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, u * (t), that minimizes J(u), we follow standard results from optimal control theory applied to systems of ordinary differential . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04. 22.20076018 doi: medRxiv preprint equations [28] . These techniques make use of Pontryagin's Maximum Principle to determine the pointwise minimum of the Hamiltonian of the system, H, 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 H with respect to the state variable corresponding to each adjoint variable. This yields

where

Eqn. (9) can be solved backwards in time with transversality conditions at time t 1 equal to zero for each adjoint variable. To find the pointwise optimal control, u * , we find the value of u that minimizes ∂H ∂u , which yields . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. 

The optimal control is subject to a lower bound of 0 and an upper bound of u max , with those values used for u * whenever the right hand side of eq. (10) yields values outside those bounds. To find u * (t) numerically, we use the forward-backward sweep method [28] , 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 eqn. (10) , subject to bounds on u(t). As is often the case for optimal control problems [28] , 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 u * (t) as the average of the 20 most recent solutions, until the algorithm was stopped after 2,000 iterations. We assessed convergence with a statistic defined as

where 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 [47] in R.

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 R 0 by solving for β as a function of R 0 and other parameters in eqn. (6) . Because estimates of R 0 for COVID-19 vary widely [39] , we chose values that span a range of estimates that may be applicable to the US.

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 US in 2018 [32] . To achieve a constant population size in the absence of COVID-19, we set δ equal to µ.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04. 22.20076018 doi: medRxiv preprint Probability of death among hospitalized cases, ∆ and h. We assume that ∆ − is equal to early estimates from the US (2.6%) [9] and that ∆ + is equal to estimates from Italy (7.2%) [38] . Estimates of how quickly ∆ + might be approached have not been made empirically, so the value of h was assumed. The intermediate value of h = 701 corresponds to an increase of H of 50% beyond H max 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 [4] . 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 [7, 8] . Consistent with that, our default value of ν resulted in 0.197% of the population being vaccinated each day.

Hospital capacity, H max . 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 US by the Institute for Health Metrics and Evaluation [50] .

Maximum effect of control, u max . A model incorporating survey data and agebased contact patterns estimated that social distancing could reduce transmission of SARS-CoV-2 by 73% [23] . 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.

To obtain realistic behavior of the model, we calibrated it to match the cumulative number of reported deaths in the US within the first 100 days of 2020 (i.e., by April 9), which we obtained from the New York Times [36] . 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 10 −12 and 10 −4 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 R 0 with low, intermediate, and

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. high values of u max and low and high values of ω. The latter represents the proportion of all deaths caused by COVID-19 that were reported. Because NPIs began going into effect in the US 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 u max . Two parameters that control the midpoint and slope of the increase from 0 to u max 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 ι.

Under each of the 18 different scenarios we considered about R 0 , u max , and ω, our model was able to reproduce the numbers and timing of reported deaths in the US reasonably well (Figure 2) . Different values of ι were required to do so under different values of those other three parameters, but all values of ι ranged 10 −7 -10 −6 ( Table 3) . Taking into account the population of the US, this equates to a range of 33 -331 imported infections per day. Higher values of R 0 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 u max was lower. The value of ω had a negligible effect on ι.

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 US began appearing before February 1 [22] , it is . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04. 22.20076018 doi: medRxiv preprint 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 [41] . The calibration of u(t) was consistent with u(t) starting at zero, reaching half of u max on April 4, 2020, and increasing by a maximum of 0.075 u max per day around that time. The timing of these changes is approximately two weeks later than changes in Google mobility data from the US [20], 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. 

For our main results, we obtained solutions to u * (t) under a total of 126 parameter combinations, crossing the 18 parameter sets in Table 3 with seven values of c spanning 10 −12 to 10 −6 by factors of ten. In one representative example (Figure 3 , left), we observed that solutions of u * (t) 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 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. ( Figure 3 , right). To assess convergence by the statistic in eqn. (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 u * (t) with convergence statistics below 10 −5 , 94% below 10 −2 , and all below 10 −1 .

Comparing values of the two components of the objective functional, t1 t0 D(t) 2 dt and t1 t0 u(t) 2 dt, 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 t1 t0 D(t)dt (cumulative deaths) and t1 t0 u(t)dt (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 10 −9 to 10 −6 ( Figure 4 ). The parameters R 0 and u max both influenced the extent to which cumulative deaths could be minimized. With a high value of u max , deaths could be kept relatively low across all values of R 0 explored, provided that c ≤ 10 −10 . With a low value of u max and a high value of R 0 , cumulative deaths could only be reduced by around 10% as c varied across its entire range from 10 −12 to 10 −6 . The parameter ω had essentially no influence on the components of the objective functional ( Figure 4 ).

We first consider the scenario where R 0 = 3 and u max = 0.9, because those are the conditions under which control has the greatest potential to influence . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. the pandemic. Under c = 10 −12 , u * (t) remains at u max until late June 2020, after which it settles down to around 75% of u max until late 2021 ( Figure 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. With a higher value of c = 10 −9 , u * (t) drops to around 50% of u max in May 2020 ( Figure 5 ). As a result, hospitalizations rebound and exceed hospital capacity by around a third in June and July before falling again, after u * (t) returns to u max . From August 2020 onward, a steady decline in u * (t) 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 c = 10 −6 , u * (t) drops almost to zero at the beginning of May 2020 ( Figure. 7) . This results in a rapid increase in hospitalizations, which is followed by an increase in u * (t). 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 u * (t) declines to zero by October 2020. During this large second wave in summer 2020, hospital capacity is ex-. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04.22.20076018 doi: medRxiv preprint 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 R 0 = 3, umax = 0.9, ω = 0.8, and c = 10 −9 . 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.

ceeded 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 c = 10 −9 and two orders of magnitude greater than when c = 10 −12 .

Under other conditions about the transmissibility of the virus and the potential for NPIs to reduce transmission, control is less capable of curtailing the epidemic. With R 0 = 3.5 and u max = 0.7, 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 (10 −12 ) (Figure 8 ). With R 0 = 4, u max = 0.5, and c = 10 −12 , hospitalizations peak in summer 2020 at thirty times hospital capacity (Figure 9 ). This results in extensive herd immunity and the relaxation of control in 2021, albeit at the cost of extensive deaths (Figure 4) . Optimal controls under a wider variety of parameter combinations can be explored interactively at http://covid19optimalcontrol.crc.nd.edu.

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 Section 3.4, with the flexibility to define u * (t) only allowed during the period from May 1, 2020 to December 31, 2021. Here, we explore how NPIs . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04.22.20076018 doi: medRxiv preprint 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 R 0 = 3.5, umax = 0.7, ω = 0.8, and c = 10 −12 . 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. 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 R 0 = 4, umax = 0.5, ω = 0.8, and c = 10 −12 . 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.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04. 22.20076018 doi: medRxiv preprint 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 c = 10 −12 , 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 u * (t) are possible, rather than changes in the balance between D(t) 2 and u(t) 2 in the optimization.

Across all combinations of R 0 , u max , and ω that we considered, cumulative deaths through 2020 and 2021 decrease when control begins earlier and increase when control begins later (Figure 10, red) . In the scenario in which a delay in the initiation of control has the smallest effect (R 0 = 3, u max = 0.5), 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 (R 0 = 4, u max = 0.9), cumulative deaths increase 28-fold with a three-week delay.

The overall amount of time spent under control throughout 2020 and 2021 increases when control begins earlier and decreases when control begins later ( Figure 10 , blue). This is the case across all combinations of R 0 , u max , 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. In addition though, 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 Figure 11 with Figure 5 ).

Under our model, the ability of NPIs to minimize deaths depends to a large extent on the maximum effect that NPIs could have on transmission, as captured by the parameter u max . If u max 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 u max 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 u max is intermediate (0.7), our model predicts that the potential impact of control is sensitive to the value of R 0 . These results emphasize the importance of careful estimation of these parameters as the pandemic progresses [18] . This includes accounting for geographic differences in both u max and R 0 [19, 21] .

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 NPIs are effective, our analysis shows that they will need to be sustained at a high level until sometime in summer 2020, at which time they can be relaxed to a small degree but will still need to be maintained at a relatively high level thereafter. Scenarios that . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. place greater weight on minimizing days under control provide insight into the possible consequences of relaxing NPIs prematurely. In a scenario in which the effect of NPIs is reduced moderately in May, resuming NPIs at high levels soon after becomes necessary to react to a second wave later in the summer. In a scenario in which NPIs 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 [13, 37, 17, 25, 51] . 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 [17, 25, 51] . 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 [14, 52] . The sustainability of either strategy would benefit from the ability to transition . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted April 23, 2020. away from heightened social distancing and more toward diagnostic testing, contact tracing, and case isolation [51, 43] .

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 [14, 45] , whereas others focus on strategies intended to have a more limited duration to begin with [40, 35, 43] . 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 u max , high R 0 ) 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 future decision making. 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 [27] .

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 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04. 22.20076018 doi: medRxiv preprint the latter purpose include our omission of subnational variation in epidemic dynamics [42, 49] , differentiation among alternative NPIs [18] , and age differences in contact patterns [23] and risk of hospitalization [9] . At least with respect to the possibility of different levels of control for different age groups, one study found that age-specific optimal controls were all relatively similar for NPIs against pandemic influenza [46] . 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 we are at a critical juncture in the pandemic, when decisions about the continuation or relaxation of NPI-based control strategies could have major implications for the possibility of keeping transmission below levels that health systems can cope with. At the same time, our analysis highlights the role that constraints play in determining optimal levels of control going forward, 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. Going forward, reducing transmission in the near term would give decision makers greater flexibility in the range of decisions available to them in the future, and gathering high-quality data could help reduce uncertainty about the consequences of those decisions.

. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04. 22.20076018 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted April 23, 2020. . https://doi.org/10.1101/2020.04. 22.20076018 doi: medRxiv preprint 

nCoV 2019 Data Working Group: Epidemiological Data from the nCoV-2019 Outbreak: Early Descriptions from Publicly Available Data

Optimal control strategies for dengue transmission in pakistan

Non-pharmaceutical public health interventions for pandemic influenza: an evaluation of the evidence base

CNN: Lead scientist says coronavirus vaccine could be ready soon

Covid-19: towards controlling of a pandemic

Backward bifurcation and optimal control in transmission dynamics of west nile virus

H1N1 Flu Vaccine

H1N1 Pandemic Timeline

Control and Prevention: Severe outcomes among patients with coronavirus disease 2019 (covid-19) -united states

The effect of travel restrictions on the spread of the 2019 novel coronavirus (covid-19) outbreak

Optimal tuberculosis prevention and control strategy from a mathematical model based on real data

Impact assessment of non-pharmaceutical interventions against coronavirus disease 2019 and influenza in hong kong: an observational study

The effect of nonpharmaceutical interventions on COVID-19 cases, deaths and demand for hospital services in the UK: a modeling study

Optimal covid-19 epidemic control until vaccine deployment. medRxiv

Further notes on the basic reproduction number

Editorial: Covid-19 in the usa: a question of time

Report 9 -Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand

Report 13 -Estimating the number of infections and the impact of nonpharmaceutical interventions on COVID-19

Preparedness and vulnerability of african countries against importations of covid-19: a modelling study

Estimation of country-level basic reproductive ratios for novel coronavirus (covid-19) using synthetic contact matrices. medRxiv

First case of 2019 novel coronavirus in the united states

Quantifying the impact of physical distance measures on the transmission of covid-19 in the uk

Modeling Infectious Diseases in Humans and Animals

Projecting the transmission dynamics of sars-cov-2 through the postpandemic period

The effect of human mobility and control measures on the covid-19 epidemic in china

Effect of non-pharmaceutical interventions for containing the covid-19 outbreak in china

Optimal Control Applied to Biological Models

An optimal control theory approach to nonpharmaceutical interventions

Developing covid-19 vaccines at pandemic speed

Optimal control applied to a seir model of 2019-ncov with social distancing. medRxiv

Births: Final data for 2018

Modeling optimal intervention strategies for cholera

Estimating the asymptomatic proportion of coronavirus disease 2019 (covid-19) cases on board the diamond princess cruise ship

Optimal, near-optimal, and robust epidemic control. arXiv

New York Times: Coronavirus (Covid-19) Data in the United States

Mathematical assessment of the impact of non-pharmaceutical interventions on curtailing the 2019 novel coronavirus

Case-Fatality Rate and Characteristics of Patients Dying in Relation to COVID-19 in Italy

Reconciling early-outbreak estimates of the basic reproductive number and its uncertainty: framework and applications to the novel coronavirus (sars-cov-2) outbreak. medRxiv

Optimal timing for social distancing during an epidemic. medRxiv (2020)

Estimating unobserved sars-cov-2 infections in the united states

Heterogeneous local dynamics revealed by classification analysis of spatially disaggregated time series data

Optimal covid-19 quarantine and testing policies. EIEF Working Papers Series

The effect of control strategies to reduce social mixing on outcomes of the covid-19 epidemic in wuhan, china: a modelling study

Control strategies to curtail transmission of covid-19. medRxiv

Optimal strategies of social distancing and vaccination against seasonal influenza

Solving differential equations in r: Package desolve

Optimal control and sensitivity analysis of an influenza model with treatment and vaccination

Forecasting covid-19 impact on hospital bed-days, icu-days, ventilator-days and deaths by us state in the next 4 months

Mathematical modelling of covid-19 transmission and mitigation strategies in the population of ontario, canada

Time-variant strategies for optimizing the performance of nonpharmaceutical interventions (npis) in protecting lives and livelihoods during the covid-19 pandemic

Acknowledgements Thanks to the University of Notre Dame's Center for Research Computing for computing resources.