A simple criterion to design optimal non-pharmaceutical interventions for epidemic outbreaks ============================================================================================ * Marco Tulio Angulo * Fernando Castaños * Rodrigo Moreno-Morton * Jorge X. Velasco-Hernández * Jaime A. Moreno ## Abstract For mitigating the COVID-19 pandemic, much emphasis is made on implementing non-pharmaceutical interventions to keep the reproduction number below one. However, using that objective ignores that some of these interventions, like bans of public events or lockdowns, must be transitory and as short as possible because of their significative economic and societal costs. Here we derive a simple and mathematically rigorous criterion for designing optimal transitory non-pharmaceutical interventions for mitigating epidemic outbreaks. We find that reducing the reproduction number below one is sufficient but not necessary. Instead, our criterion prescribes the required reduction in the reproduction number according to the desired maximum of disease prevalence and the maximum decrease of disease transmission that the interventions can achieve. We study the implications of our theoretical results for designing non-pharmaceutical interventions in 16 cities and regions during the COVID-19 pandemic. In particular, we estimate the minimal reduction of each region’s contact rate necessary to control the epidemic optimally. Our results contribute to establishing a rigorous methodology to guide the design of optimal non-pharmaceutical intervention policies. ## Introduction Since the seminal work of May and Anderson [1], the design of interventions to *eradicate* infectious diseases has the objective of achieving a basic (*R*) or effective reproduction number below one [2, 3]. The underlying assumption here is that it is possible to maintain interventions for long periods, such as long-term vaccination programs. During the COVID-19 pandemic, this same objective is guiding the design of non-pharmaceutical interventions (NPIs) [4]. However, maintaining NPIs like bans of public events or lockdowns for long periods of time is infeasible because of their substantial economic and societal costs [5, 6]. Actually, instead of aiming for eradication, NPIs aim to *mitigate* the economic and social costs of the epidemic outbreak [7]. Nevertheless, we still lack simple guidelines to design NPIs for mitigating epidemic outbreaks, analogous to the *R* < *1* condition for eradication. Here we use the classical Susceptible-Infected-Removed epidemiological model to fully characterize the design of NPIs for mitigating epidemic outbreaks. With this aim, we consider that NPIs should achieve an optimal tradeo? between two objectives [8]. First, optimal NPIs must minimize the period in which they need to be applied, consequently minimizing their associated economic and societal costs. Second, optimal NPIs must guarantee that the disease prevalence does not exceed a specified maximum level, which for example can represent the health services’ capacity for that particular disease [9]. We obtain a full analytical characterization of such optimal NPIs, specifying the optimal intervention at each state that the epidemic can be. This characterization yields the necessary and sufficient criterion for the existence of optimal NPIs for mitigation, analogous to the *R* < *1* condition for eradication. We find that reducing the reproduction number below one is sufficient but not necessary for their existence. Instead, we show that the desired maximum disease prevalence determines the necessary reduction in the reproduction number. The consequence of not reducing the reproduction number below one is that interventions must start before the disease prevalence reaches the specified maximum level. We also demonstrate numerically that the derived optimal NPIs are robust to uncertainties in the model parameters and unmodeled epidemic dynamics (e.g., undetected infections). Finally, we explore the implications of our theoretical result by analyzing the response of 16 cities and regions across the globe to the COVID-19 pandemic, finding that most regions achieved a larger-than-necessary reduction in transmission. Our results contribute to designing non-pharmaceutical interventions to respond optimally and robustly against epidemic outbreaks. ## Characterizing optimal non-pharmaceutical interventions ### Optimal epidemic mitigation using NPIs Our objective is to characterize the reduction in the disease transmission that is optimal for each *state* in which the epidemic outbreak can be. For this, we leverage on the mathematical tractability of the Susceptible-Infected-Removed (SIR) model [10], where the state can be characterized by the pair (*S, I*) ∈ [0, 1]2. Here, *S* is the proportion of the population that is susceptible to the disease, and *I* is the disease prevalence (i.e., proportion of the population that is infected), see Fig. 1a. We discuss later other more detailed epidemic models. The epidemic state changes with time *t* as the disease is transmitted, producing the trajectory (*S*(*t*), *I*(*t*)) for *t* ≥ 0. For epidemic *mitigation*, we consider that the goal is keeping the disease prevalence below a specified level *I*max 2 (0, 1]. This constant may characterize, for example, the health services’ capacity in the sense that a prevalence above *I*max causes higher mortality due to hospital saturation [11]. In general, *I*max should consider all social and economic conditions of the specific population where the outbreak occurs. To keep *I*(*t*) ≤ *I*max, we assume we can apply one or several NPIs that reduce disease transmission by the factor (1 − *u*), for some *u* ∈ [0, 1], see Fig. 1a. The NPIs achieve no reduction when *u* = 0, and they completely stop transmission when *u* = 1. Since it is unfeasible to stop transmission fully, we upper-bound the reduction by *u*max ∈ (0, 1). We say that *u* is *admissible* if *u* ∈ [0, *u*max]. ![Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F1.medium.gif) [Figure 1:](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F1) Figure 1: Optimal non-pharmaceutical interventions. **a**. Susceptible-Infected-Removed (SIR) model with non-pharmaceutical interventions (NPIs) reducing disease transmission. For the optimal NPI design problem, the objective is to design the intervention *u*∗(*t*) with minimal effective duration such that *u*∗(*t*) ∈ [0, *u*max] and *I*(t) ≤ *I*max for all *t* ≥ 0. **b and c**. Panels show the response of the SIR model for two interventions (parameters are *β* = 0.52, *γ* = 1/7, *I* = 8.855 x 10−7 and *S* = 1 − *I*). Both intervention 1 and 2 satisfy *u*(*t*) ≤ *u*max and guarantee that *I*(*t*) ≤ *I*max Actually, intervention 2 is the optimal one derived using our analysis: it is intervention with minimal effective duration satisfying *I*(*t*) ≤ *I*max. **d**. The effective duration of an intervention measures the interval between the start of the outbreak and the last time that a non-zero intervention is applied. In this example, the effective duration of intervention 1 is 120 days, while the effective duration of intervention ∈ is 69 days. Different admissible NPIs can keep the disease prevalence below *I*max. For instance, “intervention 1” in the example of Fig. 1b-c keeps this restriction and has an “effective duration” of 120 days. Here, the *effective duration* of an intervention is the interval between the start of the outbreak and the last time that a non-zero intervention is applied (Fig. 1d). “Intervention 2” of Fig. 1b-c also keeps the restriction *I*(*t*) ≤ *I*max, but its effective duration is only 69 days. To design the *optimal* NPI, we ask for the intervention with minimal effective duration. Specifically, we ask for the admissible reduction *u*∗ (*S*(*t*), *I*(*t*)) required *now* (i.e., at the current state) such that: (1) it minimizes the effective duration of the intervention; and (2) it ensures that the prevalence can be maintained below *I*max for *all future* time by using some admissible intervention. If the optimal NPI problem has a solution *u*∗, then *u*∗(*S, I*) characterizes the optimal reduction in the disease transmission that the NPIs should achieve if the epidemic state is (*S, I*). In particular, *u*∗ gives the optimal way to start and stop the NPIs. ### Optimal NPIs exist without reducing the reproduction number below one Our first main result is a complete analytical characterization of the optimal NPIs in the SIR model (see Box 1 for a summary and Supplementary Note S1 for details). To understand how the optimal NPIs work, note that the SIR model predicts a *safe zone* of states (*S, I*) where, without any further interventions, the disease prevalence will not exceed *I*max (blue zone in Fig. 2a-c). The safe zone is characterized by the inequality ![Graphic][1], where *R* is the *basic reproduction number* of the outbreak in the population, and the function *Φ**R* is defined in Eq. (2) of Box 1. The goal of an optimal NPI is thus to reach this safe zone as fast as possible without violating the restriction *I*(*t*) ≤ *I*max. The ability to achieve this goal depends on the epidemic state. That is, we can partition the plane (*S, I*) in two regions: those states from which it is possible to reach the safe zone without exceeding *I*max (*feasible* states), and those where it is impossible (*unfeasible* states). We find these two regions are characterized by the separating curve ![Graphic][2], where we call *R**c* := (1 − *u*max)*R* the *controlled reproduction number* (Fig. 2a-c). Note that *R**c* describes the maximum reduction in the basic reproduction number that (constant) admissible interventions can achieve. Therefore, *R**c* < 1 is the necessary and sufficient condition that a constant and permanent admissible intervention (i.e., *u*(*t*) ≡ const. for all *t* ≥ 0) needs to satisfy to *eradicate* a disease outbreak in the SIR model. However, for outbreak mitigation, our analysis shows that feasible states exists without achieving disease eradication (white regions in Fig. 2b-c). This result is important because it proves that optimal NPIs for epidemic mitigation do not require reducing the basic reproduction number below one. ![Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F2.medium.gif) [Figure 2:](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F2) Figure 2: Existence of non-pharmaceutical interventions in the Susceptible-Infected-Removed model. Parameters are *γ* = 1/7, *β* = 0.52, (i.e., *R* = 3.64) and *I*max = 0.1. The safe zone (in blue) consists of all states that dot not exceed *I*max without interventions. This zone is characterized by the inequality ![Graphic][3]. The plane is further divided into feasible states that can reach the safe zone without exceeding *I*max (white), and unfeasible states that cannot (gray). Feasible and unfeasible states are separated by the separating curve ![Graphic][4] (black line). **a**. For “strong” interventions with *u*max = 0.8, the controlled reproduction number is *R**c* = (1 − *u*max)*R* = 0.728 < 1. Here, the separating curve is the straight line *I*max, implying that all states below *I*max are feasible. Note this case corresponds to eradication. **b**. For “intermediate” interventions with *u*max = 0.6, the controlled reproduction number is *R**c* = (1−*u*max)*R* = 1.456 > 1. Here, the separating curve ![Graphic][5] is nonlinear, and some states below *I*max are unfeasible. **c**. For “weak” interventions with *u*max = 0.4 we obtain *R**c* = 2.184 > 1. In this case, states with *S*(0) ≈ 1 are unfeasible. **d**. For *S*(0) → 1, our design criterion for NPIs prescribe the values of *R**c* ‘s that a given *I*max can manage. ### A design criterion for optimal NPIs We demonstrated above that optimal NPIs exist even when *R**c* > 1. However, how large can *R**c* be before NPIs keeping *I*(*t*) ≤ *I*max do not exist? When *S*(0) → 1, our characterization shows that an NPIs exists if and only if ![Formula][6] The above inequality is our second main result, connecting the specified maximum disease prevalence *I*max with the outbreak’s controlled reproduction number *R**c* = (1 − *u*max)*R* (Supplementary Note S2). The inequality (1) governs the existence of NPIs for mitigating epidemic outbreaks, in analogy to how the condition *R**c* < 1 works for disease eradication. Note that *R**c* < 1 is a sufficient condition for the existence for NPIs, but the inequality (1) shows that this condition is far from necessary. If *I*max > 0, there exists *R**c* > 1 for which NPIs exist (Fig. 2d). Note also that the maximum feasible *R**c* increases with *I*max. We can use (1) to design NPIs as follows. Consider an infectious disease outbreak with a given *R* and that the specified maximum prevalence is *I*max. Then, the inequality (1) gives the criterion to design NPIs by providing the range of disease transmission reduction *u*max that the NPIs should attain. In particular, it provides the minimal reduction ![Graphic][7] in the contact rate required for the existence of NPIs. For example, if *I*max = 0.1 then ![Graphic][8] is the maximum admissible controlled reproduction number (orange point in Fig. 2d). Therefore, if an outbreak in the population has *R* = 3, then the minimal reduction is ![Graphic][9] because ![Graphic][10]. ### Optimal NPIs are simple For any epidemic state, the optimal transmission reduction takes a simple form which can be described by coloring the (*S, I*) plane, see top row of Fig. 3. Here, for all states in the white region the optimal intervention is no intervention; for all states in the yellow region the optimal intervention is *u*∗(*S, I*) = *u*max. There are regions (specifically lines) where the optimal intervention switches frequently between *u*∗ = 0 and *u*∗ = *u*max producing a so-called “singular arc” that slides along the two regions, leading to an “average” intervention *u*∗ ∈ [0, *u*max]. We find that, in general, the optimal NPIs have four phases: a first one where no intervention is needed, a second phase where interventions start with maximum strength, a third phase of gradual decrease of interventions, and a “final push” where the maximum interventions are re-applied for a short period to reach the safe zone faster. ![Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F3.medium.gif) [Figure 3:](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F3) Figure 3: Optimal non-pharmaceutical interventions in the Susceptible-Infected-Removed model. For all panels, the parameters of the SIR model are *γ* = 1/7, *β* = 0.52, (i.e., *R* = 3.64) and *I*max = 0.1. We consider a population of *N* = 8.855 x 106 individual (like in Mexico City) and *I* = 1*/N*. Panels shows trajectories for three initial proportions of the susceptible population: large *S* = 1 − *I* ≈ 1 (pink), medium *S* = 0.8 (green), and small *S* = 0.65 (purple). **a**. For *u*max = 0.8 we have *R**c* = (1 − *u*max)*R* = 0.728 ≤ 1. In this case, the optimal intervention starts when the disease prevalence reaches *I*max. Afterwards, the intervention decreases in an hyperbolic arc until reaching the point *S* = *S*∗. At that time, the intervention becomes maximum in the “final push” to reach the safe zone. **b**. For *u*max = 0.58 the controlled reproduction number is *R**c* = (1 − *u*max)*R* = 1.52 > 1. Here ![Graphic][11], implying that the epidemic still can be mitigated for initial states with *S* ≈ 1 and *I* ≈ 0 (pink trajectory). In this case, the optimal intervention starts when the initial condition hits the separating curve below *I*max at *t* = 35. At that instant the intervention starts with the maximum value *u*max, and continues in that form until the trajectory reaches *I*max. **c**. Choosing *u*max = 0.4 yields *R**c* = 2.184 > 1. In this case, the optimal intervention problem does not have a solution for all initial states *S* > 0.85. This is illustrated by pink trajectory: even when applying the maximum intervention from the start, *I*(*t*) will grow beyond *I*max. We illustrate the above behavior in three qualitatively Different cases. The first case is when the optimal intervention starts just when the disease prevalence reaches *I*max (Fig. 3a). This case occurs when the interventions are strong enough to stop the rise in prevalence at *I*max regardless of the fraction of susceptible population. Our analysis shows that this occurs if and only if *u*max is large enough to render *R**c* = (1 − *u*max)*R* ≤ 1. When the initial susceptible population is close to 1 (pink trajectory in Fig. 3a), the optimal intervention first waits until the disease prevalence reaches *I*max. At that time, the optimal NPI stops the disease prevalence exactly at *I*max, and then it gradually decreases its magnitude to ensure that the disease prevalence slides along *I*max as the susceptible population decreases. When the susceptible population reaches the threshold *S*∗, the optimal intervention is again the maximum one (Fig. 3a). This “final push” allows reaching the safe zone faster, releasing the interventions sooner. The middle and bottom panels of Fig. 3a show the resulting disease prevalence and optimal interventions as a function of time. Note that a smaller initial susceptible population yields other trajectories (green and purple in Fig. 3a). The second case is when an “early” intervention is necessary before the disease prevalence reaches *I*max (Fig. 3b). This case happens when the admissible reduction in the contact rate cannot immediately stop the disease prevalence at *I*max if the susceptible population is large at that time. We find this case occurs if and only if *u*max is small in the sense that *R**c* = (1 − *u*max)*R* > 1. Here, a trajectory may hit the yellow region before reaching *I*max (pink trajectory in Fig. 3b). When that happens, the optimal intervention starts with the maximum reduction *u*∗ = *u*max. Then it maintains this maximum reduction to “slide” the trajectory between the yellow and white regions. Once the trajectory reaches *I*max, the magnitude of the optimal intervention decreases to slide the trajectory along *I*max. Again, the final push occurs when the susceptible population reaches *S*∗. The third case is when the initial state (*S*; *I*) lies in the unfeasible region (Fig. 3c). This case occurs when *u*max is so small that, even if the maximum admissible intervention *u* = *u*max is applied from the start of the outbreak, the disease prevalence will exceed *I*max (pink trajectory in Fig. 3c). In this case the optimal intervention problem is unfeasible because it is impossible to achieve *I*(*t*) ≤ *I*max. However, note that the using *u*∗ = *u*max yields the smallest prevalence peak. ### Optimal NPIs are robust To evaluate the optimal NPIs in more realistic scenarios, we numerically analyzed their performance in three epidemic models with uncertain epidemic parameters and more detailed epidemic dynamics (see details in Supplementary Note S3). In all cases, we consider that the basic reproduction number has been estimated as ![Graphic][12] using an SIR model, and that the optimal NPIs are designed using this estimate. Then, these optimal NPIs are applied to an outbreak with possibly Different epidemic dynamics and possibly Different *R*. Note that estimation errors in *R* will affect the correct start and “final push” for reaching the safe zone. In the first scenario, we consider an outbreak with SIR dynamics where the strength of the NPIs is uncertain. We model this uncertainty replacing *u* by *ku* in the model equations, where *k* ∈ (0, 1). Then, for example, *k* = 0.9 (resp. *k* = 1.1) represents a 10% underestimation (resp. overestimation) of the NPIs strength. Across outbreaks with Different *R*’s and an uncertainty of 10% in the intervention’s strength, we find that the disease prevalence is maintained below *I*max as long as *R* is not underestimated (Fig. 4a). In the second scenario, we consider an SEIR outbreak with an incubation period for the disease. For an incubation period of 7 days as in a typical COVID-19 infection, the optimal NPIs maintain the disease prevalence below *I*max if *R* < 2.5 and its value is estimated with an error of below 30% (solid yellow and orange in Fig. 4b). For larger *R* or a larger incubation period, the disease prevalence may exceed *I*max (red in Fig. 4b). ![Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F4.medium.gif) [Figure 4:](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F4) Figure 4: Optimal non-pharmaceutical interventions are robust. For all panels, the estimated parameters used for constructing the optimal NPIs are ![Graphic][13]. We consider a population of *N* = 8.855 x 106 as in Mexico City, and the initial conditions *I*(0) = 1*/N* and *S*(0) = 1 − 1*/N*. If the models contain other state variables, they were initialized at zero. The optimal NPIs are constructed assuming ![Graphic][14], while the actual epidemic dynamics has a possibly Different *R* = *β/γ*. Panels shows results for outbreaks with three values of *R*: low (yellow), medium (orange), and large (red). **a**. A SIR model where the reduction in the disease transmission by the NPIs is uncertain. We model this case replacing *u* by *ku* in the model equations. Panel shows the results for *k* = 1.1 (dotted), *k* = 1 (solid), and *k* = 0.9 (dashed). **b**. SEIR model where exposed individuals do not transmit the infection, with *λ* > 0 the incubation period. Panel shows the results for *λ* = 1/5 (dotted), *λ* = 1/7 (solid), and *λ* = 1/11 (dashed). **c**. A SEIIR model with *λ* = 1/7 and two classes of infected individuals (symptomatic and asymptomatic). Here, *p* ∈ [0, 1] is the proportion of exposed individuals that become asymptomatic. The vertical axis denotes the disease prevalence for symptomatic individuals. The panel shows the results for *p* = 0.55 (dotted), *p* = 0.7 (solid), and *p* = 0.8 (dashed). For the final scenario, we consider an SEIIR model with an incubation period of 7 days and with a fraction *p* ∈ [0, 1] of infected individuals that are asymptomatic and thus remain hidden to the epidemic surveillance system. The goal is to maintain the prevalence of symptomatic individuals below *I*max, without knowing the fraction of asymptomatic individuals. This situation occurs during the COVID-19 pandemic, where between *p* = 0.55 and *p* = 0.8 of infections are asymptomatic [12]. For *p* < 0.7 and *R* < 3.64, the optimal NPIs maintain the disease prevalence of symptomatic individuals below or very close to *I*max if the estimation error for *R* is below 30% (dotted and solid lines in Fig. 4c). An outbreak with low *R* produces a maximum disease prevalence of symptomatic individuals below *I*max, which may result in a larger effective duration of the interventions. Overall, these numerical results shows that the optimal NPIs are robust against a wide range of parameter uncertainty and unmodeled dynamics, provided that the estimation error in the outbreak’s basic reproduction number does not exceed 30%. ## Designing optimal NPIs for the COVID-19 pandemic To explore the implications of our simple criterion for designing NPIs, we analyzed how 16 cities and regions implemented NPIs during the COVID-19 pandemic. For each region or city, we constructed *I*max using the number of available intensive care beds, considering that a fraction of the infected individuals will require them (Supplementary Note S4). The *I*max we obtain ranges from 2.87 x 10−3 for Lima (Peru) to 109.78 x 10−3 for Boston (US), reflecting the large heterogeneity of the available health services across the globe (Fig. 5a). With this information, we calculated the maximum feasible ![Graphic][15] for each region using our design criterion of inequality (1). Since ![Graphic][16] is a monotone function of *I*max, we find that ![Graphic][17] follows the same trend as *I*max (Fig. 5b). The smallest ![Graphic][18] occurs for Lima and the largest ![Graphic][19] for Boston. Note that in both cases ![Graphic][20]. This result implies that, for the *R* of a region’s disease outbreak, NPIs policies must be implemented to guarantee that at least a reduction ![Graphic][21] can be achieved such that ![Graphic][22]. ![Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F5.medium.gif) [Figure 5:](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F5) Figure 5: Minimum necessary reduction in disease transmission for NPIs in the COVID-19 pandemic. **a**. Calculated *I*max according to the proportion of available intensive care beds in each region or city and the estimated fraction of infected individuals requiring intensive care. **b**. Maximum controlled reproduction number *R**c* that each region or city can handle according to its *I*max. Larger *I*max allows a larger *R**c*. **c**. Basic reproduction number *R* per region or city before interventions started. Median (blue big dot), and 95% confidence interval (smaller dots) are shown. **d**. Minimum *u*max necessary for feasibility for each region or city (blue) according to the *R* of panel c. Grey bars denote the reported average mobility reduction in each region between March 19 and April 30. Next, we investigated the *minimal* reduction ![Graphic][23] in transmission required to achieve those upper bounds for the COVID-19 pandemic. For this, we first collected information for the *R* in each region calculated at the start of the pandemic and when the NPIs were inactive (Supplementary Note S3). We find a median nominal *R* of 2.2, with Tokyo having the smallest one (*R* = 1.3) and Madrid having the largest one (*R* = 3.11), see Fig. 5c. From these values of *R*, we calculated the minimal required reduction ![Graphic][24] per region or city (blue in Fig. 5d). For the nominal *R*’s per region or city, we find that a median reduction of ![Graphic][25] of 0.42 is necessary. However, this minimal necessary reduction is heterogeneous across regions. For example, Tokyo just requires ![Graphic][26] while Madrid requires ![Graphic][27]. These two cities have the smallest and largest *R*, respectively. If two cities have a comparable *R*, then the city with large *I*max ends requiring a smaller ![Graphic][28] (e.g., Boston with ![Graphic][29] and Lima with ![Graphic][30]). To evaluate the feasibility of achieving the minimal reduction predicted by our analysis, we collected data for the average mobility reduction in each region during the NPIs in each region (grey in Fig. 5d and Supplementary Note S4). Considering this average mobility reduction as a proxy for the reduction in disease transmission, we find that all regions achieved a greater than necessary reduction. For example, Delhi attained a mobility reduction of 0.84, while the minimal necessary reduction in transmission according to our analysis is ![Graphic][31]. Other regions are in the boundary. For example, New South Wales attained a mobility reduction of 0.48, while the minimal necessary reduction in transmission was ![Graphic][32]. Overall, across regions, we find a median excess of 0.22 in the reduction of mobility compared to the minimal reduction in transmission ![Graphic][33] predicted by our analysis. ## Discussion and concluding remarks Our choice of a simple SIR model was motivated by its epidemiological adequacy for the COVID-19 pandemic and its low dimensionality. The mathematical tractability of the SIR model gives us a complete understanding of the optimal NPIs to apply at any epidemic state. The feedback form *u*∗(*S, I*) of the optimal intervention reflects such understanding, prescribing the optimal action to perform if the epidemic is at state (*S, I*). This feedback strategy should be contrasted to most other studies applying optimal control to epidemic outbreaks, where the optimal intervention is written as an open-loop function of time *u*∗(*t*) [13–16] (see Supplementary Note S4 for details about how our work is related to existing optimal control studies). The open-loop intervention gives the optimal action at any time for a particular initial state. However, it does not tell us what the optimal action is if the epidemic is not in the exact state predicted by the model. Understanding the optimal action to perform at any state has the crucial advantage of allowing us to apply this knowledge to any model, and therefore to reality. Indeed, feedback gives control strategies the required robustness to work on real processes [17, 18], and we numerically confirm that the optimal NPIs we derived have such robustness. Future work could analyze the robustness of the optimal intervention when the state of the epidemic is not entirely known. For example, this case may happen when significative delays exist in reporting new infections, or when tests for identifying infected individuals are limited. The optimal intervention resulting from our analysis can take a continuum of values that may be infeasible to implement in practice. We can use an averaging approach to circumvent this problem. Namely, consider a time window of *T* days (e.g., a week). Suppose that the average reduction prescribed by the optimal intervention over a certain window is *ū*∗. We can realize this reduction on average by combining *d* = *T ū*∗/*u*max days of maximum reduction with (*T* − *d*) days without intervention. This approach yields an intervention similar to Karin et al. [19], with the difference that the periods of intervention and activity are optimally balanced. Our criterion to design optimal NPIs for mitigating epidemic outbreaks is obtained by characterizing the necessary and sufficient conditions for the existence of solutions to an optimal control problem. Specifically, the low-dimensionality of the SIR model allowed us to apply Green’s Theorem to compare the cost of any two interventions analytically (Supplementary Note S1.4). In this sense, the method we use to derive the optimal NPIs is closer to our previous work on optimal control for bioreactors [20]. In general, deriving such complete characterization of optimal control problems is challenging because it involves solving an infinite-dimensional optimization [21]. Indeed, computational methods cannot produce such a characterization [22], and established analytical methods like Pontryagin’s Maximum Principle only yields necessary conditions for optimality [21]. We note that there are several studies applying these and other similar methods to the SIR model [23, 24], in particular during the COVID-19 pandemic [11, 25–28]. Our results could guide a complete characterization of optimal NPIs for more detailed epidemic models or more detailed optimization objectives, but this is likely very challenging. We will inevitably face new epidemics where non-pharmaceutical interventions are the only option to control infections. Rather counter-intuitively, we find that for “ending” an epidemic outbreak as fast as possible using NPIs it is not always optimal to apply the maximum intervention. This observation illustrates the need for developing a better scientific understanding that can inform the design of optimal non-pharmaceutical interventions and plan the required health services capacity. BOX 1. **Optimal NPIs for the Susceptible-Infected-Removed (SIR) model** The SIR model with interventions *u*(*t*) ∈ [0, *u*max] reducing disease transmission takes the form ![Formula][34] Here, *S*(*t*) and *I*(*t*) are the proportion of the population that is susceptible or infected at time *t* ≥ 0, respectively. We denote by (*S*; *I*) the initial state at *t* = 0. The parameters of the SIR model are the (effective) *contact rate β* ≥ 0, and the mean *residence time* of infected individuals *γ* ≥ 0 (in units of day−1). By assuming *S* ≈ 1, these two parameters yield the *basic reproduction number R* = *β/γ*. We are interested in reaching the *safe zone* ![Formula][35] where ![Formula][36] The safe zone is the largest set with the following property: If, for any given time *t*1, the state (*S*1; *I*1) belongs to *𝒮*, we can set *u* = 0 henceforth and still have *I*(*t*) ≤ *I*max for all *t* ≥ *t*1. That is, when *𝒮* is reached, we can terminate the intervention with the assurance that a possible rebound in the disease prevalence will not exceed *I*max. Our goal is to steer an arbitrary initial state (*S*; *I*) to the safe zone *𝒮* in minimal time without violating the constraint *I*(*t*) ≤ *I*max. We say that an intervention achieving this goal is an *optimal intervention*. In Supplementary Note S1, we prove that the existence of an optimal intervention is characterized by the *separating curve* ![Graphic][37] as follows: 1. An optimal intervention exists if and only if the initial state (*S*; *I*) lies below this separating curve (i.e., ![Graphic][38]). Above, *R**c* := (1 − *u*max)*R* is the *controlled reproduction number*. Moreover: 2. If it exists, the optimal intervention *u*∗ at the state (*S, I*) is ![Formula][39] with ![Formula][40] Above, the curve *S* = *Ψ* (*I*) is defined in Supplementary Note S1, while *S* ∗ denotes the intersection of *S* = *Ψ* (*I*) and ![Graphic][41]. ## Data Availability Data is included in Supplementary Note 3 ## Supplementary Notes ### S1. Characterization of the optimal intervention in the Susceptible-Infected-Removed model The model is given by ![Formula][42] where the parameters *β* > 0, *γ* > 0 are assumed constant. Since the total population *N* = *S* + *I* + *R* remains constant all the time, the model can be reduced to that of a second order system using only the states (*S, I*). The maximal (acceptable) value of *I* is *I*max and the maximal achievable value of the control is *u*max. So the state has to belong to the following feasible sets ![Formula][43] Sometimes it will be useful to write the Differential equation in a compact form as ![Formula][44] The trajectory starting at the initial point *x* = (*S*, *I*) and subject to the control *u* : ℝ → 𝒰*F* is denoted by *ϕ* (*t, x*, *u* (·)). Let us define the function ![Formula][45] with *R**α* ∈ {*R**c*, *R*}. The optimal control problem consists in finding the control strategy *u* such that, starting from the initial point (*S*, *I*), the target set ![Formula][46] is reached in the minimal time with the state restriction *I* (*t*) ≤ *I*max satisfied for all time. Note that this set is positively invariant without control (*u* = 0), and that every trajectory that starts in this set satisfies the restriction *I*(*t*) ≤ *I*max for all *t* ≥ 0 (see Fig. S1). Now let us define the reachable set for an initial state *x* as the set of points that can be reached from the initial point *x* with feasible control, i.e., ![Formula][47] Also, we define the controllable set of the target set 𝒯 as the set of points from which some point in the target 𝒯 can be reached with a feasible control, i.e., ![Formula][48] The set 𝒞 (𝒯) can be equivalently described as ℛ (𝒯) for the system ![Formula][49] i.e., the set of points that can be reached from the set 𝒯 for the dynamics with backward time. Now, the optimal control problem has a solution if and only if ![Formula][50] Since the points of the form (*S, I*) = (*S*, 0) are equilibria for every control value, ℛ ((*S*, 0)) = (*S*, 0), we exclude them from the initial conditions for which there is a solution (except if the equilibrium is already in the target set). Now, since ![Graphic][51] for *S* > 0, *I* > 0, ![Formula][52] for every initial condition (except for initial conditions of the form (*S*, 0)). It is obvious that, for the problem to be feasible, the initial state has to be in the feasible set 𝒳*F*, i.e., ![Formula][53] ![Supplementary Figure S1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F6.medium.gif) [Supplementary Figure S1.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F6) Supplementary Figure S1. The set 𝒯 is the largest positive invariant set satisfying *I* ≤ *I*max. The figure was generated with *R* = 2, *R**c* = 1.18 and *I*max = 0.02. #### S1.1 Calculation of the orbits Although it does not seem to be possible to find the trajectories of the system explicitly, it is easy to find its orbits. For this we write (we exclude the points for which *I* = 0 since they are equilibria) ![Formula][54] which is a separable Differential equation (DE). Assuming that *u* is constant and integrating, we obtain ![Formula][55] An interesting rewriting of (S1) is ![Formula][56] This means that the quantity ![Graphic][57] remains constant along the trajectory. Note that this constant depends on the control value used. The above equation is well-known for the SIR model (see, e.g., [1]). Given an initial condition (*S*, *I*) this expression gives, for any 0 < *S* < *S* the (unique) value of *I* that is reached in future time1. Thus there exists a function *I* (*S*; (*S*, *I*)) that gives the value of *I* as a function of *S* and the initial condition. Moreover, from the first equation in the DE we obtain ![Formula][58] and, if we take the expression *I* (*S*; (*S*, *I*)), we obtain a separable DE that can be integrated, ![Formula][59] and that gives the time to reach the point (*S, I* (*S*)) from the initial point (*S*, *I*) with the (constant) control *u*. Although it does not seem possible to give an explicit expression for this integral, it is clear that *S* parametrizes uniquely the solutions (since it is monotone). #### S1.2 The number of infected people If we apply a constant control 0 ≤ *u* ≤ *u*max the infection will eventually die out, i.e., the value *I* (∞) = 0 will be reached asymptotically (otherwise *R* (*t*) would continue growing, which is impossible). We can therefore compute *S* (∞) implicitly from (S1) as ![Formula][60] or, equivalently, as ![Formula][61] Note that the final value of *S* depends on the initial values, but also on the control used. If we assume that the model is normalized, and the initial value is *S* = 1 and *I* ≈ 0, then ![Formula][62] Note that, if *u* → 1−, then *S* (∞) → 1−. So, the larger the value of *u*, the larger the value of *S* (∞). #### S1.3 Reachable set from (*S*, *I*) At each point in the state space, the directions in which the vector field points for Different values of the control are given by *F**u* (*x*) = *f* (*x*) + g (*x*) (1 − *u*). The extreme values are given by *F* (*x*) = *f* (*x*) + g (*x*) and ![Graphic][63], ![Formula][64] ![Supplementary Figure S2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F7.medium.gif) [Supplementary Figure S2.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F7) Supplementary Figure S2. Phase Plane with maximal and minimal orbits bounding the reachable set ℛ (*x*). Max corresponds to the trajectory *ϕ* (*t, x*, *u* = 0) while Min to *ϕ* (*t, x*, *u* = *u*max). The figure was generated with *R* = 2, *R*c = 1.18 and *I*max = 0.02. In the phase plane (*S, I*) both point to the “left”, since the first component (in the direction of *S*) is always negative (recall that *SI* > 0). Since for the second components of the vector fields we have ![Formula][65] it follows that *F* is “above” ![Graphic][66]. Therefore, the reachable set ℛ (*x*) is bounded by the two trajectories *ϕ* (*t, x*, *u* = 0) and *ϕ* (*t, x*, *u*max), see Fig. S2. These two bounding orbits can be easily calculated using Eq. (S1). #### S1.4 Comparing the cost of two Different trajectories In order to be able to find the orbit (trajectory) solving the optimal control problem, it is necessary to be able to compare the cost of two Different trajectories that start at the same initial point and end at the same final point. Consider two orbits *ω**i* (*x*, *x**f*, *u**i*), *i* = 1, 2, joining the (same) points *x* and *x**f* using two Different control actions, *u*1 and *u*2, respectively. The cost (*i*.*e*. time) going through *ω**i* is ![Formula][67] along the trajectory. Given two such orbits, we want to compare both costs. This can be done, for example, by subtracting them, i.e., if ![Formula][68] then the cost of *ω*1 is lower than that of *ω*2. The cost *J* (*u**i*) can be calculated as a line integral along the trajectory. We can see this in the following manner. Calculate ![Formula][69] Now, by properties of the determinant this is also the same as ![Formula][70] Therefore, ![Formula][71] which is a line integral along the orbit *ω**i*. Since the two paths have the same initial and final points, they form a closed curve, and calculating the line integral along the closed curve followed in the counterclock-wise direction we obtain the difference of the costs, i.e. ![Formula][72] where Γ is the closed path of the two orbits followed in the counterclockwise direction. For this we have to assume that: (1) the two paths (orbits) do not intersect at any points except the initial and final ones, and (ii) that Δ ≠ 0. Using Green’s theorem, the line integral can be calculated using a surface integral: ![Formula][73] where ℛ is the region enclosed by the closed curve Γ. For our problem this becomes ![Formula][74] In our case, ![Formula][75] We see that *w* < 0 everywhere, and therefore the integral is always negative, implying that the “upper” orbit has a lower cost than the “lower” orbit (in the closed path traversed in the counterclockwise direction). This observation allows us to find the optimal orbit by comparing it with others. #### S1.5 Optimal orbits From the previous results, the “upper” trajectory is the one with no control (*u* = 0) and, in terms of the cost alone, this trajectory is better than any other one joining the same two points. However, such control may be inadmissible, since the corresponding *I* can go over *I*max at some periods of time. The computation of the optimal control can be approached in two ways: * Fix the initial condition *x*, find its optimal orbit and then its associated optimal control. * Study the optimal control problem for all possible initial conditions. Although the second approach is obviously better, it is more difficult, so we will start with the first approach. In fact, both approaches should lead to the same conclusions. Now we can divide the study of the optimal orbit in several cases. ##### S1.5.1 Unfeasible trajectories This is the case if *I* > *I*max. ##### S1.5.2 Trivial trajectories This is the case in which ![Graphic][76]. That is, the case in which we start in the target set. ##### S1.5.3 Bang-bang trajectories If *x* ∉ 𝒯, it is necessary to apply some control to maintain *I* below the maximal value *I*max. Moreover, admissible trajectories necessarily cross the boundary of 𝒯 at *S* ≥ 1/*R*, that is, they enter 𝒯 at ![Formula][77] In order to find the optimal control that steers an initial state *x* to *x**f* ∈ ∂𝒯1, consider the change of coordinates ![Formula][78] with inverse ![Formula][79] ![Formula][80] Note that, given *I* = 0, we can uniquely map *μ* to *S* ≥ 1/*R*, which we will denote by *S* = *Ŝ**μ*(*μ*). Likewise, we can uniquely map *ν* to *S* ≥ 1/*R**c*, and we will denote it by *S* = *Ŝ**ν* (*ν*). In the new coordinates, the dynamic equations are ![Formula][81] ![Formula][82] Note that *μ* = const is an orbit when *u* = 0 and *ν* = const is an orbit when *u* = *u*max. In (*μ, ν*)-coordinates, the entry point at 𝒯 is the segment ![Supplementary Figure S3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F8.medium.gif) [Supplementary Figure S3.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F8) Supplementary Figure S3. An optimal bang-bang trajectory that initiates at *x* with *u* = 0 (red). When the state touches the curve 𝒮 (green), the control is set to *u* = *u*max until the state reaches *x* f ∈ ∂𝒯1. The region containing all bang and bang-bang trajectories is depicted in yellow. The figure was generated with *R* = 2, *R*c = 1.18 and *I*max = 0.02. ![Formula][83] with ![Formula][84] and ![Formula][85] (see Fig. S3). It follows from Sec. S1.4 that the fastest orbit joining an initial state (*μ*, *ν*) ∈ 𝒞(∂𝒯1) and a final state (*μ**f*, *ν**f*) ∈ ∂𝒯1 is the concatenation of a first piece connecting (*ν*, *μ*) and (*ν**f*, *μ*) with *u* = 0 and a second piece connecting (*ν**f*, *μ*) and (*ν**f*, *μ**f*) with *u* = *u*max. That is, the control is bang-bang. It is easy to verify that this control yields the fastest trajectory, as any other trajectory joining (*μ*, *ν*) and (*μ**f*, *ν**f*) is below this one. The transition times can be computed using (S3) as ![Formula][86] and ![Formula][87] so that the total time is *T* (*ν**f*, *μ**f*; *ν*, *μ*) = *T*(*ν**f*, *μ*; *ν*, *μ*) + *T**c*(*ν**f*, *μ**f*; *ν**f*, *μ*). Note that *μ**f* is fixed, but *ν**f* ∈ [*ν*min, *ν*max] is free. We will now find the closest entry point by minimizing *T* over *ν**f*. Set ![Formula][88] ![Supplementary Figure S4.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F9.medium.gif) [Supplementary Figure S4.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F9) Supplementary Figure S4. The minimal time *φ**μ* (*ν*f) it takes an initial state (*ν*min, *μ*) to reach the target point (*ν*f, *μ*f). The curves were generated with *R* = 2 and *R*c = 1.18. Note that the global minimum is well defined (it is unique). fix *μ* ∈ [*μ**f*, *μ*max] and define the map ![Formula][89] **Assumption 1**. *Global minima of* ![Graphic][90] *are unique*. By Weierstrass Theorem, global minima of ![Graphic][91] always exist. The assumption excludes the highly degenerate case in which the global minimum could occur for more than one value of *ν**f*. Figure S4 shows plots of ![Graphic][92] for various values of *μ* using the parameters *R* = 2 and *R**c* = 1.18. Note that the global minimum is unique (indeed, for large values of ![Graphic][93] is convex). We now define the function ![Formula][94] This function defines a switching curve parameterized by *μ*. In the original coordinates (S2), the switching curve takes the form ![Formula][95] Let *Ī* = max(*S,I*)∈𝒮. To simplify the exposition, we introduce Ψ : [0, *Ī*] → [0, 1], defined implicitly by ![Formula][96] We will parameterize 𝒮 using *I*, ![Formula][97] The trajectories that reach 𝒮 above *I*max are of course unfeasible, so the class of optimal bang-bang trajectories are only those that pass through 𝒮∩𝒳*F* (see the yellow region in Fig. S3). For future reference, we will denote by *S*∗ the *S* coordinate at which 𝒮 intersects the line *I* = *I*max and by *x*1 the point (*S*∗, *I*max). Summarizing, there are two possible situations: 1. **Bang**. If *x* belongs to the region delimited by ∂𝒯1, *I* = *I*max and 𝒮, then the optimal control strategy is simply ![Formula][98] 2. **Bang-Bang**. When *x* belongs to the region delimited by 𝒮, *I* = *I*max and the orbit *ϕ*(−*t, x*1, *u*max), then the optimal control strategy is ![Formula][99] ##### S1.5.4 Trajectories containing a singular arc Let us define an initial point *x* = (*S*, *I*) and the point *x*1 = (*S*∗, *I*max). We are interested in four trajectories (or orbits): 1. *ϕ* (*t, x*, *u* = 0), the trajectory without control starting at *x*. It will be useful to calculate the value ![Graphic][100] at which the orbit (first) touches *I*max. For this we solve (use (S1)) ![Formula][101] for *S* and obtain two solutions: *S*1, *S*2. Define ![Graphic][102] as the largest. Now we calculate the values of *S* ≤ *S**c* for which it is possible to achieve *İ* ≤ 0 (that is, that it is possible to stop the growth of *I*). This value can be calculated from ![Formula][103] and gives ![Formula][104] We “saturate” the value of *S**c* because *S**c* > 1 is not empidemiologically relevant. The control required to achieve the condition *I* = *I*max is the “singular” control ![Formula][105] Note that, if *S* > *S**c* at *I* = *I*max, it is no longer possible to keep *I* at *I*max because *İ* > 0. If ![Graphic][106], then the optimal control is **bang-singular arc-bang**, ![Formula][107] This case is depicted in Fig. S5 for the parameters *R* = 2, *R**c* = 1.18 and *I*max = 0.02. If *S**c* = 1, then ![Graphic][108] holds trivially and the optimal strategy is again (S4). 2. *ϕ* (*t, x*, *u*max), the trajectory with maximal control starting at *x*. 3. *ϕ* (−*t, x*1, *u* = 0), the trajectory without control that passes through *x*1. If *x* is at the left of this trajectory, the optimal orbit is bang-bang, as shown in the previous section. Optimal trajectories starting at the right have singular arcs. 4. *ϕ* (−*t, x*1, *u*∗), the trajectory with control ![Supplementary Figure S5.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F10.medium.gif) [Supplementary Figure S5.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F10) Supplementary Figure S5. An optimal trajectory that initiates at *x* with *u* = 0 (red) and such that ![Graphic][109]. When the state touches the curve *I* = *I*max, the control is set to *u* = *u*sing until the state reaches *x*1. The control is finally switched to *u* = *u*max until the state reaches the target set. The region containing all bang-singular arc-bang trajectories is depicted in yellow. The figure was generated with *R* = 2, *R*c = 1.18 and *I*max = 0.02. ![Formula][110] that passes through *x*1. The control *u*∗ is such that this trajectory does not violate the restriction *I* ≤ *I*max. For values of *S* ≥ *S**c*, it is equal to *u*max, and for *S* ≤ *S**c* it is the control for the singular arc, i.e., it maintains *I* = *I*max until *x**f* is reached. When ![Graphic][111] then it is necessary to start with the control strategy before reaching the maximal value of *I* = *I*max. Otherwise, this limit will be surpassed. However, this is only feasible if, moving backwards from the point (*S**c*, *I*max) with the maximal control *u*max it is possible to reach a point (*S*, *I**c*) such that *I**c* ≥ *I*. The value of *I**c* can be calculated from (S1), ![Formula][112] If *I**c* = *I*, the optimal control is ![Formula][113] When *I**c* > *I*, the control is **bang-bang-singular arc-bang**, ![Formula][114] where (*S**s*, *I**s*) is a switching point. It is characterized as follows: the trajectory *ϕ* (*t, x*, *u* = 0) intersects the trajectory *ϕ* (−*t, x*1, *u*max) at (*S**s*. *I**s*). Such point can be calculated from (S1) as ![Formula][115] By substituting the first into the second we get ![Formula][116] Solving for *S**s* in the second we arrive at ![Formula][117] This case is depicted in Fig. S6 again for the parameters *R* = 2, *R**c* = 1.18 and *I*max = 0.02. If *I**c* < *I*, then it is not possible to solve the optimal problem, since any strategy will surpass the maximal value *I*max. This is the case if, e.g., *u*max is reduced and *R**c* increases to 1.27 (see Fig. S7). #### S1.6 A feedback control strategy The previous “open loop” strategy can be implemented as a state feedback control. This strategy is rather simple, since there are basically only two switching curves: *ϕ* (−*t, x*1, *u*∗) and 𝒮. Another switch takes place when the target region has been attained and the control is switched off, but this happens in a “natural” manner. The switching curve *ϕ* (−*t, x*1, *u*∗) can be written as ![Formula][118] We can further define the “waiting” set ![Formula][119] ![Supplementary Figure S6.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F11.medium.gif) [Supplementary Figure S6.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F11) Supplementary Figure S6. An optimal trajectory that initiates at *x* with *u* = 0 (red) and such that ![Graphic][120]. When the state arrives at (*S**s*, *I**s*), the control is set to *u* = *u*max. When the state touches the curve *I* = *I*max, the control is set to *u* = *u*sing until the state reaches *x*1. The control is finally switched to *u* = *u*max until the state reaches the target set. The region containing all bang-bang-singular arc-bang trajectories is depicted in yellow. The figure was generated with *R* = 2, *R**c* = 1.18 and *I*max = 0.02. ![Supplementary Figure S7.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F12.medium.gif) [Supplementary Figure S7.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F12) Supplementary Figure S7. For the initial condition *x* the problem is unfeasible, *I*max will be surpassed, no matter which strategy is used. The figure was generated with *R* = 2, *R**c* = 1.27 and *I*max = 0.02. ![Supplementary Figure S8.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F13.medium.gif) [Supplementary Figure S8.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F13) Supplementary Figure S8. The optimal feedback strategy (S5). The target and waiting sets, 𝒯 and 𝒲, are illustrated for *R* = 2, *R**c* = 1.18 and *I*max = 0.02. The optimal control feedback is thus given by ![Formula][121] Such strategy is summarized in Fig. S8. Alternatively, we can implement a pure switching control since the “equivalent control”[2] will realize the singular control on the singular arc, ![Formula][122] Note that this control strategy extends the control action beyond the region where the optimal control is feasible. ### S2. Necessary and sufficient conditions for the existence of optimal NPIs Let (*S*, *I*) denote the initial state of the SI model. As shown in Supplementary Note S1, the necessary and sufficient condition for the existence of NPIs is that ![Graphic][123] where ![Graphic][124] is the separating curve. To characterize a condition that is independent of the initial state, we consider the limit case of *S* = 1 and *I* = 0. Under this assumption, the necessary and sufficient condition of existence is that ![Graphic][125]. In other words, the boundary of existence of NPIs is when the separating curve exactly crosses *I* = 0 at *S* = 1. Substituting *S* = 1 in the separating curve we obtain the condition ![Formula][126] which is precisely the inequality (1) of the Main Text.. ### S3. Robustness of the optimal intervention Here we describe the models used to evaluate the robustness of the optimal intervention. #### S3.1 Robustness to the presence of demography and an incubation period To evaluate the robustness of the optimal intervention to the presence of an incubation period of the disease, we considered the SEIR dynamics ![Formula][127] Above, *E*(*t*) denotes the fraction of individuals in the population exposed to the disease, but which are not yet infectious, at time *t*. The parameter 1/*λ* ≥ 0 denotes the *incubation period* of the disease in units of days. The parameter *μ* ≥ 0 denotes the *recruitment rate* in units of days−1. For the result of our paper, we choose *μ* = 1/(365 · 75) corresponding to a life expectancy of 75 years. For this model, the intervention we apply is *u*(*t*) = *u*∗(*S*(*t*), *I*(*t*)) with *u*∗(*S, I*) as in Eq. (S5). #### S3.2 Robustness to the presence of hidden infected individuals To evaluate the robustness of the optimal intervention to hidden infected individuals, consider that that infections can be symptomatic or asymptomatic. We assume that all asymptomatic infections do not require hospital care, and hence remain undetected by the epidemic surveillance system. To model this scenario, we consider the dynamics ![Formula][128] Above, *I**s* denotes the fraction of symptomatic infections and *I**a* the fraction of asymptomatic ones. The model assumes that a fraction *p* ∈ [0, 1] of exposed individuals result in symptomatic infections, and the rest (1− *p*) in asymptomatic ones. We assume that infectious period 1/*γ* is the same for both symptomatic and asymptomatic individuals. For the results of our paper, we choose *λ* = 1/7. Since we assume that only symptomatic individuals end up requiring hospital care, we consider that the objective is to keep *I**s*(*t*) ≤ *I*max only. The control applied is *u*(*t*) = *u*∗(*S*(*t*), *I**s*(*t*)) where *u*∗(*S, I*) is given by Eq. (S5). ### S4. Application to the COVID-19 pandemic #### S4.1 Estimate for the fraction of infected individuals requiring intensive care For COVID-19 pandemic by the SARS-CoV-2 virus, we estimated the fraction *f* of infected individuals requiring intensive-care under the following assumptions: 1. Current estimates for the fraction *p* ∈ [0, 1] of infected individuals that are symptomatic show a large variability [3], ranging from a 20/100 in a report of the World Health Organization, to 96/100 in a study of 328 adults in Shanghai[4]. We take the nominal value of *p* = 60/100. 2. Following Kremer et al.[5], we assume that from the individuals that are symptomatic, a fraction 15/100 develop severe symptoms. 3. Finally, following Li et al. [6], from the individuals that develops severe symptoms, we assume that the fraction 28/100 will require intensive care. Under the above assumptions, the fraction of infected individuals requiring intensive care has a nominal value ![Formula][129] #### S4.2 Data used in our analysis Supplementary Fig. S9 shows the data used for our analysis. Data was collected using the following methodology: 1. **Number of intensive care beds in each city**. This was obtained from official statements when possible (e.g., the Massachusetts Department of Public Health for Boston). In other cases, this number was obtained from public statements of authorities of each city. A complete list of the references appears in the Supplementary Fig. S9. 2. **Population in each city**. Data was obtained from Wikipedia. 3. **Reduction of mobility in each city**. This was obtained from Google Community Mobility Reports [https://www.google.com/covid19/mobility/](https://www.google.com/covid19/mobility/). For our analysis, we considered three categories of mobility: retail & recreation, transit stations, and workplaces. To estimate an overall mobility reduction, we averaged the mobility reduction in these three categories from March 19 to April 30. Data was accessed on May 7, 2020. 4. **Basic reproduction number**. We estimated this quantity from the value of the effective time-varying reproduction number *R**t* at the start of the pandemic around March 8, 2020. These estimates were obtained from the website [https://epiforecasts.io/covid/](https://epiforecasts.io/covid/). ![Supplementary Figure S9.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/10/05/2020.05.19.20107268/F14.medium.gif) [Supplementary Figure S9.](http://medrxiv.org/content/early/2020/10/05/2020.05.19.20107268/F14) Supplementary Figure S9. Table with the response of 16 cities during the COVID-19 pandemic. ### S5. Related work For the control of infectious diseases, there is a large body of work using optimal control methods to design interventions, including vaccination and quarantines[7, 8], drug treatments[9], or dispersal of insecticides and education campaigns[10]. The standard tool to solve these optimal control problem is the celebrated Pontryagin’s Maximum Principle[11]. However, note that the Maximum Principle only gives necessary conditions for optimality. The gap between the necessary and sufficient conditions for optimality needs to be closed using additional arguments, often relying on assuming that the control appears multiplying an affine function of the state variables. This assumption is not satisfied in our formulation of optimal NPIs. We emphasize that the optimal interventions obtained from this approach result in *open loop* strategies which only depend on time. By contrast, our analysis gives a feedback optimal strategy that characterizes the optimal action to make according to the actual state of the epidemic. Indeed, our characterization of optimal NPIs does not rely on the Maximum Principle. Instead, the low dimensional of our model allows us to apply Green’s Theorem to compare the cost of two Different interventions. The consequence of our approach is that we obtain a feedback or *closed loop* strategy that corrects itself based on the actual state of the epidemic. The COVID-19 pandemic has stirred much interest on designing non-pharmaceutical interventions. This has led to strategies like interspacing mitigation with brief periods of activity[12]. Optimal control methods have been also applied, for example to minimize the peak of infection[13], minimize the number of infections[14], minimize the economic costs[15], or maximize welfare[16]. Compared to these studies, our analytical characterization of optimal NPIs provides gives us a complete understanding of the optimal decisions that need to be made. For example, no intervention is needed before reaching the separating curve. ## Footnotes * Revised calculations * 1 If we select *S* > *S* the obtained value of *I* is reached in a past time (*t* < 0). * Received May 19, 2020. * Revision received October 5, 2020. * Accepted October 5, 2020. * © 2020, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NonCommercial-NoDerivs 4.0 International), CC BY-NC-ND 4.0, as described at [http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/) ## References 1. 1.Anderson, R. M. & May, R. M. Vaccination and herd immunity to infectious diseases. Nature 318, 323–329 (1985). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/318323a0&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=3906406&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F10%2F05%2F2020.05.19.20107268.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1985AVC5500047&link_type=ISI) 2. 2.Anderson, R. M. & May, R. M. Infectious diseases of humans: dynamics and control (Oxford university press, 1992). 3. 3.Fraser, C., Riley, S., Anderson, R. M. & Ferguson, N. M. Factors that make an infectious disease outbreak controllable. Proceedings of the National Academy of Sciences 101, 6146–6151 (2004). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMToiMTAxLzE2LzYxNDYiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMC8xMC8wNS8yMDIwLjA1LjE5LjIwMTA3MjY4LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 4. 4.Kupferschmidt, K. Ending coronavirus lockdowns will be a dangerous process of trial and error. Science| AAAS (2020). 5. 5.Ferguson, N. M. et al. Strategies for mitigating an influenza pandemic. Nature 442, 448–452 (2006). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nature04795&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=16642006&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F10%2F05%2F2020.05.19.20107268.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000239278900042&link_type=ISI) 6. 6.Hollingsworth, T. D., Klinkenberg, D., Heesterbeek, H. & Anderson, R. M. Mitigation strategies for pandemic influenza A: balancing conflicting policy objectives. PLoS computational biology 7 (2011). 7. 7.Ferguson, N. et al. Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand. Imperial College COVID-19 Response Team 2020. 8. 8.Anderson, R. M., Heesterbeek, H., Klinkenberg, D. & Hollingsworth, T. D. How will country-based mitigation measures influence the course of the COVID-19 epidemic? The Lancet 395, 931–934 (2020). 9. 9.Elston, J., Cartwright, C., Ndumbi, P. & Wright, J. The health impact of the 2014–15 Ebola outbreak. Public Health 143, 60–70 (2017). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.puhe.2016.10.020&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28159028&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F10%2F05%2F2020.05.19.20107268.atom) 10. 10.Kermack, W. O. & McKendrick, A. G. A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character 115, 700–721 (1927). 11. 11.Djidjou-Demasse, R., Michalakis, Y., Choisy, M., Sofonea, M. T. & Alizon, S. Optimal COVID-19 epidemic control until vaccine deployment. medRxiv (2020). 12. 12.For Evidence-Based Medicine, O. C. COVID-19: What proportion are asymptomatic? <[https://www.cebm.net/covid-19/covid-19-what-proportion-are-asymptomatic/](https://www.cebm.net/covid-19/covid-19-what-proportion-are-asymptomatic/)>. 13. 13.Behncke, H. Optimal control of deterministic epidemics. Optimal control applications and methods 21, 269–285 (2000). 14. 14.Yan, X. & Zou, Y. Optimal and sub-optimal quarantine and isolation control in SARS epidemics. Mathematical and Computer Modelling 47, 235–245 (2008). 15. 15.Rowthorn, R. E., Laxminarayan, R. & Gilligan, C. A. Optimal control of epidemics in metapopulations. Journal of the Royal Society Interface 6, 1135–1144 (2009). 16. 16.Caetano, M. A. L. & Yoneyama, T. Optimal and sub-optimal control in Dengue epidemics. Optimal control applications and methods 22, 63–73 (2001). 17. 17.Åström, K. J. & Murray, R. M. Feedback systems: an introduction for scientists and engineers (Princeton university press, 2010). 18. 18.Angulo, M. T. & Velasco-Hernandez, J. X. Robust qualitative estimation of time-varying contact rates in uncertain epidemics. Epidemics 24, 98–104 (2018). 19. 19.Karin, O. et al. Adaptive cyclic exit strategies from lockdown to suppress COVID-19 and allow economic activity. medRxiv (2020). 20. 20.Moreno, J. Optimal time control of bioreactors for the wastewater treatment. Optimal Control Applications and Methods 20, 145–164 (1999). 21. 21.Lenhart, S. & Workman, J. T. Optimal control applied to biological models (CRC press, 2007). 22. 22.Bonnans, F., Martinon, P. & Grélard, V. Bocop-A collection of examples (2012). 23. 23.Hansen, E. & Day, T. Optimal control of epidemics with limited resources. Journal of mathematical biology 62, 423–451 (2011). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20407775&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F10%2F05%2F2020.05.19.20107268.atom) 24. 24.Di Giamberardino, P. & Iacoviello, D. Optimal control of SIR epidemic model with state dependent switching cost index. Biomedical signal processing and control 31, 377–380 (2017). 25. 25.Alvarez, F. E., Argente, D. & Lippi, F. A simple planning problem for covid-19 lockdown tech. rep. (National Bureau of Economic Research, 2020). 26. 26.Piguillem, F. & Shi, L. Optimal covid-19 quarantine and testing policies (2020). 27. 27.Tsay, C., Lejarza, F., Stadtherr, M. A. & Baldea, M. Modeling, state estimation, and optimal control for the US COVID-19 outbreak. arXiv preprint arxiv:2004.06291 (2020). 28. 28.Morris, D. H., Rossine, F. W., Plotkin, J. B. & Levin, S. A. Optimal, near-optimal, and robust epidemic control. arXiv preprint arxiv:2004.02209 (2020). ## References 1. [1]. H. R. Thieme. Mathematics in population biology, vol. 12 (Princeton University Press, 2018). 2. [2]. L. Fridman, J. Moreno and R. Iriarte. Sliding modes after the first decade of the 21st century. Lecture notes in control and information sciences 412, pp. 113–149 (2011). 3. [3].O. C. for Evidence-Based Medicine. Covid-19: What proportion are asymptomatic? URL [https://www.cebm.net/covid-19/covid-19-what-proportion-are-asymptomatic/](https://www.cebm.net/covid-19/covid-19-what-proportion-are-asymptomatic/). 4. [4]. X. Zhou et al. Follow-up of the asymptomatic patients with sars-cov-2 infection. Clinical Microbiology and Infection. 5. [5]. M. U. Kraemer et al. The effect of human mobility and control measures on the covid-19 epidemic in china. Science 368 no. 6490, pp. 493–497 (2020). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzNjgvNjQ5MC80OTMiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMC8xMC8wNS8yMDIwLjA1LjE5LjIwMTA3MjY4LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 6. [6]. R. Li et al. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (sars-cov2). Science. 7. [7]. H. Behncke. Optimal control of deterministic epidemics. Optimal control applications and methods 21 no. 6, pp. 269–285 (2000). 8. [8]. X. Yan and Y. Zou. Optimal and sub-optimal quarantine and isolation control in sars epidemics. Mathematical and Computer Modelling 47 no. 1-2, pp. 235–245 (2008). 9. [9]. R. E. Rowthorn, R. Laxminarayan and C. A. Gilligan. Optimal control of epidemics in metapopulations. Journal of the Royal Society Interface 6 no. 41, pp. 1135–1144 (2009). 10. [10]. M. A. L. Caetano and T. Yoneyama. Optimal and sub-optimal control in dengue epidemics. Optimal control applications and methods 22 no. 2, pp. 63–73 (2001). 11. [11]. D. Liberzon. Calculus of variations and optimal control theory: a concise introduction (Princeton University Press, 2011). 12. [12]. O. Karin et al. Adaptive cyclic exit strategies from lockdown to suppress covid-19 and allow economic activity. medRxiv. 13. [13]. D. H. Morris et al. Optimal, near-optimal, and robust epidemic control. arXiv preprint arxiv:2004.02209. 14. [14]. C. Tsay et al. Modeling, state estimation, and optimal control for the us covid-19 outbreak. arXiv preprint arxiv:2004.06291. 15. [15]. F. E. Alvarez, D. Argente and F. Lippi. A simple planning problem for covid-19 lockdown. Tech. rep., National Bureau of Economic Research (2020). 16. [16]. F. Piguillem and L. Shi. Optimal covid-19 quarantine and testing policies. [1]: /embed/inline-graphic-1.gif [2]: /embed/inline-graphic-2.gif [3]: F2/embed/inline-graphic-3.gif [4]: F2/embed/inline-graphic-4.gif [5]: F2/embed/inline-graphic-5.gif [6]: /embed/graphic-3.gif [7]: /embed/inline-graphic-6.gif [8]: /embed/inline-graphic-7.gif [9]: /embed/inline-graphic-8.gif [10]: /embed/inline-graphic-9.gif [11]: F3/embed/inline-graphic-10.gif [12]: /embed/inline-graphic-11.gif [13]: F4/embed/inline-graphic-12.gif [14]: F4/embed/inline-graphic-13.gif [15]: /embed/inline-graphic-14.gif [16]: /embed/inline-graphic-15.gif [17]: /embed/inline-graphic-16.gif [18]: /embed/inline-graphic-17.gif [19]: /embed/inline-graphic-18.gif [20]: /embed/inline-graphic-19.gif [21]: /embed/inline-graphic-20.gif [22]: /embed/inline-graphic-21.gif [23]: /embed/inline-graphic-22.gif [24]: /embed/inline-graphic-23.gif [25]: /embed/inline-graphic-24.gif [26]: /embed/inline-graphic-25.gif [27]: /embed/inline-graphic-26.gif [28]: /embed/inline-graphic-27.gif [29]: /embed/inline-graphic-28.gif [30]: /embed/inline-graphic-29.gif [31]: /embed/inline-graphic-30.gif [32]: /embed/inline-graphic-31.gif [33]: /embed/inline-graphic-32.gif [34]: /embed/graphic-7.gif [35]: /embed/graphic-8.gif [36]: /embed/graphic-9.gif [37]: /embed/inline-graphic-33.gif [38]: /embed/inline-graphic-34.gif [39]: /embed/graphic-10.gif [40]: /embed/graphic-11.gif [41]: /embed/inline-graphic-35.gif [42]: /embed/graphic-12.gif [43]: /embed/graphic-13.gif [44]: /embed/graphic-14.gif [45]: /embed/graphic-15.gif [46]: /embed/graphic-16.gif [47]: /embed/graphic-17.gif [48]: /embed/graphic-18.gif [49]: /embed/graphic-19.gif [50]: /embed/graphic-21.gif [51]: /embed/inline-graphic-36.gif [52]: /embed/graphic-22.gif [53]: /embed/graphic-23.gif [54]: /embed/graphic-24.gif [55]: /embed/graphic-25.gif [56]: /embed/graphic-26.gif [57]: /embed/inline-graphic-37.gif [58]: /embed/graphic-27.gif [59]: /embed/graphic-28.gif [60]: /embed/graphic-29.gif [61]: /embed/graphic-30.gif [62]: /embed/graphic-31.gif [63]: /embed/inline-graphic-38.gif [64]: /embed/graphic-32.gif [65]: /embed/graphic-34.gif [66]: /embed/inline-graphic-39.gif [67]: /embed/graphic-35.gif [68]: /embed/graphic-36.gif [69]: /embed/graphic-37.gif [70]: /embed/graphic-38.gif [71]: /embed/graphic-39.gif [72]: /embed/graphic-40.gif [73]: /embed/graphic-41.gif [74]: /embed/graphic-42.gif [75]: /embed/graphic-43.gif [76]: /embed/inline-graphic-40.gif [77]: /embed/graphic-44.gif [78]: /embed/graphic-45.gif [79]: /embed/graphic-46.gif [80]: /embed/graphic-47.gif [81]: /embed/graphic-48.gif [82]: /embed/graphic-49.gif [83]: /embed/graphic-51.gif [84]: /embed/graphic-52.gif [85]: /embed/graphic-53.gif [86]: /embed/graphic-54.gif [87]: /embed/graphic-55.gif [88]: /embed/graphic-56.gif [89]: /embed/graphic-58.gif [90]: /embed/inline-graphic-41.gif [91]: /embed/inline-graphic-42.gif [92]: /embed/inline-graphic-43.gif [93]: /embed/inline-graphic-44.gif [94]: /embed/graphic-59.gif [95]: /embed/graphic-60.gif [96]: /embed/graphic-61.gif [97]: /embed/graphic-62.gif [98]: /embed/graphic-63.gif [99]: /embed/graphic-64.gif [100]: /embed/inline-graphic-45.gif [101]: /embed/graphic-65.gif [102]: /embed/inline-graphic-46.gif [103]: /embed/graphic-66.gif [104]: /embed/graphic-67.gif [105]: /embed/graphic-68.gif [106]: /embed/inline-graphic-47.gif [107]: /embed/graphic-69.gif [108]: /embed/inline-graphic-48.gif [109]: F10/embed/inline-graphic-49.gif [110]: /embed/graphic-71.gif [111]: /embed/inline-graphic-50.gif [112]: /embed/graphic-72.gif [113]: /embed/graphic-73.gif [114]: /embed/graphic-74.gif [115]: /embed/graphic-75.gif [116]: /embed/graphic-76.gif [117]: /embed/graphic-77.gif [118]: /embed/graphic-78.gif [119]: /embed/graphic-79.gif [120]: F11/embed/inline-graphic-51.gif [121]: /embed/graphic-83.gif [122]: /embed/graphic-84.gif [123]: /embed/inline-graphic-52.gif [124]: /embed/inline-graphic-53.gif [125]: /embed/inline-graphic-54.gif [126]: /embed/graphic-85.gif [127]: /embed/graphic-86.gif [128]: /embed/graphic-87.gif [129]: /embed/graphic-88.gif