A simple criterion to design optimal nonpharmaceutical interventions for epidemic outbreaks =========================================================================================== * Marco Tulio Angulo * Fernando Castaños * Jorge X. Velasco * Jaime A. Moreno ## Abstract To mitigate the COVID-19 pandemic, much emphasis exists on implementing non-pharmaceutical interventions to keep the reproduction number below one. But 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. 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 maximum health services’ capacity. To explore the implications of our theoretical results, we study the non-pharmaceutical interventions implemented in 16 cities during the COVID-19 pandemic. In particular, we estimate the minimal reduction of the contact rate in each city that is necessary to control the epidemic optimally. We also compare the optimal start of the intervention with the start of the actual interventions applied in each city. Our results contribute to establishing a rigorous methodology to guide the design of non-pharmaceutical intervention policies. Since the seminal work of May and Anderson1, the design of interventions to *eradicate* infectious diseases has the objective of achieving a basic reproduction number *R* below one.2,3 The underlying assumption here is that interventions like vaccination programs can be maintained for long periods, producing near *permanent* changes to the epidemic dynamics. This same objective is guiding the design of non-pharmaceutical interventions (NPIs) for the COVID-19 pandemic4. Yet, maintaining NPIs like bans of public events or lockdowns for long periods is infeasible because of their substantial economic and societal costs5,6. In this sense, the objective of NPIs cannot be to eradicate a disease; rather, they aim to *mitigate* the economic and social costs of the epidemic outbreak7. However, as evidenced by the world’s controversies during the COVID-19 pandemic, 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. For this, we consider that NPIs should achieve an optimal tradeoff between two objectives8. 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 number of infections does not exceed the health services’ capacity, avoiding unnecessary deaths due to a saturated health system9. We obtain a full analytical characterization of such optimal NPIs, specifying the reduction in the disease transmission that is optimal for each state in which an epidemic can be. Furthermore, this characterization yields the necessary and sufficient criterion for the existence of optimal NPIs, analogous to the *R* < 1 condition. We find that reducing the reproduction number below one is sufficient but not necessary for their existence. Instead, our criterion shows that the maximum health services’ capacity 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 maximum health services’ capacity. We explore the implications of our theoretical result by analyzing the response of 16 cities across the globe to the COVID-19 pandemic. We find that most cities achieved a larger-than-necessary reduction in transmission. We also compare the actual start of the NPIs in each city with the optimal start obtained from our analysis, finding that most cities responded before it was necessary. Our results contribute to understanding how to best respond to an epidemic outbreak using non-pharmaceutical interventions. ## RESULTS Our objective is to characterize the reduction in the disease transmission that is optimal for each *state* in which the epidemic outbreak can be. In the classical Susceptible-Infected-Removed model, the state can be characterized by the pair (*S*, *I*) ∈ [0,1]2, where *I* is the disease prevalence (i.e., proportion of the population that is infected) and S the proportion of the population that is susceptible to the disease (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. To mitigate the epidemic outbreak, consider we can apply one or several NPIs that reduce disease transmission by the factor (1 − *u*), for some *u* ∈ [0,1]. Here, *u* = 0 if the NPIs achieve no reduction, and *u* = 1 if they completely stop transmission. Since it is unfeasible to stop transmission fully, we upper-bound the reduction by *umax* ∈ (0, 1). We say that the intervention u is *admissible* if *u* ∈ [0, *umax*]. To describe the health services’ capacity, we consider they can adequately manage a maximum prevalence *Imax* ∈ (0,1]. For example, hospitals may saturate when the disease prevalence exceeds *Imax*, causing higher mortality10. In this sense, NPIs must maintain the disease prevalence below *Imax*. ### Characterizing optimal non-pharmaceutical interventions In principle, several admissible interventions with different duration can keep the disease prevalence below *Imax*. For instance, the intervention 1 in the example of Fig. 1b-c keeps this restriction and has a duration of 63 days. Intervention 2 also keeps this restriction, but it only lasts 32 days. To design the *optimal* NPI, we ask for the intervention with minimal duration. More precisely, 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 duration of the intervention; and (2) it ensures that the prevalence can be maintained below *Imax* 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 given that the epidemic is in the state (*S*, *I*). The optimal intervention also gives the optimal state to start the NPIs, and optimal way to stop them. ![Figure 1](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F1.medium.gif) [Figure 1](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F1) Figure 1 Optimal non-pharmaceutical interventions **a**. A Susceptible-Infected-Removed (SIR) model with the non-pharmaceutical interventions (NPIs) that reduce disease transmission. For the optimal NPI design problem, the objective is to design the intervention *u**(*t*) with minimal duration such that *u**(*t*) ∈ [0, *umax*] and *I*(*t*) < *Imax* 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 × 10-7 and *S* = 1 − *I*). Intervention 1 guarantees that *I*(*t*) < *Imax* and has a duration of 63 days (dashed line). Intervention 2 also guarantees that *I*(*t*) ≤ *Imax* but its duration is only 32 days (solid line). Actually, intervention 2 is the optimal one derived using our analysis: it is the shortest intervention satisfying *I*(*t*) ≤ *Imax*. For comparison, panel b also shows the response of the model without intervention (i.e., *u*(*t*) ≡ 0). Our first main result is a complete analytical characterization of the optimal NPIs in the SIR model (Supplementary Note S1). We prove that the solution to the optimal intervention is fully characterized by the *separating curve* ![Formula][1] where ![Graphic][2]. Above, *R* is the basic reproduction number of the outbreak without intervention. We define *Rc*:= (1 − *umax*)*R* as the *controlled* reproduction number. Note that *Rc* describes the maximal reduction in *R* that (constant) admissible interventions can achieve. Therefore, *Rc* < 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 in the SIR model. The separating curve is key because it partitions the (*S*, *I*) plane in two regions (top row in Fig. 2). All states below or on the separating curve are feasible, meaning that the optimal NPI problem has a solution. By contrast, all states above the separating curve are unfeasible because the admissible interventions cannot maintain the disease prevalence below *Imax*. Note that the shape of the separating curve depends on *Imax* and *Rc*. If *Rc* ≤ 1, the separating curve is the straight line Ф(*S*) = *Imax*. When *Rc* > 1, the separating curve becomes nonlinear. ![Figure 2](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F2.medium.gif) [Figure 2](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F2) Figure 2 Optimal interventions in the Susceptible-Infected model For all panels, the parameters of the SIR model are *γ* = 1/7, *β* = 0.52, (i.e., *R* = 3.64) and *Imax* = 0.1. We consider a population of *N* = 8.855 × 106 individuals, like in Mexico City, and an initial state *I* = 1/*N* and S = 1 − *I* (white dot with red boundary). a. For *umax* = 0.8 we have *Rc* = (1 − *umax*)*R* = 0.728 ≤ 1. This yields *S** = 1 (green point in the top panel) and the straight separating curve Ф(*S*; *Imax*, *Rc*) = *Imax* (blue line). In this case, the optimal intervention waits until *I*(*t*1) = *Imax* at time *t*1 = 37 (red point). At that time the intervention starts with roughly *u**(*t*1) = 0.7, then decreases in an hyperbolic arc 1 − (*R**S*)−1, finishing at time *t* = 78. The total duration was ![Graphic][3]. b. For *umax* = 0.58 the controlled reproduction number is *Rc* = (1 − *umax*)*R* = 1.52 > 1. This yields *S** < 1 (green point in top panel) and a curved separating curve Ф(*S*; *Imax*, *Rc*) (blue in top panel). Note here that Ф(1) > 0, implying that the epidemic still can be controlled for *S* ≈ 1 and *I* ≈ 0. In this case, the optimal intervention for the initial condition hits the separating curve below *Imax* at time *t* = 35 (red point in top panel). At that instant the intervention starts with the maximal value *umax*, and continues in that form until the trajectory reaches *Imax* (green point). At that instant, the intervention decreases in an hyperbolic arc until vanishing at *t* = 80. The duration of the intervention was ![Graphic][4]. **c**. Choosing *umax* = 0.4 yields *Rc* = 2.184 > 1, leading to *S** < 1 and a separating curve that reaches zero before near *S* = 0.8. In this case, the optimal intervention problem does not have a solution for all initial states *S* > 0.8. This is illustrated by the red initial state. The problem does not have a solution because even by applying the maximum intervention the infected *I*(*t*) will go beyond *Imax*. ### The optimal intervention The separating curve also characterizes the optimal transmission reduction at any state of the epidemic (Box 1). For all states (*S*, *I*) below the separating curve, the optimal intervention is *u**(*S*, *I*) = 0. With this intervention, the trajectory of some states will hit the separating curve at some point (*S*1, *I*1), where *I*1 = Ф(*S*1). When that happens the intervention starts with its maximal value *u**(*S*1, *I*1) = *umax* (i.e. a, maximal reduction of transmission). Then, the optimal intervention “slides” the trajectory along the separating curve. Finally, the optimal intervention stops when the susceptible *S* reaches the region ![Graphic][5]. Note that once the trajectory reaches this region, prevalence will decrease without further intervention. To better understand the optimal intervention, we illustrate its behavior in its three qualitative different cases. The first is when the intervention starts just when the prevalence reaches *Imax* (Fig. 2a). This case happens when interventions are strong enough to maintain *I*(*t*) constant the first time *t*1 they reach *I*(*t*1) = *Imax* (top in Fig. 2a). Our analysis shows that this occurs if and only if *umax* is large enough to render *Rc* = (1 − *umax*)*R* ≤ 1, yielding the constant separating curve Ф(*S*) = *Imax*. In this case, the optimal intervention is “do nothing” until *I*(*t*1) = *Imax*. The optimal reduction *u**(*t*1) in the contact rate is maximal at that time, and then decreases to zero in a hyperbolic arc (bottom in Fig. 2a). The second is when an “early” intervention is necessary before the prevalence reaches *Imax* (Fig. 2b). This case happens when no admissible reduction in the contact rate can immediately stop *I*(*t*) at *Imax* if the susceptible population is large at that time. We find this occurs if and only if *umax* yields *Rc* > 1, leading to *S** < 1 and a separating curve that decreases to the right of *S** (top in Fig. 2b).separating curve Here, the optimal intervention starts when the trajectory hits the separating curve (red point in Fig. 2b). Since the separating curve points is nonlinear, this hitting time marks the “early start” of the intervention. The intervention starts with the maximum reduction *u** = *umax*. Then, it maintains this maximum reduction to slide the trajectory over the separating curve. Once the trajectory reaches *Imax*, the optimal reduction decreases in a hyperbolic arc (bottom in Fig. 2b). The third is when the initial state (*S*, *I*) lies above or to the right of the separating curve (Fig. 2c). This case occurs when *umax* is so small that, even if the maximum admissible reduction maintained from the start, the prevalence will exceed *Imax* (top in Fig. 2c). In this case the optimal intervention problem is unfeasible because it is impossible to achieve *I*(*t*) ≤ *Imax*. Note, however, that the optimal intervention still yields the smallest prevalence peak. ### A simple criterion to design optimal non-pharmaceutical interventions We demonstrated above that the optimal transitory intervention exists even when *Rc* > 1. However, how large can *Rc* be before an optimal NPI does not exist because it is impossible to maintain the prevalence below *Imax*? In the case *I*(0) → 0 and *S*(0) → 1, our characterization shows that an optimal NPIs exists for all (*Rc*, *Imax*) such that the separating curve satisfies Ф(1) ≥ 0 (Supplementary Note S2-1). This observation yields the necessary and sufficient condition for the existence of NPIs: ![Formula][6] The above inequality is our second main result, connecting the health services’ capacity *Imax* with the controlled reproduction number *Rc* = (1 − *umax*)*R* that it can successfully maintain. This inequality governs the existence of optimal NPIs for mitigating epidemic outbreaks, in analogy to how the condition *Rc* < 1 works for the eradication of diseases. Note that the condition *Rc* < 1 is sufficient for NPIs, but the inequality (2) shows that this condition is far from necessary. If *Imax* > 0, there exists *Rc* > 1 for which optimal NPIs exist (Fig. 3a). In these cases, it is possible to maintain *I*(*t*) < *Imax* using an admissible intervention. However, this requires an “early start” of the intervention, as demonstrated in Fig. 2b. Note also that the maximum feasible *Rc* increases with *Imax*. ![Figure 3](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F3.medium.gif) [Figure 3](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F3) Figure 3 Minimum necessary reduction in transmission for the COVID-19 pandemic **a**. The optimal NPI problem is feasible if it is possible to maintain *I*(*t*) ≤ *Imax* for all t ≥ 0. The panel shows the feasibility region determined by Eq. (2). The problem is feasible even when the interventions cannot reduce the basic reproduction below one (i.e., even when *Rc* = (1 − *umax*)*R* > 1). **b**. Duration of the optimal intervention in units of τ = *γt*. Lower *Imax* results in longer duration of the optimal intervention. **c**. Calculated *Imax* according to the proportion of available intensive care beds in each city and the estimated fraction of infected individuals requiring intensive care. **d**. Maximum controlled reproduction number *Rc* that each city can handle according to its *Imax*. Larger *Imax* allows a larger *Rc*. **e**. Initial basic reproduction number *R* per city, before any intervention started. Median (blue big dot), and 95% confidence interval (smaller dots) are shown. **f**. Minimum *umax* necessary for feasibility for each city (blue) according to the *R* of panel c. Grey bars denote the reported average mobility reduction in each city between March 19 and April 30. To understand how to apply (2) for designing NPIs, consider an infectious disease outbreak with a given *R* and that the maximum prevalence that the health services can manage is *Imax*. Then, the inequality (2) gives the criterion to design NPIs by providing the range of contact rate reduction *umax* that the NPIs should attain. In particular, it provides the minimal reduction ![Graphic][7] in the contact rate required for the existence of optimal NPIs. By construction, optimal NPIs have the minimum possible duration. But their duration might be more than is feasible in practice. We investigated how *Rc* and *Imax* change the duration of the optimal intervention (Supplementary Note S2-2). In general, the duration of the optimal NPIs increases as *Imax* decreases (Fig. 3b). This result makes sense as a smaller health system’s capacity requires a stronger flattening of the prevalence curve. The duration also increases as the pair (*Rc*, *Imax*) approaches the feasibility boundary of the inequality (2). In particular, close to this feasibility boundary, the duration is very sensitive to changes in these two parameters. Hence, the pair *umax* and *Imax* should be designed to remain sufficiently far from the feasibility boundary. ### Application to the COVID-19 pandemic We explored the implications of our simple rule for designing NPIs by analyzing how 16 cities implemented NPIs during the COVID-19 pandemic. To estimate the proportion *Imax* of each city, we first collected information about the available intensive care beds in each city during the pandemic. We then used available knowledge of the disease to estimate the fraction of infected individuals who require intensive care (Supplementary Note S3). The *Imax* we obtain ranges from 2.87 × 10-3 for Lima (Peru) to 109.78 × 10−3 for Boston (US), reflecting the large heterogeneity of health services across the globe (Fig. 3c). We calculated the maximum feasible ![Graphic][8] for each city from these quantities using our design criterion (2). Since ![Graphic][9] is a monotone function of *Imax* we find the same trend as in *Imax* (Fig. 3d). The smallest ![Graphic][10] occurs for Lima and the largest ![Graphic][11] for Boston. Note that in both cases ![Graphic][12]. For the *R* of a city’s disease outbreak, NPIs policies must be implemented to guarantee that reduction *umax* can be achieved such that ![Graphic][13]. Next, we investigated the *minimal* reduction ![Graphic][14] in transmission required to achieve those upper bounds for the COVID-19 pandemic case. For this, we first collected information for the *R* in each city 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), Fig. 3e. From these values of *R*, we calculated the minimal required reduction ![Graphic][15] per city (blue in Fig. 3f). For the nominal *R*’s per city, we find that a median reduction of ![Graphic][16] of 0.42 is necessary. However, this minimal necessary reduction is heterogeneous across cities. For example, Tokyo just requires ![Graphic][17] while Madrid requires ![Graphic][18]. These two cities have the smallest and largest *R*, respectively. If two cities have a comparable *R*, then the city with large *Imax* end ups requiring a smaller ![Graphic][19] (e.g., Boston with ![Graphic][20] and Lima with ![Graphic][21]. To evaluate the feasibility of achieving the minimal reduction predicted by our analysis, we collected data for the average mobility reduction in each city during the NPIs in each city (grey in Fig. 3d and Supplementary Note S3). Considering this average mobility reduction as a proxy for the reduction in disease transmission, we find that all cities 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 0.42. Other cities are in the boundary. For example, New South Wales attained a mobility reduction of 0.48, while the minimal necessary reduction in transmission was 0.44. Overall, across cities, we find a median excess of 0.22 in the reduction of mobility compared to the minimal reduction in transmission predicted by our analysis. Finally, we compared the start of the optimal NPI with the actual start of the NPIs implemented in each city (Fig. 4). Most cities (10/16) started their NPIs before the optimal start date, with a median of 40 days of anticipation. Here, Tokyo is an outlier with almost 150 days of anticipation. Four cities (New South Wales, New York, British Columbia, and Boston) started their NPIs at almost the optimal time. Only one city (Madrid) started its NPIs after the optimal start date calculated by our analysis (about ten days of delay). We also found significant heterogeneity in the duration of the calculated optimal NPIs across cities, ranging from about 15 days for Boston to more than one year in Lima. This heterogeneity comes from the heterogeneity in the estimated health system capacity *Imax*. In Lima, for example, the ratio of intensive care units to habitants is 517/8.575 × 106 = 6.02 × 10−5 (i.e., *Imax* = 2.87 × 10−3), while in Boston the ratio is two orders of magnitude larger 1600/694×103 = 2.3×10−3 (i.e., *Imax* = 109.78×10−3). Cities with smaller *Imax* yields optimal NPIs that are longer, up to the point when their duration makes them impossible to implement. Assuming, as illustrative example, that 150 days is the maximal duration of feasible NPIs, our analysis predicts that *Imax* needs to be at least as in Mexico City (12.63 × 10−3). ![Figure 4](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F4.medium.gif) [Figure 4](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F4) Figure 4 Optimal interventions for the COVID-19 pandemic Day 0 corresponds to the date of the first confirmed case in each city. Figure shows the optimal interventions for the SI model using two different proportions of symptomatic individuals (solid and dashed), assuming the median *R* in each city and considering *γ* = 1/5. Black lines denotes the start of the actual interventions in each city. ## DISCUSSION AND CONCLUSIONS Obtaining our criterion to design optimal NPIs was possible only because we characterized the necessary and sufficient conditions for the existence of solutions to an optimal control problem. Deriving such complete characterization is challenging —especially for high-dimensional systems— because it involves solving an infinite-dimensional optimization problem11. Indeed, computational methods12 cannot produce such characterization, and established analytical methods like Pontryagin’s Maximum Principle only yields necessary conditions for optimality11. We note there are several studies applying these and other similar methods to the SIR model13,14, in particular during the COVID-19 pandemic10,15–17. Our choice of a simple SIR model was motivated by its epidemiological adequacy for the COVID-19 pandemic and its low dimensionality, allowing us to apply Green’s Theorem to compare the cost of any two interventions analytically. In this sense, the method we use is closer to our previous work on optimal control for bioreactors18. Our results could guide a similar complete characterization for more detailed epidemic models or more detailed optimization objectives, but this is likely very challenging. Our choice of the simple SIR model also gives us a complete understanding of the optimal intervention at any state that the epidemic can take. The feedback form *u**(*S*, *I*) of the optimal intervention reflects such understanding, telling us the optimal action to perform if the epidemic is in the state (*S*, *I*). Contrast this feedback strategy to most other studies applying optimal control to epidemic diseases, where the optimal intervention is written in open-loop *u**(*t*) (see Supplementary Note S4 for details about how our work is related to existing optimal control studies). For a particular initial state, the open-loop intervention gives the optimal action at any time. However, it does not tell us what is the optimal action if we would have started at a different initial state. 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, it is well-established that feedback gives the control strategies the required robustness to work on real processes19. Our feedback optimal intervention has also this robustness (Supplementary Note S5). More precisely, we show that the optimal intervention works despite the presence of an incubation period, or when a fraction of the infected individuals remain hidden to the epidemic surveillance (e.g., because they are asymptomatic, as in COVID-19). Future work could also analyze the robustness of the optimal intervention when the current state of the epidemic is not entirely known. For example, this case may happen when significative delays exist in reporting newly infected cases, or when test for identifying infected individuals are limited. The optimal intervention resulting from our analysis takes 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ū**/*umax* days of maximum reduction with (*T* − *d*) days without intervention. This approach yields an intervention similar to Karin et al.20, with the difference that the periods of intervention and activity are optimally balanced. We will inevitably face new epidemics where non-pharmaceutical interventions are the only option to control infections. Our analysis suggests that, during the COVID-19 pandemic, less than half of the cities we studied attained enough health resources to respond adequately. Indeed, in some cities like Lima and Dehli, their health services’ capacity is so low that they are practically condemned to either exceed their capacity or live with near-permanent interventions. We must develop a scientific understanding that can inform the design of non-pharmaceutical interventions and plan the required health services capacity. We hope our work helps to catalyze the efforts to bring that understanding into practice. ## Data Availability Data is included in Supplementary Note 3 ## Funding M.T.A. gratefully acknowledges the financial support from CONACyT project A1-S-13909, México. J. X.V.H. acknowledges the financial support from UNAM IN115720. ## Correspondence Correspondence and requests for materials should be addressed to M.T.A. (email: mangulo@im.unam.mx). #### BOX 1. Optimal NPIs for the Susceptible-Infected-Removed (SIR) model The SIR model is a keystone in our understanding of infectious diseases, capturing the most essential features of the epidemiological dynamics for the mitigation or eradication of epidemic outbreaks2,21. The SIR model with interventions *u*(*t*) ∈ [0, *umax*] reducing disease transmission takes the form ![Formula][22] Here, *S*(*t*), *I*(*t*), and *R*(*t*) are the proportion of the population that is susceptible, infected, or removed (i.e., recovered or dead) at time *t* ≥ 0, respectively. Because *S*(*t*) + *I*(*t*) + *R*(*t*) = 1 for all *t* ≥ 0, we can just consider the (*S*, *I*) dynamics. We denote by (*S*, *I*) the initial state of the model at *t* = 0. The parameters are the *contact rate β* ≥ 0 and the the mean *residence time* of infected individuals *γ* ≥ 0 (in units of day−1). By assuming *S*(0) ≈ 1, these two parameters yield the *basic reproduction number R* = *β*/*γ*. In Supplementary Note S1, we prove that the optimal intervention is fully characterized by the *separating curve* of Eq. (1). This separating curve characterizes the optimal intervention as follows: (1) an optimal intervention exists if and only if the initial state (*S*, *I*) lies below this separating curve (i.e., *I* ≤ Ф(*S*)); (2) if it exists, the optimal intervention *u** takes the feedback form ![Formula][23] In words, the optimal intervention starts when *I*(*t*) reaches Ф(*S*(*t*)), and then it slides *I*(*t*) along Ф(*S*) until reaching the region where ![Graphic][24]. ## Supplementary Notes **Contents** * S1. Characterization of the optimal intervention in the Susceptible-Infected-Removed model 2 * S1.1 Calculation of the orbits 3 * S1.2 The number of infected people 4 * S1.3 Reachable set from (*S*, *I*) 4 * S1.4 Comparing the cost of two different trajectories 5 * S1.5 Determining the optimal orbit 6 * S1.5.1 Trivial 7 * S1.5.2 Unfeasible 7 * S1.5.3 The uncontrolled orbit is optimal 7 * S1.5.4 Singular arc 7 * S1.6 The optimal path 9 * S1.7 A feedback control strategy 11 * S2. Necessary and sufficient conditions for the existence of optimal NPIs 15 * S2.1 Necessary and sufficient conditions for existence 15 * S2.2 Calculating the duration of the optimal intervention 15 * S3. Application to the COVID-19 pandemic 16 * S3.1 Estimate for the fraction of infected individuals requiring intensive care 16 * S3.2 Data used in our analysis 16 * S4. Related work 18 * S5. Robustness of the optimal intervention 19 * S5.1 Robustness to the presence of an incubation period 19 * S5.2 Robustness to the presence of hidden infected individuals 19 ## S1. Characterization of the optimal intervention in the Susceptible-Infected-Removed model The model is given by ![Formula][25] 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*), which is what we will do. The maximal (acceptable) value of *I* is *Imax* and the maximal achievable value of the control is *umax*. So the state has to belong to the following feasible sets ![Formula][26] Sometimes it will be useful to write the differential equation in a compact form as ![Formula][27] The trajectory starting at the initial point *x* = (*S*, *I*) and subject to the control *u*: ![Graphic][28] is denoted by *ϕ* (*t*, *x*, *u* (·)). The optimal control problem consists in finding the control strategy *u* such that, starting from the initial point (*S*, *I*) the target set1 ![Formula][29] is reached in the minimal time with the state restriction *I*(*t*) ≤ *Imax* satisfied for all time. Now let us define the reachable set for an initial state x0 as the set of points that can be reached from the initial point *x* with feasible control, i.e., ![Formula][30] Also, we define the controllable set of the target set ![Graphic][31] as the set of points from which some point in the target ![Graphic][32] can be reached with a feasible control, i.e., ![Formula][33] The set ![Graphic][34] can be equivalently described as ![Graphic][35] for the system ![Formula][36] i.e., the set of points that can be reached from the set ![Graphic][37] for the dynamics with backward time. Now, the optimal control problem has a solution if and only if ![Formula][38] Since the points of the form (*S*, *I*) = (*S*, 0) are equilibria for every control value, *R* ((*S*, 0)) = (*S*, 0), so 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 *Ṡ* < 0 for *S* > 0, *I* > 0, ![Formula][39] 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 *XF*, i.e., ![Formula][40] ### S1.1 Calculation of the orbits Although it does not seem to be possible to find explicitly the trajectories of the system, it is easy to find its orbits. For this we write (we exclude the points for which *I* = 0 since they are equilibria) ![Formula][41] which is a separable differential equation (DE). Assuming that *u* is constant and integrating we obtain ![Formula][42] An interesting rewriting of (S1) is ![Formula][43] This means that the quantity ![Graphic][44] remains constant along the trajectory. Note that this constant depends on the control value used. The above equation is well-known for the SIR model21. Given an initial condition (*S*, *I*) this expression gives, for any 0 < *S* < *S* the (unique) value of *I* that is reached in future time.2 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][45] and, if we take the expression *I* (*S*; (*S*, *I*)), we obtain a separable DE that can be integrated, ![Formula][46] 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* ≤ *umax* 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][47] or, equivalently, as ![Formula][48] 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* = *I* and *I* ≈ 0, then ![Formula][49] Note that, if *u* → 1−, then *S* (∞) → 1− so, the larger the value of *u*, the larger the value of *S* (∞) also is. ### 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 *Fu* (*x*) = *f* (*x*) + *g* (*x*) (1 − *u*). The extreme values are given by *F* (*x*) = *f* (*x*) + *g* (*x*) and ![Graphic][50], ![Formula][51] 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][52] it follows that *F* is “above” ![Graphic][53]. Therefore,the reachable set *R* (*x*) is bounded by the two trajectories *ϕ* (*t*, *x*, *u* = 0) and *ϕ* (*t*, *x*, *u*max) (see Fig. S1). These two bounding orbits can be easily calculated using (S1). In particular, we can calculate the maximal value achieved by *I* for every (constant) control action, ![Formula][54] which is achieved at ![Formula][55] ![Supplementary Figure S1](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F5.medium.gif) [Supplementary Figure S1](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F5) Supplementary Figure S1 Phase Plane with maximal and minimal orbits bounding the reaching set *R* (*x*). Max corresponds to the trajectory *ϕ* (*t*, *x*, *u* = 0) while Min to *ϕ* (*t*, *x*, *u* = *umax*) and gives the maximum incidence as a function of the control strategy *u*: ![Formula][56] ### S1.4 Comparing the cost of two different trajectories In order to be able to find the optimal 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*, *xf*, u), *i =* 1, 2, joining the (same) points *x* and *xf* using two different control actions, *u*1 and *u*2, respectively. The cost going through *ωi* is (recall that our cost is time) ![Formula][57] 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][58] 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][59] Now, by properties of the determinant this is also the same as ![Formula][60] Therefore, ![Formula][61] 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][62] where ![Graphic][63] 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][64] where ![Graphic][65] is the region enclosed by the closed curve ![Graphic][66]. For our problem this becomes ![Formula][67] In our case, ![Formula][68] 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 Determining the optimal orbit 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 *Imax* 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. #### 51.5.1 Trivial This is the case in which ![Graphic][69]. That is, the case in which we start in the target set. #### 51.5.2 Unfeasible This is the case if *I* > *Imax*. #### 51.5.3 The uncontrolled orbit is optimal This will be the case if, and only if, the maximum *Ī* given by (S3) with u = 0 is such that 1. ![Graphic][70], or 2. ![Graphic][71] and ![Graphic][72] as given by (S2) with *u* = 0 satisfies ![Graphic][73] and ![Graphic][74], that is, 1 ≤ *S* and ![Graphic][75]. #### 51.5.4 Singular arc If *Ī* given by (S3) with *u* = o is larger than *I*max and ![Graphic][76], it is necessary to apply some control to maintain *I* below the maximal value *I*max. Now we calculate the value of *S* = *Sc* at which the orbit (first) touches *I*max. For this we solve (use (S1)) ![Formula][77] for *S* and obtain two solutions: *S*1, *S*2. Define *Sc* = max {*S*1, *S*2} as the largest. Now we calculate the value of *S* = *S** at which it is possible to achieve *İ* ≤ o (that is, that it is possible to stop the growth of *I*). This value can be calculated from ![Formula][78] and gives ![Formula][79] There are two possible situations: 1. Bang plus singular arc. If ![Graphic][80], then the optimal control strategy consists in ![Formula][81] where *using* is the control required to maintain *I* = *Imax* constant, i.e., *İ* = 0, ![Formula][82] 2. Bang-Bang plus singular arc. When *Sc* > *S** then it is necessary to start (if possible) with the control strategy before reaching the maximal value of *I* = *Imax*. Otherwise, this limit will be surpassed. However, this is only feasible if, moving backwards from the point (*S**, *Imax*) with the maximal control *umax* it is possible to reach a point ![Graphic][83] such that ![Graphic][84]. If ![Graphic][85], then it is not possible to solve the optimal problem, since any strategy will either surpass the maximal value *umax* or it will not reach the target. The value of ![Graphic][86] can be calculated from (S1), ![Formula][87] so that ![Formula][88] If ![Graphic][89], the optimal control is ![Formula][90] When ![Graphic][91], the control is given by ![Formula][92] where the value of (*Ss*, *Is*) is a switching point. It is characterized as follows: the trajectory starting at (*S*, *I*), i.e., *ϕ* (*t*, (*S*, *I*), *u* = 0) intersects the trajectory that starts at (*S**, *Imax*) but goes backwards in time, i.e., *ϕ* (−*t*, (*S**, *Imax*), *umax*). Such point (*Ss*, *Is*) can be calculated from (S1) as ![Formula][93] Substituting the first into the second we get ![Formula][94] Solving for *Ss* in the second we arrive at ![Formula][95] ![Supplementary Figure S2](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F6.medium.gif) [Supplementary Figure S2](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F6) Supplementary Figure S2 Phase Plane with the four key trajectories and points ## S1.6 The optimal path For obtaining the optimal path it is useful to construct some system trajectories (not all of them are really required for the final calculation, though). Let us define the initial point *x* = (*S*, *I*) and the final point as ![Graphic][96].*xf* is the point of the target set at the upper right corner. The trajectories (or orbits) we want to find are four: 1. *ϕ* (t, *x*, *u* = 0). In words, it is the trajectory without control starting at *x*. 2. *ϕ* (*t*, *x*, *umax*). In words, it is the trajectory with maximal control starting at *x*. 3. *ϕ* (−*t*, *xf, u =* 0. In words, it is the trajectory without control that ends in *xf*. 4. *ϕ* (−*t*, *xf*, *u**): In words, it is the trajectory with control ![Formula][97] that ends at *xf*. The control *u** is such that this trajectory does not violate the restriction *I* ≤ *I*max. For values of *S* ≤ *S**, it is equal to *umax*, and for *S* ≤ *S** it is the control for the singular arc, i.e., it maintains *I* = *I*max until *xf* is reached. These trajectories are presented in Fig. S2 for the parameters: *β* = 0.52, ![Graphic][98], *I*max = 0.1, *S* = 0.99, *I* = 0.01 and *umax =* 0.4. There are also some key points to find, besides the initial *x* and final *xf* ones. These are: * The “blue” point ![Graphic][99] in Fig. S2 where *ϕ* (*t*, *x*, *u* = 0) attains its maximum value, ![Formula][100] This is the same point at which this trajectory crosses the critical value ![Graphic][101], after which *İ <* 0 with zero control. When *Ī* ≤ *I*max then the control strategy is simply “do nothing” all the time, i.e. ![Formula][102] * The “green” point (*S**, *I*max). *S** corresponds to the maximal value of *S* for which it is possible to keep *I*max constant, given by ![Formula][103] We “saturate” the value of *S**, since *S** > 1 does not make sense. The control required to achieve this is the “singular” control ![Formula][104] Note that, if *S* > *S**, it is not possible to keep I at *I*max and *İ >* 0. If *S** = 1, then the optimal control is ![Formula][105] * The “black” point (*S*max, *I*max) is the point at which the trajectory *ϕ* (*t*, *x*, *u* = 0) reaches the value *I*max and can be calculated from the equation ![Formula][106] This equation has a unique solution in the interval ![Graphic][107] only if ![Graphic][108]. If *S** > *S*max, then the optimal control is (as in the previous case) ![Formula][109] * (see Fig. S3). * The “red” point (*Ss*, *Is*). At this point the trajectory *ϕ* (−*t*, *xf*, *u**) crosses the trajectory *ϕ* (*t*, *x*, *u* = 0). Note that if this crossing does not exist, then the optimal control problem is unfeasible. We have three possible cases: * – If *S** ≥ *S*max this point exists and it coincides with (*S*max, *I*max), i.e. (*Ss*, *Is*) = (*S*max, *I*max). The optimal control is as in the previous case. * – If *S** < *S*max the point (*Ss*, *Is*) may exist or not. We calculate ![Formula][110] * If *S** < *S** < *S* then (*Ss*, *Is*) exists and (*Ss*, *Is*) is given by ![Formula][111] The optimal control is a bang-bang-singular arc-bang strategy: ![Formula][112] This is the case of Fig. S2. * If *S* < *S** then (*Ss*, *Is*) does not exist. We fix its value (arbitrarily) as (*Ss*, *Is*) = (1, 0). In this case the optimal control problem is unfeasible. This is the case if in our example we set *u*max ≤ 0.35. See Fig. S4. ![Supplementary Figure S3](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F7.medium.gif) [Supplementary Figure S3](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F7) Supplementary Figure S3 Phase Plane for *u*max = 0.57 ## S1.7 A feedback control strategy The previous “open loop” strategy can be implemented as a state feedback control. This strategy is rather simple, since there is basically only one switching curve: *ϕ* (−*t*, *xf*, *u**) (this is the discontinuous blue line in the Figures). There is a second switch: when the target region has been attained, the control is switched off, but this happens in a “natural” manner. The switching curve is defined as ![Formula][113] where ![Formula][114] ![Supplementary Figure S4](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F8.medium.gif) [Supplementary Figure S4](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F8) Supplementary Figure S4 Phase Plane with *u*max = 0.35 The optimal control feedback is thus given by ![Formula][115] Alternatively, we can implement a pure switching control since the “equivalent control”22 will realize the singular control on the singular arc, ![Formula][116] Note that this control strategy extends the control action beyond the region where the optimal control is feasible. This extension is not strictly based on the value function, and therefore there is not a unique way to do so. In our case, for example, the zero control region is extended to the (non feasible) region defined by ![Graphic][117]. But another possible (and maybe better) extension is to set *u* = *u*max in this region, since then the limit *I*max will be reached faster than without control action. The resulting controller is then given by ![Formula][118] with ![Formula][119] where we have slightly changed the switching function to include this region. Some results can be seen in Figs. S5 and S6.6 From these figures one also observes that, an extra benefit of applying some control compared to not doing anything is that, when the infection dies, the total number of infected people is larger if no action is taken than if some control action is applied. We see this in the figures by noting that *S* (∞) is larger with control than without it. This number can be further increased if instead of taking no control once ![Graphic][120] one still applies some control action (of course, the number of infected people is minimized by applying *u*max). ![Supplementary Figure S5](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F9.medium.gif) [Supplementary Figure S5](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F9) Supplementary Figure S5 Optimal feedback control with *umax* = 0.25, so that the problem is unfeasible. ![Supplementary Figure S6](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F10.medium.gif) [Supplementary Figure S6](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F10) Supplementary Figure S6 Optimal feedback control with *umax* = 0.6, so that the problem is feasible. ## S2. Necessary and sufficient conditions for the existence of optimal NPIs ### S2.1 Necessary and sufficient conditions for existence 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 *I* ≤ Ф(*S*) where Ф(*S*) is the separating curve. To characterize a condition that is independent of the initial state, we consider the limit case of *S* = l and *I* = 0. Under this assumption, the necessary and sufficient condition of existence is that Ф(1) >0. 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 of Eq. (1), we obtain the condition ![Formula][121] which is precisely the inequality (2). ### 52.2 Calculating the duration of the optimal intervention In the SI model ![Formula][122] the duration of the optimal intervention depends on the parameters *β*, *γ*, *u*max and *I*max To describe duration only in terms of *R* = *β*/*γ*, *u*max and *I*max, we note the above system can be rewritten as ![Formula][123] where τ = *γt*,. In this last rewriting, time is in units of mean residence time of infected individuals. We simulated the above normalized SIR model on a time interval *t* ∈ [0, *tf*] with *u*(*t*) = *u**(*t*) = *u**(*S*(*t*), *I*(*t*)) given by the optimal feedback intervention of Eq. (S4). Then, the duration i of the optimal intervention is given by the length of time where *u**(*t*) > 0, i.e., ![Formula][124] ## S3. Application to the COVID-19 pandemic ### S3.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 variability23, ranging from a 20/100 in a report of the World Health Organization, to 96/100 in a study of 328 adults in Shanghai24. We take the nominal value of *p* = 60/100. For the results of Fig. 3 we consider the interval of uncertainty *p* ∈ [35/100, 50/100]. 2. Following Kremer et al.25, we assume that from the individuals that are symptomatic, a fraction 15/100 develop severe symptoms. 3. Finally, following Li et al.26, 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][125] ## S3.2 Data used in our analysis Supplementary Fig. S7 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. S7. 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 *Rt* 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/). 5. **Start of NPIs for each city**. We used data from the Oxford Coronavirus Government Response Tracker27. For our analysis, we considered only the start dates for school closing and workplace closing. Here we assumed that the NPIs started in each city at the same time they started in the country. The time of start of the NPIs was calculated with respect to the date of the first confirmed case in each country. ![Supplementary Figure S7](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F11.medium.gif) [Supplementary Figure S7](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F11) Supplementary Figure S7 Table with the response of 16 cities during the COVID-19 pandemic. ## S4. 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 quarantines28,29, drug treatments30, or dispersal of insecticides and education campaigns31. The standard tool to solve these optimal control problem is the celebrated Pontryagin’s Maximum Principle32. 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 activity20. Optimal control methods have been also applied, for example to minimize the peak of infection33, minimize the number of infections17, minimize the economic costs15, or maximize welfare16. 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. ## S5. Robustness of the optimal intervention Here we evaluate the robustness of the optimal intervention agains dynamics that were not considered in its derivation. ### S5.1 Robustness to the presence of 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][126] 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. Supplementary Fig. S8 shows the result of applying the optimal intervention to the above SEIR with different values of the incubation period. Namely, we apply *u*(*t*) = *u**(*I*(*t*), *S*(*t*)) with *u**(*I*, *S*) as in Eq. (S4). When the incubation period is small compared to the period of the disease (i.e., 1 /*γ*), the response of the optimal intervention in the original SIR model and in the SEIR are very similar. These incubation periods are reasonable for influenza-like diseases. We do find that prevalence exceeds *I*max, but this excess is less than 1% (Supplementary Fig. S8a). This confirms the robustness of the optimal intervention to the presence of incubation period. The response of SEIR with the optimal intervention remains acceptable up to an incubation period that is about 75% of the disease period (Supplementary Fig. S8b-c). For this last case, the incubation period is similar to what has been observed for the SARS-CoV-2 virus. ### S5.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][127] Above, *Is* denotes the proportion of symptomatic infections and *Ia* the fraction of asymptomatic ones. The model assumes that a fraction *p* ∈ [0,1] of infections result in symptomatic cases, and the rest (1 − *p*) in asymptomatic ones. We assume that infectious period 1/*γ* is the same for both symptomatic and asymptomatic individuals. Since we assume that only symptomatic individuals end up requiring hospital care, we consider that the objective is to keep *Is*(*t*) ≤ *I*max only. The control applied is *u*(*t*) = *u**(*Is*(*t*), *S*(*t*)) where *u**(*I*, *S*) is given by Eq. (S4). Supplementary Fig. S9 shows the result of applying the optimal intervention to the above model with different values of the proportion of symptomatic cases *p*. For the wide range *p* ∈ [55/100,80/100], we find that effect of hidden infections is that the optimal intervention becomes over cautious, in the sense that now *Is* never reaches *I*max. ![Supplementary Figure S8](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F12.medium.gif) [Supplementary Figure S8](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F12) Supplementary Figure S8 Robustness of the optimal intervention to the cpresence of an incubation period Parameters of the model not specified in the panels are: *β* = 0.52, *u*max = 0.6, *I*max = 1/8, *I* = 1.129 × 10−7, *S* = 1 − *I*, and *E* = 0. Solid line denotes the proportion of infected individuals under the optimal intervention. Dashed line denotes the proportion of infected individuals without intervention. **a**. Small incubation period. **b**. Medium incubation period. **c**. Large incubation period. ![Supplementary Figure S9](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2020/05/21/2020.05.19.20107268/F13.medium.gif) [Supplementary Figure S9](http://medrxiv.org/content/early/2020/05/21/2020.05.19.20107268/F13) Supplementary Figure S9 Robustness of the optimal intervention to the presence of hidden infections Parameters of the model not specified in the panels are: *β* = 0.52, *γ* = 1/7, *u*max = 0.6, *u*max = 1/8, *IS0* = *Ia0* = 1.129 × 10−7, *S* = 1 − 2*Is*. Pinkdenotes the proportion of symptomatic infected individuals *Is*(*t*). The grey area denotes the proportion of asymptomatic infected individuals *Ia*(*t*). **a**. Large proportion of symptomatic individuals *p*. **b**. Medium proportion of symptomatic individuals *p*. **c**. Small proportion of symptomatic individuals *p*. ## Acknowledgements We thank Luis Rodrigo Moreno Morton, Pablo Barberis, Leonid Fridman, Saul Hernandez, Alejandro Vargas, Sebastian Michel Mata for helpful discussions. ## Footnotes * 1 Note that this set is positively invariant without control (*u* = 0), and when the trajectory is in this set it will die out without control action. * 2 If we select *S* > *S* the obtained value of *I* is reached in a past time (*t* < 0). * Received May 19, 2020. * Revision received May 19, 2020. * Accepted May 21, 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]. R. M. Anderson and R. M. May. Vaccination and herd immunity to infectious diseases. Nature 318 no. 6044, pp. 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%2F05%2F21%2F2020.05.19.20107268.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1985AVC5500047&link_type=ISI) 2. [2]. R. M. Anderson and R. M. May. Infectious diseases of humans: dynamics and control (Oxford university press, 1992). 3. [3]. C. Fraser et al. Factors that make an infectious disease outbreak controllable. Proceedings of the National Academy of Sciences 101 no. 16, pp. 6146–6151 (2004). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMToiMTAxLzE2LzYxNDYiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMC8wNS8yMS8yMDIwLjA1LjE5LjIwMTA3MjY4LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 4. [4]. K. Kupferschmidt. Ending coronavirus lockdowns will be a dangerous process of trial and error. Science| AAAS. 5. [5]. N. M. Ferguson et al. Strategies for mitigating an influenza pandemic. Nature 442 no. 7101, pp. 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%2F05%2F21%2F2020.05.19.20107268.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000239278900042&link_type=ISI) 6. [6]. T. D. Hollingsworth et al. Mitigation strategies for pandemic influenza a: balancing conflicting policy objectives. PLoS computational biology 7 no. 2. 7. [7]. N. Ferguson 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]. R. M. Anderson et al. How will country-based mitigation measures influence the course of the covid-19 epidemic?. The Lancet 395 no. 10228, pp. 931–934 (2020). 9. [9]. J. Elston et al. The health impact of the 2014–15 ebola outbreak. Public Health 143, pp. 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%2F05%2F21%2F2020.05.19.20107268.atom) 10. [10]. R. Djidjou-Demasse et al. Optimal covid-19 epidemic control until vaccine deployment. *medRxiv*. 11. [11]. S. Lenhart and J. T. Workman. Optimal control applied to biological models (CRC press, 2007). 12. [12]. F. Bonnans, P. Martinon and V. Grélard. Bocop-a collection of examples. 13. [13]. E. Hansen and T. Day. Optimal control of epidemics with limited resources. Journal of mathematical biology 62 no. 3, pp. 423–451 (2011). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20407775&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2020%2F05%2F21%2F2020.05.19.20107268.atom) 14. [14]. P. Di Giamberardino and D. Iacoviello. Optimal control of sir epidemic model with state dependent switching cost index. Biomedical signal processing and control 31, pp. 377–380 (2017). 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. 17. [17]. C. Tsay et al. Modeling, state estimation, and optimal control for the us covid-19 outbreak. *arXiv preprint arXiv:2004.06291*. 18. [18]. J. Moreno. Optimal time control of bioreactors for the wastewater treatment. Optimal Control Applications and Methods 20 no. 3, pp. 145–164 (1999). 19. [19]. K. J. Åström and R. M. Murray. Feedback systems: an introduction for scientists and engineers (Princeton university press, 2010). 20. [20]. O. Karin et al. Adaptive cyclic exit strategies from lockdown to suppress covid-19 and allow economic activity. *medRxiv*. 21. [21]. H. R. Thieme. Mathematics in population biology, vol. 12 (Princeton University Press, 2018). 22. [22]. 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). 23. [23].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/). 24. [24]. X. Zhou et al. Follow-up of the asymptomatic patients with sars-cov-2 infection. Clinical Microbiology and Infection. 25. [25]. 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/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzNjgvNjQ5MC80OTMiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMC8wNS8yMS8yMDIwLjA1LjE5LjIwMTA3MjY4LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 26. [26]. R. Li et al. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (sars-cov2). Science. 27. [27]. T. P. S. W. B. K. Thomas Hale, Anna Petherick and N. Angrist. Oxford coronavirus government response tracker. URL [https://www.bsg.ox.ac.uk/research/research-projects/coronavirus-government-response-trackerl](https://www.bsg.ox.ac.uk/research/research-projects/coronavirus-government-response-trackerl). 28. [28]. H. Behncke. Optimal control of deterministic epidemics. Optimal control applications and methods 21 no. 6, pp. 269–285 (2000). 29. [29]. 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). 30. [30]. 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). 31. [31]. 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). 32. [32]. D. Liberzon. Calculus of variations and optimal control theory: a concise introduction (Princeton University Press, 2011). 33. [33]. D. H. Morris et al. Optimal, near-optimal, and robust epidemic control. *arXiv preprint arXiv:2004.02209*. [1]: /embed/graphic-2.gif [2]: /embed/inline-graphic-1.gif [3]: F2/embed/inline-graphic-2.gif [4]: F2/embed/inline-graphic-3.gif [5]: /embed/inline-graphic-4.gif [6]: /embed/graphic-4.gif [7]: /embed/inline-graphic-5.gif [8]: /embed/inline-graphic-6.gif [9]: /embed/inline-graphic-7.gif [10]: /embed/inline-graphic-8.gif [11]: /embed/inline-graphic-9.gif [12]: /embed/inline-graphic-10.gif [13]: /embed/inline-graphic-11.gif [14]: /embed/inline-graphic-12.gif [15]: /embed/inline-graphic-13.gif [16]: /embed/inline-graphic-14.gif [17]: /embed/inline-graphic-15.gif [18]: /embed/inline-graphic-16.gif [19]: /embed/inline-graphic-17.gif [20]: /embed/inline-graphic-18.gif [21]: /embed/inline-graphic-19.gif [22]: /embed/graphic-7.gif [23]: /embed/graphic-8.gif [24]: /embed/inline-graphic-20.gif [25]: /embed/graphic-9.gif [26]: /embed/graphic-10.gif [27]: /embed/graphic-11.gif [28]: /embed/inline-graphic-21.gif [29]: /embed/graphic-12.gif [30]: /embed/graphic-13.gif [31]: /embed/inline-graphic-22.gif [32]: /embed/inline-graphic-23.gif [33]: /embed/graphic-14.gif [34]: /embed/inline-graphic-24.gif [35]: /embed/inline-graphic-25.gif [36]: /embed/graphic-15.gif [37]: /embed/inline-graphic-26.gif [38]: /embed/graphic-16.gif [39]: /embed/graphic-17.gif [40]: /embed/graphic-18.gif [41]: /embed/graphic-19.gif [42]: /embed/graphic-20.gif [43]: /embed/graphic-21.gif [44]: /embed/inline-graphic-27.gif [45]: /embed/graphic-22.gif [46]: /embed/graphic-23.gif [47]: /embed/graphic-24.gif [48]: /embed/graphic-25.gif [49]: /embed/graphic-26.gif [50]: /embed/inline-graphic-28.gif [51]: /embed/graphic-27.gif [52]: /embed/graphic-28.gif [53]: /embed/inline-graphic-29.gif [54]: /embed/graphic-29.gif [55]: /embed/graphic-30.gif [56]: /embed/graphic-32.gif [57]: /embed/graphic-33.gif [58]: /embed/graphic-34.gif [59]: /embed/graphic-35.gif [60]: /embed/graphic-36.gif [61]: /embed/graphic-37.gif [62]: /embed/graphic-38.gif [63]: /embed/inline-graphic-30.gif [64]: /embed/graphic-39.gif [65]: /embed/inline-graphic-31.gif [66]: /embed/inline-graphic-32.gif [67]: /embed/graphic-40.gif [68]: /embed/graphic-41.gif [69]: /embed/inline-graphic-33.gif [70]: /embed/inline-graphic-34.gif [71]: /embed/inline-graphic-35.gif [72]: /embed/inline-graphic-36.gif [73]: /embed/inline-graphic-37.gif [74]: /embed/inline-graphic-38.gif [75]: /embed/inline-graphic-39.gif [76]: /embed/inline-graphic-40.gif [77]: /embed/graphic-42.gif [78]: /embed/graphic-43.gif [79]: /embed/graphic-44.gif [80]: /embed/inline-graphic-41.gif [81]: /embed/graphic-45.gif [82]: /embed/graphic-46.gif [83]: /embed/inline-graphic-42.gif [84]: /embed/inline-graphic-43.gif [85]: /embed/inline-graphic-44.gif [86]: /embed/inline-graphic-45.gif [87]: /embed/graphic-47.gif [88]: /embed/graphic-48.gif [89]: /embed/inline-graphic-46.gif [90]: /embed/graphic-49.gif [91]: /embed/inline-graphic-47.gif [92]: /embed/graphic-50.gif [93]: /embed/graphic-51.gif [94]: /embed/graphic-52.gif [95]: /embed/graphic-53.gif [96]: /embed/inline-graphic-48.gif [97]: /embed/graphic-55.gif [98]: /embed/inline-graphic-49.gif [99]: /embed/inline-graphic-50.gif [100]: /embed/graphic-56.gif [101]: /embed/inline-graphic-51.gif [102]: /embed/graphic-57.gif [103]: /embed/graphic-58.gif [104]: /embed/graphic-59.gif [105]: /embed/graphic-60.gif [106]: /embed/graphic-61.gif [107]: /embed/inline-graphic-52.gif [108]: /embed/inline-graphic-53.gif [109]: /embed/graphic-62.gif [110]: /embed/graphic-63.gif [111]: /embed/graphic-64.gif [112]: /embed/graphic-65.gif [113]: /embed/graphic-67.gif [114]: /embed/graphic-68.gif [115]: /embed/graphic-70.gif [116]: /embed/graphic-71.gif [117]: /embed/inline-graphic-54.gif [118]: /embed/graphic-72.gif [119]: /embed/graphic-73.gif [120]: /embed/inline-graphic-55.gif [121]: /embed/graphic-76.gif [122]: /embed/graphic-77.gif [123]: /embed/graphic-78.gif [124]: /embed/graphic-79.gif [125]: /embed/graphic-80.gif [126]: /embed/graphic-82.gif [127]: /embed/graphic-83.gif