key: cord-0107900-e04w4ota
authors: Ma, Qianqian; Liu, Yang-Yu; Olshevsky, Alex
title: Optimal Lockdown for Pandemic Control
date: 2020-10-24
journal: nan
DOI: nan
sha: b9d943b8613c1ad34de28cae0794f185d67dd3d1
doc_id: 107900
cord_uid: e04w4ota

As a common strategy of contagious disease containment, lockdowns will inevitably weaken the economy. The ongoing COVID-19 pandemic underscores the trade-off arising from public health and economic cost. An optimal lockdown policy to resolve this trade-off is highly desired. Here we propose a mathematical framework of pandemic control through an optimal stabilizing non-uniform lockdown, where our goal is to reduce the economic activity as little as possible while decreasing the number of infected individuals at a prescribed rate. This framework allows us to efficiently compute the optimal stabilizing lockdown policy for general epidemic spread models, including both the classical SIS/SIR/SEIR models and a new model of COVID-19 transmissions. We demonstrate the power of this framework by analyzing publicly available data of inter-county travel frequencies to analyze a model of COVID-19 spread in the 62 counties of New York State. We find that an optimal stabilizing lockdown based on epidemic status in April 2020 would have reduced economic activity more stringently outside of New York City compared to within it, even though the epidemic was much more prevalent in New York City at that point. Such a counterintuitive result highlights the intricacies of pandemic control and sheds light on future lockdown policy design.

1. Introduction. The COVID-19 pandemic has resulted in more than 92.3M confirmed cases and 2.0M deaths (up to Jan 13th, 2021) [44] and has impacted the lives of more than 90% global population [18, 79] . Curbing the spread of the pandemic like COVID-19 depends critically on the successful implementation of non-pharmaceutical interventions such as lockdowns, social distancing, shelter in place orders, contact tracing, isolation, and quarantine [20, 64, 27, 30] . However, these interventions can also lead to substantial economic damage, motivating us to investigate the problem of curbing pandemic spread while minimizing the induced economic losses.

We consider the problem of designing an optimal stabilizing lockdown that minimizes the economic damage while reducing the number of new infections to zero at a prescribed rate. Such a lockdown should be non-uniform, because shutting down different locations has different implications both for the economic cost and for pandemic spread. The difficulty is whereas a uniform lockdown can be found through a search over a single parameter, a non-uniform lockdown is parametrized by many parameters associated with different locations. Despite of its significance and implications, a computationally efficient framework to design stabilizing lockdown strategies is still lacking.

Here we propose such a framework by mapping the design of optimal lockdown policy to a classical problem in control theory -design an intervention that affects the eigenvalues of a matrix governing the dynamics of a dynamical system. It turns out that, even though general epidemic spreading dynamics are highly nonlinear, an eigenvalue bound for a linear approximation of the spreading dynamics nevertheless forces the number of infections at each location to go to zero asymptotically at a prescribed rate for all time. We provide two polynomial-time algorithms that design the optimal lockdown to achieve such an eigenvalue bound.

We apply these algorithms to design an optimal stabilizing lockdown on both synthetic and real data (using data from SafeGraph [74] to fit a county-level model of New York State) for epidemic spread models of COVID-19 using disease parameters from the literature [31, 7, 8] . Unsurprisingly, we find that the heterogeneous lockdown is far more economical than a homogeneous lockdown. However, we find additional features of the optimal stabilizing lockdown that are counter-intuitive. For example, we find that in models of random graphs, degree centrality and population do not affect the strength of the lockdown of a location unless its population (or degree centrality) takes extremely smaller (or larger) values than others. Most surprisingly, we show that an optimal stabilizing lockdown based on the epidemic status in April 2020 would

Here s i (x a i or x s i ) stands for the proportion of susceptible (asymptomatic or symptomatic infected, respectively) population at location i, a i j captures the rate at which infection flows from location j to location i, β a (or β s ) is the transmission rate of asymptomatic (or symptomatic) infected individuals, r a (or r s ) is the recovery rate of asymptomatic (or symptomatic) infected individuals. We assume infected individuals are asymptomatic at first and is the rate at which they develop symptoms. We use different parameters for symptomatic and asymptomatic individuals because a recent study [49] reported that asymptomatic individuals have viral load that drops more quickly, so they not only recover faster, but also are probably less contagious.

Note that our model of COVID-19 spreading can be considered as a generalization of the classical SIR model and the SEIR model of epidemic spread. Indeed, by setting β s = = r s = 0, we recover the SIR model; and by setting β a = r a = 0, we recover the SEIR model. However, neither the SIR nor the SEIR model captures the existence of two classes of individuals who transmit infections at different rates as above. Our model can also be considered as a simplification of existing models in studying COVID-19 spreading [8, 31] . For example, in [8] , asymptotic stability was considered in a slightly more general model including both births and deaths. Here for simplicity in our model we consider a fixed population size. In [31] eight classes of patients (instead of two) were introduced, depending on whether the infection is diagnosed, whether the patient is hospitalized, as well as other factors.

In matrix form, we can write our model as

It turns out that, if we want the number of infections at each location (or a linear combination of those numbers) to go to zero at a prescribed rate α, we just need to ensure that the linear eigenvalue condition λ(M(t 0 )) ≤ −α holds (see SI Sec. 4.3 for a formal proof). Note that this is quite different from what usually happens in nonlinear systems when we pass to an eigenvalue bound of at a point: here as long as the eigenvalue condition λ(M(t 0 )) ≤ −α is satisfied, we obtain that infections go to zero at rate α over all times t ≥ t 0 .

In the remainder of this paper, we will attempt to design strategies that enforce decay of infections with a prescribed rate by modifying the matrix A through lockdowns to satisfy such an eigenvalue bound. This is different from the more traditional approach of optimal control of network epidemic processes [9, 1, 26, 2] in a number of ways. First, this gives rise to a fixed lockdown, whereas a traditional optimal control approach would result in a lockdown that is different at every time t, which is obviously unrealistic. If the timevarying lockdown is approximated through a series of infrequently changing fixed lockdowns, the optimality guarantee is lost. Second, the optimal control approach results in lockdowns that relax in strength as the number of infections decreases. If policymakers are tasked with repeated lockdown relaxations, a potential danger is that political considerations will result in lockdowns that are too loose, leading cases to increase again. For example, a recent CDC report found that pre-mature relaxations of restrictions drove an increase of cases throughout the United States in 2020 [38] . Finally, as we will see later, one of the main benefits of our approach is the guaranteed scalability: the main result of this paper is a nearly linear time algorithm. By contrast, optimal control of epidemic processes is based on methods which are either known to be nonscalable or which sometimes fail to converge at all. We discuss this at more length in Section I of the Supplementary Information.

Lockdown model. Methods of constructing A capturing spatial heterogeneity have been well studied [36, 55, 75, 3, 5, 66, 16, 24] . Here we follow a recent work [8] that is particularly well-suited to model the lockdowns in curbing COVID-19. Denote the fixed population size at location i as N i , and assume people travel from location i to location j at rate τ i j . It is well accepted that such travel rates determine the evolution of an epidemic. For example, it has been reported that regional progress of influenza is much more correlated with the movement of people to and from their workplaces rather than geographic distances [81] ; in the context of COVID-19, mobility based on cell-phone data has been predictive as a measure of epidemic spread [19, 32] . The quantities a i j can then be determined as (see SI Sec. 4.2 for details):

It is intuitive that a i j is the sum of the terms involving τ il τ jl since this product captures the interactions between people from locations i and j through visits to location l. Eq. (2.3) can also be written in matrix form as A = CB with

where τ = (τ i j ), D 1 = diag( k N k τ kl ) −1 while D 2 = diag(N 1 , . . . , N n ). When a lockdown is ordered heterogeneously across different locations, this has two consequences. First, the transmission rates will be altered. For instance, ensuring that all buildings have a maximum enforced density limits the rate at which people can interact, as do mandatory face-covering, and other measures, resulting in a number of transmissions that is a fraction of what they otherwise would have been. We may account for this as follows. From Eq. (2.3), we have that β a a i j = n l=1 β a τ il τ jl N j n k=1 N k τ kl . 3 This manuscript is for review purposes only.

The effect of the lockdown is to replace β a in each term of the sum by β a f l , for some location-dependent f l ∈ [0, 1]. The effect on β s a i j is similar. Secondly, travel rates to location l are also a fraction of what they were before since there is now reduced inducement to travel, i.e., τ il should be replaced by τ il g l with g l ∈ [0, 1] for each location l. To avoid overloading the notation, we will not change the definitions of β a , β s or the travel rates τ il but instead achieve the same effect by changing the definition of a i j as:

where z l = f l g l ∈ [0, 1]. In matrix notation, the post-lockdown A matrix is

The quantities z 1 , . . . , z n can be thought of as measuring the intensity of the lockdown at each location.

Lockdown cost. Clearly, setting z l = 1 corresponds to doing nothing and should have a zero economic cost. On the other hand, choosing z l = 0 corresponds to a complete lockdown and should be avoided. We will later apply our framework to real data collected from counties in New York State; shutting down a county entirely would result in people being unable to obtain basic necessities, and thus the economic cost should approach +∞ as z l → 0. With these considerations in mind, a natural choice of lockdown cost is (2.6) c(z 1 , . . . , z n ) = n i=1 c l 1 z l − 1 .

Here c l captures the relative economic cost of closing down location l. Throughout this paper, we will choose c l to be the employment at local l, but other choices are also possible (e.g., c l could be the GDP generated at location l).

Besides the cost function in Eq. (2.6), we will also consider cost functions that blow up with different exponents as n i=1 c l z −k l − 1 , as well as costs which threshold as n i=1 min c l (z −1 l − 1), C l which alter our cost function by saturating at some location-dependent cost C l rather than blowing up as z l → 0.

One of the advantage of these cost functions is that they inherently discourage extreme disparities among nodes. Indeed, a lockdown that places all the burden on a small collection of nodes by setting their z i close to zero will have cost that blows up. By contrast, some previous works such as [8, 17] used cost functions certain locations contribute little to disease spread but have very high relative economic cost of lockdown, one could even increase activity in these locations to allow for a harsher lockdown elsewhere with the same overall cost. While our methods work for both variations, all of our simulations and empirical results will consider the constrained version. Our approach of stabilizing the system by forcing the eigenvalues to have negative real part is a standard heuristic in control theory [52, 62, 33] . This approach comes with a caveat-if pushed to the extreme by moving the eigenvalues further and further towards negative infinity, the asymptotically better performance will start coming at the expense of the non-asymptotic behavior of the system.

We next turn to describe our main results. Our first step is to discuss an assumption required by one of our algorithms, i.e., the recovery rate γ has to be small relative to the entries in the matrices C and B. We call this "high-spread assumption" (see Assumption 4.5 in SI Sec. 4.4 for the formal description). We will later show that, under this assumption, the constrained and unconstrained lockdown problems are equivalent. This is quite intuitive: if the epidemic spreads sufficiently fast everywhere, the unconstrained shutdown will never choose to increase the activity of any location.

Main theoretical contribution. With the above preliminaries in place, we can now state our main theoretical contribution. Our main theorem provides efficient algorithms for both the unconstrained and constrained lockdown problems (see Theorem 4.6 in SI Sec. 4.4 for the formal description of this theorem). In particular, we first prove that the unconstrained lockdown problem for both SIS and COVID-19 models can be exactly mapped to the classical matrix balancing problem (see SI Sec. 3.2 for details) and solved with nearly linear time complexity. Moreover, we prove that if the "high-spread" assumption holds, then the constrained lockdown problem is equivalent to the unconstrained lockdown problem and consequently is also reducible to matrix balancing. Even if the "high-spread" assumption does not hold, we prove that under certain conditions the constrained lockdown problem for the SIS and COVID-19 models can still be solved by applying the covering semi-definite program with polynomial time complexity ofÕ(n 3 ), where the tilde hides factors logarithmic in model parameters. To summarize, we give three separate cases that cover all possible scenarios. In two of these cases, the optimal stabilizing lockdown problem is solvable in nearly linear time, and in the remaining case, it is solvable inÕ(n 3 ).

In practice, we find the optimal lockdown problem is solvable in linear time in the vast majority of the cases. Specifically, when we fit the models to New York State data, in 23 experiments out of 27, the linear time algorithm gave the correct answer.

We now apply the algorithms we've developed to design an optimal stabilizing lockdown policy for the 62 counties in the State of New York (NY). Our goal is to reduce activity in each county in a non-uniform way to curb the spread of COVID-19 while simultaneously minimizing economic cost. The data sources we employed are presented in SI Sec. 6 .

We consider only the constrained lockdown problem here. When the "high-spread" assumption is satisfied, we will apply the matrix-balancing algorithm, otherwise we will apply the covering semi-definite program. To provide valid estimation results, we employed three different sets of disease parameters provided in literature [7] , [8] , [31] (See Supplementary Table 2) .

Comparison with other lockdown policies. We used the data of the 62 counties in NY on April 1st, 2020 as initialization and estimated the number of active cases over 300 ∼ 800 days and the number of cumulative cases over 500 ∼ 1, 500 days with different lockdown policies. Fig.2a -c show the estimated active cases over times, and Fig.2d -e show the estimated cumulative cases over time. Here results from different columns of Fig.2 were calculated by using different sets of disease parameters adopted from literature [8, 31, 7] .

We compared the optimal stabilizing lockdown policy calculated by our methods with several other benchmark policies: (1) no lockdown, i.e., z l = 1 for all locations; (2) random lockdown, z l is randomly chosen from a uniform distribution U[a, b] with the lower-and upper-bounds a and b chosen such that the overall cost of this policy is the same as that of our policy; (3) uniform lockdown, where z l = z is the same for all the locations and z is chosen such that the overall economic cost is the same as that of our policy; (4) uniformly-bounded-decline locdown, where the "decline" is uniformly bounded across locations, i.e., the decay rate of the infections in each location is bounded by a constant α, where α is chosen such that the economic cost of this policy is the same as that of our policy. Note that among the four benchmark policies both the uniformly-bounded-decline lockdown and the random lockdown are heterogeneous locationwise.

From Fig.2 , we can see that our optimal stabilizing lockdown policy outperforms all other lockdown polices in terms of the total final number of cumulative cases. Similar findings for the SIS and SIR models are reported in Fig.1 and Fig.2 .

Optimal stabilizing lockdown rate z * l for each county. Fig.3 shows lockdown-rate profiles z l calculated by various policies. First, we found that the optimal stabilizing lockdown profile is quite sensitive to the disease parameters. Second, surprisingly, the values of z * l for counties in NYC are relatively higher (corresponding to a less stringent lockdown) than that of counties outside NYC, regardless of the disease parameters. This is a counter-intuitive result: even though the epidemic was largely localized around NYC on the date we used to initialize the infection rates, the calculated optimal stabilizing lockdown profile indicates that it is cheaper to reduce the spread of COVID-19 by being harsher on neighboring regions with smaller populations. It can be observed from Fig. 3 that this pattern only appears in our heterogenous optimal stabilizing lockdown policy.

To see why this is counter-intuitive, consider the case of a single-node (i.e., non-network) model. It is easy to see that the optimal stabilizing lockdown is insensitive to population. Intuitively, doubling population doubles the cost of the lockdown and also doubles the benefits in terms of lives saved. In terms of our model, the stabilization constraint is on the proportion of infected, so doubling the population may change the optimal cost but does not change the optimal solution. Furthermore, the strength of the optimal stabilizing lockdown is increasing in s(t 0 ): harsher restrictions are needed to achieve the same decay rate if more people are infected. Thus it is surprising that when we consider a network model of New York State, the region with the highest population and highest proportion of infected is treated the lightest under the optimal stabilizing shutdown.

In SI Sec. 8, we further replicate the same finding in a much simpler city-suburb model: we consider a city with large population and a neighboring suburb with small population and observe that the optimal stabilizing lockdown will choose to shut down the suburb more stringently. In SI Sec. 11, we further confirmed the same counterintuitive phenomenon using other cost functions. In SI Sec. 12, we checked the robustness of this counterintuitive phenomenon with respect to the uncertainty of the travel rate matrix τ by adding noises or removing part of the travelling data. It turns out that this phenomenon is quite robust against the uncertainty of the matrix τ . In particular, perturbing each τ i j by noise with variance up to 10τ 2 i j , preserves the result, as does randomly setting half of the τ i j to zero (See Fig. 18 ).

We investigate this finding further in SI Sec. 9, where we let the metric to be optimized be the total number of infections. As a counterpoint, we do a greedy search over all two-parameter lockdowns that shut down NYC harder than the rest of New York State. Our results show that our lockdown has a smaller number of infections than the best two-parameter lockdown of the same cost.

This phenomenon likely occurs because shutting down a suburb with small population yields benefits proportional to the much larger population of the city, since a shutdown in the suburb affects the rate at which infection spreads in the city as city residents can infect each other through the suburb. Similarly, a possible explanation for the phenomenon we observe on New York State data is that shutdowns outside of 6 NYC may be a cheaper way to curb the spread of infection within NYC. We stress that this effect is due to the network interactions. In particular, this counterintuitive phenomenon does not occur in a hypothetical model of New York State where residents always stay within their own county: in that case, the lockdown problem reduces to a collection of single-node models which do not interact, and optimal stabilizing shutdown will be increasing in the proportion of infected (and insensitive to population), thus hitting NYC harder than the rest of New York State.

Additional observations. To fully understand the effect of the disease-related parameters to the optimal stabilizing lockdown-rate profile {z * l } and the economic cost, we implemented additional numerical experiments to analyze the sensitivity. The results are shown in SI Fig.6 . It can be seen that the value of z * l and the corresponding economic cost are sensitive to recovery rate and the initial growth rate but not to other parameters.

We also studied the relationship of the value of z * l and the structure of the underlying graph. We plotted the obtained z * l with respect to degree, home-stay rate, population, and employment in Fig.7 . We also implemented random permutation experiments (where we randomly permute one parameter while fixing everything else) in terms of degree, home-stay rate, population, employment and initial susceptible rate, and the results are shown in Fig.8 and Fig.9 . From these experiments, we found the distribution of z * l can be strongly affected by permutations of centrality, population, and the home stay rate. However, the distribution of z * l is not altered much by permuting employment and the initial susceptible rate. More details about these experiments are presented in SI Sec. 7.

From the empirical analysis of data in New York State, we hypothesize that the home-stay rate, degree centrality, and population are three major parameters that impact the optimal lockdown rate z * l of county-l. However, no inferences can be made about the effect of these parameters from empirical data because all of them vary together. To study how the value of z * l is related to these parameters, we implement experiments on synthetically generated data. We describe the experiments and results next, with full details provided in the supplementary SI Sec. 7.

Impact of degree centrality. To study the impact of degree centrality, we considered geometric random graphs (see SI Sec. 7.5 for other graphs). The population, the home-stay rate, and the initial susceptible rate of different nodes are set as the same values across all the nodes. The simulation results are presented in Fig.4a -b. We found that degree centrality only matters for the value of z * l when there exist hotspots (i.e., hubs node with very high degrees). Beyond such hotspots, the effect of degree centrality is essentially ignorable.

Impact of population. To study the impact of population, we fix all other model parameters and vary the population of the nodes. We again considered geometric random graphs, where node degrees are similar. The simulation results are presented in Fig.4c . We found that nodes with small populations are assigned smaller values of z * l , but once the population is large enough, z * l is almost independent of the population size. Impact of home-stay rate. To study the impact of home-stay rate, we fixed all other model parameters and tuned the home-stay rate of the nodes. Our simulation results based on the geometric random graphs were presented in Fig.4d . We found that z * l increases with increasing home-stay rate, which agrees well with our intuition.

3. Discussion. The main contribution of this paper is two-fold. Our first contribution is methodological: we give a modeling framework that gives rise to efficient methods for pandemic control through fixed lockdowns, with our main algorithm taking nearly linear time. The linear time nature of the methods allows us to scale up in a way that is not known for any other method. For example, at present data is not available to design lockdowns at the city level, but the method presented here can be scaled up to design a lockdown even at the neighborhood level for the entire United States. Although this is highly unlikely to happen for the COVID-19 pandemic, it may be of use in future outbreaks.

Our second contribution is to use our algorithms to observe counter-intuitive properties of lockdowns. In particular, we observe that a model of epidemic spread in New York State will tend to shut down outside of NYC more stringently that NYC itself, even if the epidemic is largely localized to NYC. We compared the lockdown found our model against exhaustive search of all two-parameter lockdowns which shut down NYC harder than the rest of New York State to verify that indeed it outperforms. We further found that this result is robust against significant perturbations to the travel matrix, the epidemic parameters, as well as the epidemic model.

While we have focused on simple models in this work, our methods can be applied to more complex models, such as the SIDHARTE model [31] , which contains more than two classes (asymptomatic/symptomatic) of infected people. If additional data sets become available on variations of activity by age and location in the future, it is possible to incorporate this as in [13] ; one could, for example, split each city into multiple nodes, with each node corresponding to a different age group residing in that city.

We next briefly discuss future directions. The main limitations of research on lockdown at the present time is the lack of available data. This limitation drove a number of the modeling choices made in the present work, as we describe next.

For example, it is tempting to divide trips into several different types (e.g., work, family, entertainment, etc), and argue that different types have different likelihood of leading to infections. Unfortunately, we are not aware of any data source which either counts such trips for different locations through the United States or estimates the infectivity differentials across types. One might further attempt to fit different transmissibility parameters to different counties, but again such data does not appear to be available at this level of granularity. In general, there are many ways to build sophisticated models but the primary limitation is lack of data.

Even the SafeGraph data we have relied on here is not without limitations. Indeed, SafeGraph estimates are based on cell phone data, and by definition do not sample people without cell phones. Further, we do not even know how many of the minutes recorded as travel outside might have been spent alone in a vehicle, with no possibility for transmission. However, as a counterpoint, mobility from cell phone data has been highly predictive in modeling COVID-19 spread as reported in [19, 32] . Finally, appropriately anonymized public data sets would allow us to better understand how people respond to lockdown and estimate a model with a more heterogeneous response compared to our model here. Future work is likely to be driven by fitting finer models as more data becomes available.

We conclude by mentioning that, because of all of these considerations, we have focused on qualitative patterns of the lockdown which hold across different models and parameter values. For example, as mentioned earlier, our finding that it is better to have a tighter lockdown in NYC holds when each entry of the travel matrix τ i j is perturbed with noise of variance up to 10τ 2 i j (see Supplementary Figure 18 , and a detailed explanation of this experiment in Section 10 of the Supplementary Information). The same finding is also unaffected by setting half of τ i j randomly to zero. Thus the precise values in the travel matrix do not appear to be important for this finding as long as the zero-nonzero structure remains very broadly similar. Finally, the same finding remains true in the SIS/SIR models. Thus our main empirical result is a robust property of the optimal stabilizing shutdown that holds even in the presence of severe data and model mis-specifications. 8 This manuscript is for review purposes only. Fig. 1 Framework of the optimal stabilizing lockdown design. a, the COVID-19 model we consider, which corresponds to Eq. (2.1). b, the travel pattern of a three-nodes network, where A represents a city with large population, B and C represent two suburbs with small population, τ i j represents the travel rate from location i to location j in a day. c, the infection flow pattern of the three-nodes network before the lockdown, a i j represents the infection flow as described in Eq.(2.1). Note that although no travel occurs between nodes B and C, the corresponding entries a 23 and a 32 are nonzero since people from these locations can meet each other in location A. d, the infection flow pattern of the same network after the lockdown. The widths of the edges in c, d are propotional to the value of a i j . When lockdown policy are implemented, the value of a i j will decrease, which lead to the control of the epidemics. e, the estimated number of cumulative cases of different lockdown polices. In this figure, "ours" represents the heterogeneous optimal stabilizing lockdown calculated by our method, "uniform lockdown" represents a uniform policy with the same economic cost as our lockdown. Random lockdown rate z l is randomly chosen from a uniform distribution U[a, b] with the lowerand upper-bounds a and b chosen such that the overall cost of this policy is the same as that of our policy. Uniformlybounded-decline lockdow implies a policy where the decay rate of the infections in each county is bounded by a constant α, where α is chosen such that the economic cost of this policy is the same as that of our policy. For this small network, the lockdown rates of our heterogeneous policy are z 1 = 0.21, z 2 = 0.06, z 3 = 0.06. By contrast, the lockdown rates of the uniform policy with the same cost are z 1 = 0.16, z 2 = 0.16, z 3 = 0.16. It can be seen that our policy outperforms all the other lockdown policies. 9 This manuscript is for review purposes only. In a, d, the disease parameters are set as in [7] , the decay rate α is chosen as 0.0231 which corresponds to halving every 30 days. In b, e, the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In c, f, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. Uniform lockdown, random lockdown, and uniformly-bouded-decline lockdown are defined as in . a-c, optimal lockdown rate z * l given by our method . d-f, uniform lockdown rate z l . g-i, random lockdown rate z l . j-l, uniformly-bounded-decline lockdown rate z l . Uniform lockdown, random lockdown, and uniformly-bouded-decline lockdown are defined as in Fig 1. In a, d, g, j, the disease parameters are set as in [7] , the decay rate α is chosen as 0.0231 which corresponds to halving every 30 days. In b, e, h, k the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In c, f, i, l, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. It can seen from a-c that the value of z * l for counties in NYC are relatively higher than other counties in New York State, which implies we should shutdown the outside of NYC harder than itself. Besides, it can be seen that such counter-intuitive phenomenon does not appear in any other lockdown polices. 11 This manuscript is for review purposes only. 10 Fig. 4 Numerical results of optimal lockdown rates on synthetic networks. a, b, the impact of the degree centrality on the geometric random graph. Each point represents a node in the network. It can be observed that centrality only matters for the value of z * l when there exist hotspots in all the three models. Without such hotspots, the effect of degree centrality is essentially ignorable. c, the impact of population. Each point represents a node in the network. Nodes with small populations are assigned smaller values of z * l . Surprisingly, once the population is big enough, it doesn't affect z * l too much. d, the impact of home-stay rate. Each point represents a node in the network. It can be observed that z * l increases as the home-stay rate increases, and in fact the home-stay rate has by far the largest influence on z * l compared to the other parameters we consider. 2. Related work. Our work is related to a number of recent papers motivated by the spread COVID-19, as well as some older work. Indeed, spatial spread of epidemic admits a natural network representation, where nodes represent different locations and edges encode traveling of residents between the locations. Such spatial epidemic network model has received attention in the studies of COVID-19 recently [45, 20, 11, 82, 67] . We begin by discussing several papers most closely related to our work.

Our paper builds on the results of [8] , which proposed a spatial epidemic transmission model and consider the effect of lockdowns; we use the same model of lockdown of [8] in this work. The major difference between this work in [8] is two-fold. First, we do not consider asymptotic stability in a model with births and deaths as our focus is on a shorter scale. Second, we propose new algorithms with improved running times; in particular, our main contribution is a linear time method that is applicable to the vast majority of cases we have considered. Similarly, the main difference of this work relative to [59] and the references therein are new algorithms (though the lockdown models differ somewhat), as well as the new observations on counter-intuitive phenomena satisfied by the optimal lockdown. Our work has some similarities with literature [13] , which divided the population into 18 compartments and found that population heterogeneity could significantly impact disease-induced herd immunity.

Control. An alternative approach would be to approach the lockdown problem using the techniques of optimal control. This approach is explored by a number of papers [47, 9, 17, 1, 26, 2, 12] .

As mentioned in the main body of the paper, our work has several main advantages over the traditional optimal control approach. The first is that we are looking for a fixed lockdown, whereas an optimal control based approach would offer a lockdown which varies for all time t. The second is that an approach based on optimal control would ask policymakers to repeatedly design lockdown relaxations when cases begin to decrease. As the public grows more impatient with lockdowns, political constraints could easily result in poor decision-making; one could argue this is what happened in the United States in 2020 [38] . A singlefixed lockdown implemented when the number of cases is growing and maintained until the epidemic is extinct does not have this problem. Finally, scalability is central to our results: our main result is a nearly linear time algorithm. This is not the case for the optimal control approach. For example, Khanafer and Ba∫ ar [47] wrote down the optimal control formulation for the SIS case and remarked that solving the resulting equations "is intractable." As Khanafer and Ba∫ ar pointed out, there are no methods that are guaranteed to find the optimal control efficiently. Nevertheless, in a number of recent works promising numerical results are obtained. For a direct formulation of the problem, one can turn to [17] , which is also in the same spirit as our work, in that it studies control of COVID-19 using, among other things, movement restrictions. This results in a mixed-integer non-linear programming problem. Solving such problems is generally intractable, so [17] used a genetic algorithm as a heuristic.

Another possibility might be to use indirect methods (i.e., relying on the maximum principle) to solve the optimal control problem. The same complexity considerations, however, come up in this context again. Most of the literature on the optimal control of epidemics over networks seems to use variants of the forwardbackward sweep method [53, 26, 50, 10] to solve the equations arising from the maximum principle, but it is known that this method can diverge even for simple examples [56] . Nevertheless, on many examples the method converges in reasonable time: for example, [53] reports excellent results on random networks of large size.

Recent papers studying control of COVID-19 using numerical optimal control methods are [1, 26, 2, 12] . These papers assumed a cost of a human life lost (sometimes set based on the average lifetime earnings) and considered the discounted total cost over the entire epidemic (alternatively, [1] considered the frontier of possible strategies over all possible ways to value life). This framework is conceptually different from ours: when modified to fit into our framework, strategies derived in this way will begin to relax the lockdown once the epidemic drops below a certain threshold as the number of lives lost comes to balance the cost of the lockdown, ultimately driving the epidemic to an endemic state; by contrast, our framework is designed to send the number of infections to zero.

More generally, the scalability of these approaches is unclear, due to their reliance on numerical methods without a clear convergence theory. For example, solving for the solution of this nonlinear optimal control problem as in [26] requires an iterative method closely related to forward-backward sweeping. Each step requires the numerical solution of a system of differential equations, a matrix inversion, and a maximization of the Hamiltonian, without any a-priori bounds on the total number of steps the procedure will take, or any guarantee that the procedure will converge.

3. Mathematical Background. The supplementary information will provide the proofs of the main results of the paper, as well as give details of many of our empirical and numerical results that were summarized in the main text. We begin with some definitions.

A matrix is called continuous time stable if all of its eigenvalues have nonpositive real parts. A matrix is called discrete time stable if all of its eigenvalues are upper bounded by one in magnitude. A central concern of this paper is to get certain quantities of interest (e.g., number of infected individuals) to decay at prescribed exponential rates. We will say that y(t) decays at rate α beginning at t 0 if y(t) ≤ y(t 0 )e −αt for all t ≥ t 0 . Note that the decay in this definition is not asymptotic but results in a decrease starting at time t 0 .

We will associate to every matrix A ∈ R n×n the graph G(A) corresponding to its nonzero entries: the vertex set of G(A) will be {1, . . . , n} while (i, j) ∈ G(A) if and only if A ji 0. Informally, (i, j) is an edge in G(A) when the variable j "is influenced by" variable i. We will say that A is strongly or weakly connected if the graph G(A) has this property.

Here c i are nonnegative scalars and C i are positive semi-definite matrices. It turns out that covering semidefinite programs can be solved considerably faster than general semi-definite programs. Indeed, the recent 17 This manuscript is for review purposes only.

paper [43] showed that to compute a fixed-accuracy additive approximation of the optimal solution takes

where ω is the exponent of matrix multiplication. With these preliminaries in place, we next turn to justifying the analytical claims made in the main body of the paper.

Balancing. The matrix balancing problems plays a fundamental role in our main results, and we briefly introduce it here. Given a nonnegative matrix P ∈ R n×n , we say it is balanced if it has identical row and column sums. The matrix balancing problem is, given a nonnegative P, to find a nonnegative diagonal matrix D such that DPD −1 is balanced.

The problem of matrix balancing is quite old; for example, an asymmetric version of this problem was introduced in the classic work of Sinkhorn and Knopp in the 1960s [76] . It is impossible to survey all the literature on matrix balancing and related problems, though we refer the reader to [42] . Recently, a powerful algorithm for matrix balancing was given in [22] . It was shown in that work that this problem can be solved in linear time, understood as follows: solving the problem to accuracy requires onlyÕ(nnz(P) log κ log −1 ) where nnz(P) is the number of nonzero entries in the matrix A, κ = D * max /D * min is the imbalance of the optimal solution, and the O hides logarithmic factors. Thus matrix balancing problems can be solved in nearly the same time as it takes to simply read the data, provided κ is bounded away from zero. In the event that we do not have an a-prior bound on κ, [22] give complexity bounds ofÕ(nnz(P) 3/2 ) andÕ(nnz(P)diam(A)) where diam(A) is the diameter of the graph corresponding to the matrix A.

We will summarize this complexity by saying that the running time "explicitly scales nearly linearly in the number of nonzero entries." The "nearly" comes from the logarithmic terms; the word "explicit" comes because the scaling also depends on the imbalance of the optimal lockdown κ, and one can construct families of examples where the κ will have some kind of scaling with network size.

In this section, we first present the details of the SIS model [70] as well as how the matrix A is constructed. Next, we justify the main theoretical achieved via eigenvalue bounds; that optimal lockdown for these models can be reduced to a covering semi-definite program (which can be solved in matrix multiplication time); and that, under the high spread condition, optimal lockdown for these models reduces to a matrix balancing problem (which can be solved in linear time).

4.1. Network SIS Model. is described by the following set of ordinary differential equations

Here β denotes the transmission rate, which captures the rate at which an infected individual infects others, γ denotes the recovery rate, and a i j captures the rate at which infection flows from the population at location j to location i. Becauseẋ i scales with 1 − x i , the SIS model assumes that everyone who is not infected is susceptible.

For simplicity of notation, we can stack up the coefficients a i j into a matrix as as A = [a i j ]. Then we can write the network SIS model asẋ

where 1 denotes the vector of all-ones while diag(u) makes a diagonal matrix out of the vector u.

It is desirable to have x i (t) → 0 for all i = 1, . . . , n, i.e., to have the infection die out. It is mathematically convenient to encode this into the following equivalent condition: we will require that there exists some linear combination of the quantities x i (t) with positive coefficients which approaches zero, which happens if and only if the matrix βA − γI is continuous-time stable [51, 4, 37, 48, 57] . To achieve an exponential decay rate α of each x i (t), we require that there exists a positive linear combination of x i (t) decaying at that rate, which is guaranteed if λ max (βA − γI) ≤ α (see formal proof in SI Sec. 4). Note that even though the network SIS dynamics is nonlinear, the asymptotic convergence nevertheless reduces to a linear eigenvalue problem.

A. We next describe how the matrix A is constructed. Our discussion will only be for the SIS case, as the COVID-19 case is similar.

Observe that the susceptible individuals at location i can be infected at location i as well as in other locations, the flow of susceptible population from location i to location l is (1 − x i )τ il . Besides, the rate of infection at location i is proportional to the fraction of infected people in the total population of location i. Then we can rewrite the SIS model aṡ

Let m(l) = n k=1 N k τ kl , then (4.2) can be written aṡ

As already remarked, this approach is not original to our work and is taken from [8] . Note that in COVID-19 case, via similar process, we can obtain the same a i j as in Eq. (4.3).

Recall that the network SIS model is given by the system of equations

where n is the number of nodes in the underlying graph and x i (t) is the proportion of infected individuals at node i. By contrast, the network COVID-19 model is given by

where now x a i , x s i are the proportion of infected/asymptomatic and infected/symptomatic individuals at node i. This is a system of 3n equations, with three equations per node of the network. Recall also the notation 19 This manuscript is for review purposes only. M(t) used to denote the bottom 2n × 2n submatrix of the above matrix. In the main text, we stated that stability and decay rate of the network SIS model is equivalent to the eigenvalues of the matrix A − γI, while the stability of the network COVID-19 model can be ensured by bounding the eigenvalues of M(t). We next give a pair of propositions formally justifying these assertions. Proposition 4.1. Suppose the matrix A is strongly connected. If λ max (βA − γI) ≤ −α, then for the network SIS dynamics of Eq. (4.4), there exists a positive linear combination of the quantities x i (t) that decays to zero at rate α starting at any time t 0 . Conversely, if λ max (βA − γI) > −α, then there exists an initial condition in x(0) > 0 so that every positive linear combination of the quantites x i (t) fails to decay at rate α. The idea behind these propositions is standard in control theory: to get your system to decay at a rate of e −αt , make sure your eigenvalues have real parts that are at most −α. For linear systems, this is guaranteed to work after a transient time, and for nonlinear systems, the situation is complicated. Fortunately, the nonlinear systems corresponding to network SIS and COVID dynamics have favorable properties, so that it is possible to draw conclusions about global behavior from the eigenvalues of a linear approximation at a point.

To get each of the quantities x a i (t), x s i (t) to converge to zero at an asymptotic rate of e −αt , it suffices to have an arbitrary positive combination of them decay at rate α. Interestingly, for the SIS case, this happens if and only if an eigenvalue condition is satisfied. In the COVID-19 case, the eigenvalue condition merely suffices to establish this.

We now turn to the proof of these propositions. Our first step is to restate the Perron-Frobenius theorem in a form that will be particularly useful to us.

Let us suppose P is a strongly connected matrix whose off-diagonal elements are nonnegative. Then there exists a real eigenvalue of P which is as large as the real part of any other eigenvalue of P. This eigenvalue is simple and the eigenvector corresponding to it is positive.

This lemma follows by observing that A + αI is nonnegative for a large enough choice of α, so we can apply the Perron-Frobenius theorem (in the form of Theorem 2 in Section 8.2 of [29] ) to A + αI.

Definition 4.4. We will use λ max (P) and v max (P) to denote the eigenvalue/eigenvector described in Lemma 4.3.

We next give a proof of Proposition 4.1.

Proof of Proposition 4.1. Suppose P = βA − (γ − α)I. Let λ max (P) and v max (P) be the corresponding eigenvalue/eigenvector pair. Note that, by our assumptions, λ max ≤ 0. We thus have that

where the second line used that β, v max , A, x are nonnegative while 1 − x ∈ [0, 1] n ; and the last line used that λ max ≤ 0. We conclude that v max x decays at a rate of α starting at any time. 20 This manuscript is for review purposes only.

On the other hand, suppose λ max (βA − γI) > −α. Observe that the linearization of Eq. (4.1) around the origin isẋ = (βA − γI)x. Let us write the Jordan normal form,

where J is the upper-triangular matrix λ 1 , . . . , λ n are the eigenvalues of βA − γI, with λ 1 being the Perron-Frobenius eigenvalue. Because the Perron Frobenius argument is simple, we have that the top Jordan block is 1 × 1; so that J and all of its powers have zero entries in the first row. Using the standard formula for the matrix exponential of a Jordan block, we next decompose as

where U t is an upper triangular matrix, depending on t, but having only zero entries in its first row. We now choose the initial condition x 0 = lT e 1 , where e 1 is the first basis vector, and l is a small-enough positive scalar; we'll discuss how small l has to be later. We then have T −1 x 0 = le 1 and therefore under the floẇ

This establishes the property we want for the flow of the linear systemẋ = (βA − γI)x, but we need to establish that the network SIS dynamics has the same property. This can be done via a particular form of the Hartman-Grobman theorem. Indeed, let y(t) be the network SIS trajectory starting from y 0 . The Harman-Grobman theorem, in the form of Theorem 3 of [39] , guarantees the existence of a homeomorphism h : R n → R n such that

and, as proved in [39] , because the network SIS dynamics are infinitely differentiable, we can further take h to be differentiable at the origin with the derivative at the origin equalling identity:

This implies that

which we rearrange as

Morever, for small enough ||x|| 2 , we have

This further implies that for small enough ||x|| 2 ,

This manuscript is for review purposes only.

With these observations in mind, we now choose the initial condition y 0 = h −1 (x 0 ), where x 0 = lT e 1 as above. We then have that

appealing in the last step to Eq. (4.6). We conclude the proof by arguing that our initial condition y 0 is nonnegative provided we choose l small enough. Indeed, we first argue that T e 1 is strictly positive. Indeed, observe that this is the first column of T , and since βA − γI = T DT −1 can be rewritten as (βA − γI)T = T D, we obtain that the first column of T is the Perron-Frobenius eigenvector of A, which is positive by Lemma 4.3. Finally,

by Eq. (4.7); and for small enough l, this has to be strictly positive by strict positivity of T e 1 . This concludes the proof.

We next give the proof of Proposition 4.2. Since, unlike in the SIS case, only one direction must be proven, the proof is straightforward.

Proof of Proposition 4.2. By Lemma 4.3, the matrix M(t 0 ) has a left-eigenvector v max which is positive. Let λ max be the corresponding eigenvalue; by assumption λ max ≤ −α. Let us define

where the second equation uses that v max and p(t) are nonnegative, and as a consequence of the nonincreasing of s(t), v max M(t)p(t) ≤ v max M(t 0 )p(t); while the final equation used λ max ≤ −α. We conclude that v max p(t) decreases at rate α starting from any time.

We now turn to the algorithmic question of designing an optimal lockdown. We will first present our main results. Next we will present a string of lemmas and observations which will culminate in the proof of the main Theorem. It is here that we will perform the reduction from the problem of computing the optimal lockdown to matrix balancing and covering semi-definite programs.

Our first step is to discuss an assumption required by one of our algorithms. The formal statement of the assumption is as follows.

Assumption 4.5 (High spread assumption).

1. In the network SIS model, we have diag(B C) ≥ γ 2. In the network COVID-19 model, we must have

This manuscript is for review purposes only.

To see why this condition is satisfied in a "high spread" regime, note that A = CB and, given our choices of C and B in Eq. (2.4) , the entry C ii B ii corresponds to the spread of the epidemic in location i purely from the same-location trips of residents of location i. This has to be bigger than the natural rate γ at which people recover. In other words, the natural rate of spreading of the epidemic has to be high everywhere. Assumption 4.5 is actually somewhat looser than this, as what must be bounded below by γ is diag(B C), which is a sum of O(n) products, only one of which is C ii B ii . Note that this is not quite the same as requiring that A ii ≥ γ since A +CB , and we've flipped the order of multiplication on the condition. In the COVID-19 case, the interpretation is similar: the recovery rates r a , r s need to be small relative to s i (t 0 )C ii B ii as well as the spread parameters β a , β s , , though the relation is now more involved.

Main theoretical contribution. Our main theorem provides algorithms for the unconstrained and constrained lockdown problems in the cases when Assumption 4.5 does and does not hold. Our key contribution is to give an algorithm for optimal stabilizing lockdown whose complexity has an explicit scaling which is nearly linear in the number of nonzero entries of the matrix A. That is to say, not only can the optimal heterogeneous lockdown be computed exactly, but doing so takes nearly as much time as just reading through the data. Theorem 4.6. Suppose the graph corresponding to positive entries of the matrix A is strongly connected and s(t 0 ) > 0. Then:

1. The unconstrained lockdown problem for both SIS and COVID-19 models can be reduced to matrix balancing.

2. Suppose further Assumption 4.5 holds. Then the constrained lockdown problem is equivalent to the unconstrained lockdown problem and consequently is also reducible to matrix balancing. Next, we will prove Theorem 4.6. Our starting point will be the following lemma on a "splitting" of a positive matrix. 

This lemma is a small variation on a well-known fact: usually, parts 1 and 2 are stated for strictly stable matrices, in which case all the inequalities need to be strict (see [14] , Theorem 15.17 for the strict version of part (i) and Proposition 1 of [71] for the strict version of part (ii)). The non-strict version additionally requires that the matrices be strongly connected, which is not needed for the nonstrict version of this problem. Note that we do not claim that any part of this lemma is novel. For completeness, we nevertheless give a proof next. 23 Proof of Lemma 4.7. For part (1) , observe that Pd ≤ 0 if the same as Pdiag(d)1 ≤ 0, which is equivalent to diag(d) −1 Pdiag(d) being a matrix with nonpositive row sums. Since we have assumed that P has nonnegative off-diagonal elements, this implies that the diagonal elements of diag(d) −1 Pdiag(d) are non-positive. By Gershgorin circles, diag(d) −1 Pdiag(d) must be continuous-time stable. Since its eigenvalues are the same as the eigenvalues of P, we conclude that P is also continuous-time stable. Similarly, for part (2), if such a d exists, then diag(d) −1 Bdiag(d) is a nonnegative matrix whose row sums are upper bounded by one, and must be discrete-time stable again by Gershgorin circles. This proves the "if" part of parts (1) and (2).

For the "only if" parts, suppose P is continuous-time stable and strongly connected with non-negative off-diagonal entries. By Lemma 4.3, we have that there exists a positive vector d such that Pd = λd and λ ≤ 0; this proves the "only if" of part 1. For part 2, if B is discrete-time stable and strongly connected, we have that Bd = λd for the Perron-Frobenius eigenvalue d of P, which is positive. Since now λ < 1, this proves the "only if" statement of part 2.

For part (3), suppose first that P = L − D is continuous-time stable. Since P = L − D has nonnegative off-diagonal elements, we can apply part (1) to observe that continuous time stability of P is equivalent to existence of d satisfying

We now to multiply both sides by D −1 and obtain that there exists a d satisfying

Note that we used that D −1 is nonnegative to multiply both sides by D −1 . Since D −1 L is strongly connected, the last equation is, by part (2), exactly the statement that

Conversely, suppose D −1 L is discrete-time stable (note that we cannot simply reverse the above chain of implications since we do not assume that D is elementwise nonnegative; thus we can multiply a linear inequality by D −1 but not necessarily by D). Let λ max , v max denote the Perron-Frobenius eigenvalue/eigenvector pair of D −1 L, guaranteed to exist by Lemma 4.3; note that v max is strictly positive. We then have

where the last step used that λ max ≤ 1. It follows that Pv max = (L − D)v max ≤ 0, and since, as already observed, v max is positive, we obtain that P is continuous-time stable by applying part (1).

We will later need to interchange the order of products while still preserving the condition of being strongly connected. To that end, the following lemma will be useful. Lemma 4.8. Suppose U, V are two nonnegative n × n matrices with no zero rows or columns such that VU is strongly connected. Then UV is strongly connected.

Proof. Consider a directed bipartite graph G on 2n vertices, with vertices l 1 , . . . , l n and r 1 , . . . , r n denoting the two sides of the bipartition, defined as follows. If U i j > 0, then we put an edge from l j to r i ; and if V i j > 0 we put an edge from r j to l i . Let G 1 be the graph on l 1 , . . . , l n where we put the directed edge from l i to l j if there is a path of length two from l i to l j in G. Likewise, let G 2 be the directed graph on r 1 , . . . , r n such that we put an edge from r i to r j whenever there is a path of length two in G. 24 This manuscript is for review purposes only.

Then the strong connectivity of VU is equivalent to having G 1 be strongly connected: indeed, (VU) ab > 0 if and only if there exists a link from b to a in G 1 . Similarly, the strong connectivity of UV is equivalent to having G 2 be strongly connected. We will show that if G 1 is not strongly connected, neither is G 2 . The converse is established via a similar argument.

Indeed, suppose G 1 is not strongly connected. That means there exists a proper subset of the vertices, say L 1 = {l 1 , . . . , l k }, with no edges outgoing to L c 1 = {l k+1 , . . . , l n }. Let r 1 be the set of out-neighbors of L 1 in G. Then we must have that r 1 is a proper subset of {r 1 , . . . , r n } (for otherwise, the assumption that V has no zero rows/columns would contradict no edges going from L 1 to L c 1 in G 1 ) and the set of out-neighbors of r 1 in G is contained in L 1 .

But since (i) r 1 is a proper subset of the right-hand side (ii) the out-neighbors of r 1 in G are contained in L 1 (iii) the out-neighbors of L 1 in G are contained in r 1 , we obtain that there are no edges leading from r 1 to r c 1 in G 2 . This proves G 2 is not strongly connected.

With these preliminary lemmas in place, we now turn to core of our reduction. We begin with the network SIS dynamics, where we will reduce the task of finding an (unconstrained) optimal lockdown to the problem of matrix balancing defined earlier. The reduction will go through several "intermediate" problems, the first of which as follows. Definition 4.9. We will refer to the following as the stability scaling problem: given a nonnegative strongly connected matrix P and positive diagonal matrix D, find positive scalars q 1 , . . . , q n minimizing

The utility of this definition should become clear after the following lemma.

Lemma 4.10. Suppose A = CB is strongly connected. The minimum cost lockdown problem for the SIS network can be written as stability scaling. Under Assumption 4.5, the constrained lockdown problem for the network SIS problem can also be written as stability scaling.

Proof. We consider the unconstrained lockdown problem first. Recall that, in the SIS model, we are looking for a minimum cost positive vector z such that

Since the nonzero eigenvalues of a product of two matrices do not change after we change the order in which we multiply them, this is the same as requiring that

is continuous-time stable. This is exactly the stability scaling problem provided we have two additional conditions. The first condition is that α < γ (because the diagonal matrix subtracted needs to be positive).

The second condition is that B C should be strongly connected. But by Lemma 4.8, the second condition is true because we assumed that A = CB is strongly connected.

Observe that the reduction resulted in an instance of stability scaling with matrix P = B C. Since we have assumed A = CB is strongly connected, we can apply Lemma 4.8 to obtain that P is strongly connected as needed. Under the assumption α < γ, the matrix (γ − α)I is a positive diagonal matrix as required.

We thus have the reduction we want, from optimal lockdown to stability scaling, assuming α < γ. But what if α ≥ λ? In that case, we claim that the minimum cost lockdown problem does not have a solution. 25 This manuscript is for review purposes only.

Indeed, since at the optimal solution we must have all z * i > 0, we have that Cdiag(z * )B is an irreducible nonnegative matrix and its Perron-Frobenius eigenvalue is strictly positive by the Perron-Frobenius theorem [29] . Consequently, Cdiag(z * )B − (γ − α)I has a positive eigenvalue and cannot be continuous-time stable.

Finally, we consider the constrained version. We argue that, under Assumption 4.5, the optimal solution to the constrained lockdown problem will have z * i ≤ 1 for all i, so we can simply drop the constraint. Indeed, first observe that we can assume α ≤ γ, else the problem does not have a solution as explained above. Now suppose that z * j > 1; then diag(z * )B C − diag(γ − α)I has a nonnegative j'th row with at least one positive entry in that row. Indeed, the off-diagonal entries in the j'th row are clearly nonnegative, while the diagonal entry is nonnegative by Assumption 4.5.

By Lemma 4.7, diag(z * )B C − diag(γ − α)I cannot be continuous-time stable. As matrix Cdiag(z)B − (γ − α)I has the same nonzero eigenvalues as the matrix diag(z * )B C − diag(γ − α)I, it can not be stable either.

Lemma 4.11. In the unconstrained case, the minimum lockdown model for COVID-19 dynamics can be reduced to stability scaling provided A is strongly connected and s(t 0 ) > 0. In the constrained case, the same holds under Assumption 4.5.

Proof. Let us begin by assuming that (4.9) α < min(r s , + r a ).

We will revisit this assumption later. We need to make the matrix (4.10)

where we have introduced the notation that A z = Cdiag(z)B . Let us write

We next apply part 3 of Lemma 4.7 to get that A 0 is continuous-time stable if and only if D −1 L is discrete time stable. To do this, however, we need to verify that D −1 is elementwise nonnegative, and that both P and D −1 L are strongly connected. This easily follows from observing that

and recalling that s(t 0 ) > 0 as well as Eq. (4.9).

To summarize, we have shown that equivalently we need to make sure that B = D −1 L is discrete-time stable. Since the eigenvalues of a matrix which is the product of two matrices do not change after changing 26 This manuscript is for review purposes only.

the order of the product, this is the same as having LD −1 be discrete-time stable. But we have the simple expression (4.11)

. But the eigenvalues of LD −1 are just the eigenvalues of its top n × n block. Thus we obtain that LD −1 is stable if and only if the matrix diag(s(t 0 ))Ab 1 is discrete-time stable. Plugging in A = Cdiag(z)B , we now need that

We next appear to part (3) of Lemma 4.7 again, using D = diag(s(t 0 )) −1 b −1 I to obtain that we need

We can apply part (3) of Lemma 4.7 as strong connectivity follows from z being elementwise positive and strong connectivity of A = CB . But the last condition is equivalent to Equation (4.8), which we've already shown how to reduce to stability scaling.

We conclude the reduction by observing that, once again, this results in an instance of stability scaling with the matrix P = B C, which is strongly connected because A = CB is strongly connected, allowing us to apply Lemma 4.8. Further, the matrix diag(s(t 0 )) −1 b −1 1 I is a positive diagonal matrix due to the assumptions that s(t 0 ) > 0 and the assumption of Eq. (4.9).

For the constrained case, we argue that under Assumption 4.5, the optimal solution will have all z * i ≤ 1, so we can equivalently consider the unconstrained case. Indeed, suppose e.g., that z * j > 1. In that case, the matrix diag(z * )B diag(s(t 0 ))C has its ( j, j)'th entry at least b −1 1 . Since the j'th row of diag(z * )B diag(s(t 0 ))C is nonnegative, and since B diag(s(t 0 ))C is strongly connected by Lemma 4.8, we have that the j'th row of diag(z * )B diag(s(t 0 ))C is not zero. By by the same argument as we made in the SIS case, the matrix diag(z * )B diag(s(t 0 ))C − b −1 1 I cannot be continuous-time stable. This implies that Eq. (4.13) cannot be satisfied.

Finally, we conclude the proof by revisiting the assumption we made in the very beginning, namely the assumption of Eq. (4.9). We now argue the problem of optimal lockdown has no solution if that equation fails. Indeed, in that case, the matrix A 0 has either its first n × n block nonnegative, or its last n × n block nonnegative. We show that the problem has no solution in the first case; the other case has a similar proof.

We want to argue that having d > 0 such that A 0 d ≤ 0 cannot occur if the top n × n block of A 0 is nonnegative; by Lemma 4.7 part (1) this is enough to show A cannot be stable regardless of z. We can partition d = [d 1 , d 2 ], and we immediately see that we must have that d 1 = 0, since, by the strong connectivity of diag(s(t 0 ))A z , every entry of d 1 is multiplied by some entry in the top left n × n submatrix of A 0 . Applying the same argument to the "top right" n × n block of A 0 , get that d 2 = 0, and this is a contradiction. This concludes the proof.

To recap where we are: we have just reduced the unconstrained lockdown problem for the SIS and COVID-19 models to stability scaling; further, the constrained lockdown problem was reduced to the same under the high-spread condition.

We need to make the following remark: provided the matrix A = CB was strongly connected and s(t 0 ) > 0, both of these reductions arrived at a stability scaling problems Our next step is to reduce the stability scaling problem to a new problem, defined next, which has a somewhat simpler form. 27 This manuscript is for review purposes only. Definition 4.12. Given a nonnegative matrix A, the problem of finding scalars l 1 , . . . , l n such that A − diag(l 1 , . . . , l n ) is continuous-time stable while minimizing n i=1 c i l i will be called the diagonal subtraction problem.

Lemma 4.13. Stability scaling can be reduced to diagonal subtraction.

Proof. Indeed, starting from min q 1 >0,...,q n >0 n i=1 c i q −1 i such that diag(q 1 , . . . , q n )P − D is continuous time stable, we apply Lemma 4.7, part (c) to obtain that the constraint is the same as

Here we used crucially that the matrix P is nonnegative and strongly connected and that the diagonal matrix D is positive, ensuring that the assumptions of Lemma 4.7, part 3 are satisfied. We next use the same Lemma 4.7 part (c) again to obtain that this the last constraint is identical to P − diag(q 1 , . . . , q n ) −1 D is continuous time stable For simplicity, it is convenient to introduce the notation δ i = q −1 i , i = 1, . . . , n. In terms of these new variables δ i , we have the problem With all these preliminaries in place, we are now ready to turn to the proof of our main theoretical resu.t

Proof of Theorem 4.6, parts 1 and 2. In Lemma 4.10, we have shown how to reduce the optimal lockdown problem for SIS and COVID-19 dynamics to stability scaling. In the subsequent Lemma 4.13, we showed how to reduce stability scaling to diagonal subtraction. We now prove Theorem 4.6 by reducing diagonal subtraction to matrix balancing.

Indeed, we need to solve min l 1 >0,...,l n >0 n i=1 c i l i such that A − diag(l 1 , . . . , l n ) is continuous time stable, 28 This manuscript is for review purposes only.

where A is a nonnegative matrix. We use part (a) of Lemma 4.7 to write this ass (4.14) min

The difficulty here is that we are optimizing over l and d, and the constraint includes a product of the variables. We try to make this into a simpler problem by introducing the variables µ i = l i d i , i = 1, . . . , n. We change variables from (l, d) to (µ, d) to obtain the equivalent problem:

(4.15) min

In this reformulation, the constraint is linear, and the nonlinearity has been moved to the objective.

Next, consider what happens when we fix some d > 0 and consider the best µ for that particular d. Because c i > 0 for all i, and d > 0 , we want to choose each element µ i as small as possible. The best choice is clearly µ = Ad; no component of µ can be lower than that by our constraints. Because A is a nonnegative matrix without zero rows by strong connectivity, this results in a feasible µ > 0.

Thus we can transform this into a problem which optimizes over just the variables d 1 , . . . , d n :

Our next step is to change variables one more time. Since d is a positive vector, it is natural to write d i = e g i . In terms of the new variables g 1 , . . . , g n , we just need to minimize the function f (g 1 , . . . , g n ) defined as min f (g 1 , . . . , g n ) := min g 1 ,...,g n n i=1 n j=1 c i a i j e g j −g i .

We now have an unconstrained minimization of a convex function. In particular, if we can find a point where ∇ f (g 1 , . . . , g n ) = 0, we will have found the global optimum. Let us consider the k'th component of the gradient of f (g 1 , . . . , g n ). We have the equation

c i a ik e g k −g i (4.16)

Let D g = diag(e g 1 , . . . , e g n ) and D c = diag(c 1 , . . . , c n ). Observing that [D −1 g D c AD g ] uv = e −g u c u a uv e g v = c u a uv e g v −g u , 29 This manuscript is for review purposes only.

we have that Eq. (4.16) can be written as

In other words, the matrix D −1 g D c AD g needs to have it's k'th column sum equals to its k'th row sum. Thus provided we can find a balancing of the matrix D c A, we will have found the minimum of f (g 1 , . . . , g n ) as needed.

However, a matrix X can be balanced if and only if the underlying graph G(X) is strongly connected [46] . Since we have assumed that c i are all positive and A is strongly connected, we have that D c A can be balanced. We conclude that the unique minimum of f (g 1 , . . . , g n ) can be recovered from the balancing of this matrix. This concludes the reduction of diagonal subtraction to matrix balancing.

We next turn to the third part of Theorem 4.6. We will not be using any of the reductions used to prove the first two parts of the theorem. The first steps of our proof are very similar to the proof of the main result of [8] , which gave a semidefinite formulation of a lockdown problem ensuring asymptotic stability in a model involving births and deaths. We diverge from [8] when we write the problem as a covering semi-definite problem as described in Section 3.1 and analyze its complexity.

Proof of Theorem 4.6, part 3. We can simply repeat the proof of Lemma 4.10 to reduce the constrained lockdown problem for network SIS and COVID-19 processes to stability scaling; except that, without Assumption 4.5, we now have to keep the constraint that z ∈ [0, 1] n . Indeed, note that the only place where Assumption 4.5 was used in those two proofs was to argue that we can omit that constraint.

In this way, we can reduce the constrained lockdown problem to the following: In this SIS case, a = 1 while q = γ − α. In the COVID-19 case, from Eq. (4.13) we have that a = s(t 0 ) while q = b −1 1 . Note that this optimization problem is only over the variables z 1 , . . . , z n ; all other quantities appearing in it are parameters.

Suppose next that matrices C and B are chosen according to Eq. (2.4). In other words, we have the problem

where, recall, the matrices D 1 , D 2 are defined immediately after Eq. (2.4).

Since the nonzero eigenvalues of a product of two matrices do not change when we flip the order in which they are multiplied, the second constraint is equivalent to diag(z 1 , . . . , z n )D 1 τ D 2 diag(a)τ − qI is continuous time stable . 30 This manuscript is for review purposes only.

Applying part (3) of Lemma 4.7, we get that this is the same as q −1 diag(z 1 , . . . , z n )D 1 τ D 2 diag(a)τ is discrete time stable .

Here we used the strong connectivity and positive diagonal of the matrix τ. Next, applying the same lemma again, we further get the equivalence to

This can be simplified further by observing that a symmetric matrix is stable if and only if it is non-positivedefinite. Indeed, let us introduce the notation u i = (D 1 ) −1 ii qz −1 i , and Q = τ D 2 diag(a)τ. We can therefore write our problem as As mentioned earlier, this is known as a "covering SDP." In the recent paper, [43] , it was shown that it can be be solved, up to various logarithmic factors, in matrix multiplication time. We next discuss in detail the results in [43] and how they are applicable to Eq. (4.17). Specifically, in [43] , an algorithm for checking whether the system of equations

is feasible was given; here C i , P i are arbitrary nonnegative definite matrices. To begin applying this to our problem, let us choose δ = 1/2 and P i = 1 2n u −1 c i 11 for all i. In that case, we are checking the feasibility of (4.18)

The running time of the algorithm from [43] for checking whether this is feasible is O( i nnz(C i ) + n ω ), where the O hides a multiplicative factor which is a power of logarithm is O n 2 max i λ max (C i )/λ max (P i ) (see discussion under "Main Claim" in [43] , and note that, in our case, the number of matrices and the dimension are proportional). In our case, we will choose C i = e i e i + e n+i e n+i so that

This manuscript is for review purposes only.

and plugging the definition of c i , we have that the factor being hidden in the tilde is a power of logarithm of

To summarize, we have discussed how long it takes to check feasibility of Eq. (4.18). The conclusion is that the running time O( i nnz(C i ) + n ω ), where the the tilde notation hides a factor that is logarithmic in the problem parameters n, u −1 , min i j c −1 i N j τ ji . Now let us transition from just checking feasibility of Eq. (4.18) to minimizing n i=1 c i x i subject to the constraint i x i C i 0, x ≥ 0. We can, of course, find an -additive approximation of this problem with a lot of feasibility checks.

Let us suppose that we actually know that x * i ≥ l i for some l i (we do, in fact, know this because Eq. (4.17) forces

If c * is the optimal cost, then we need to do O log(c * / i c i l i ) + log −1 feasibility checks to get an additive approximation to the optimal cost. Upper bounding c * ≤ nc max max i x * i , we see that we need to do

feasibility checks. Now we cannot quite apply this directly to Eq. (4.17), since one issue remains. That issue is that we have discussed an algorithm to handle the constraint i u i B i ≥ I; but in Eq. (4.17), the right-hand side is not the identity matrix.

However, the right-hand side of Eq. (4.17) is symmetric, so it can be decomposed as UΛU for orthogonal U and diagonal Λ. Since A B implies ZAZ ZBZ , by choosing Z = Λ −1/2 U we can equivalently rewrite Eq. Letting C i = Λ −1/2 U B i UΛ −1/2 , we thus obtain a problem to which the algorithm of Eq. ( [43] ) is applicable. Unfortunately, in the transformation the number of nonzero entries in the matrix C i changes. Nevertheless, 32 This manuscript is for review purposes only. as explained in [43] (in footnote 6 in the arxiv version), despite the change of variables, the scaling in the running time remains with the number of nonzero matrices in the matrices B i . However, we now need to add the complexity of computing the eigendecomposition of the right-hand side of Eq. (4.17) to out complexity bounds.

We can now put together everything we have observed above. The total complexity scales as the complexity of matrix multiplication and computing the eigendecomposition of a symmetric matrix, which we takes O(n 3 ) exact arithmetic operations [65] . This is multiplied by a power of a logarithm in the quantities

All of these are listed in the theorem statement, except max i x * i , u −1 ,

Let's upper bound each of these four quantities. Since After Eq. (4.17), we switched from u to x as we began discusisng the results of [43] , we have that

Next, the smallest u will be is

Thus we can "kill two birds with one stone" by first upper bounding u −1 as

and this also upper bounds i c i N visit i by a polynomial in the same quantities and q −1 . Next,

, which is in turn a polynomial in 1 λ min (Q)

,

Finally, since Q = τ D 2 diag(a)τ, we have that

This manuscript is for review purposes only.

To summarize, the arithmetic complexity of computing an additive approximation to the optimal solution is O(n 3 ) times a polylog in the variables

Finally, this is to compute and -approximate solution of the modified cost function in Eq. (4.17) which omits a factor of q −1 . To account for this, we can simply divide by a factor depending on q; but since q = b −1 1 and both b 1 , b −1 1 are already in the list of quantities under the poly log term, this does not change anything.

We now make an explicit a connection between the constraint we use on the growth rate of the epidemic, namely the constraint λ max (M(t 0 )) ≤ −α, and the basic reproduction number R(t 0 ), defined as the expected number of secondary cases produced in a completely susceptible population by a typical individual at time t 0 . In particular, we show that our constraint is completely equivalent to an upper bound on the reproduction number. Proof. The argument essentially reprises the proof of Lemma 4.11. Indeed, in that lemma we considered the matrix M(t 0 ) = β a diag(s(t 0 ))A z − − r a β s diag(s(t 0 ))A z −r s .

In Eq. (4.12), it was shown that the constraint λ max (M(t 0 ) ≤ α was equivalent to the constraint λ max (diag(s(t 0 ))A z b 1 (α)) ≤ 1.

Here we write the expression b 1 (α), unlike in Eq. (4.12), this way to highlight the dependence on α. Using the definition of b 1 (α) from that lemma, the last condition can be rewritten as

We thus have that an upper bound the condition λ max (M(t 0 )) ≤ −α is equivalent to the upper bound ρ(diag(s(t 0 ))A z ) ≤ f (α), where f (α) is a monotonically decreasing function of α. We next argue that a similar finding holds for the constraint R(t 0 ) ≤ r. 34 This manuscript is for review purposes only.

To do this, we need to express the reproduction number in terms of the matrix M(t 0 ). We use an expression from [80] in the form discussed in the recent paper [77] . Indeed, Definition 3 of [77] shows 4 , referring to [80] , that if we split

However, we have already computed this spectral radius in Eq. (4.11):

We have thus shown that the constraint R(t 0 ) ≤ r is equivalent to the condition ρ(diag(s(t 0 )A z ) ≤ g(r), where g(r) is a monotonic function of r. Since both constraint on R(t 0 ) and λ max (M(t 0 )) are equivalent to constraints on ρ(diag(s(t 0 )A z ), they are equivalent to each other. Finally, inspecting Eq. (5.1) and Eq. (5.2), we see that α = 0 corresponds to r = 1, which is as one might expect. On the other hand, we see that as α approaches the endpoint of its domain, min(r s , + r a ), we have that the corresponding r approaches zero. This concludes the proof.

We now describe how our data was obtained and how our simulations were conducted in a higher level of detail compared to the discussion in the main body of the paper.

Disease parameters. The disease parameters (β s , γ, r a , r s , , β a ) are essential to derive the optimal lockdown rate z. To provide valid results and fully understand the role of the disease parameters, we use three different groups of estimates provided in recent literature [8] [31] [7] to construct the disease parameters, respectively. However, none of those models quite match our two-state COVID-19 model, which is why one needs to be careful in re-using the parameters estimated in those models.

We directly use the recovery rate γ, r a , and r s from these literature. The reason is that these quantities have the interpretation that infected people recover at a rate γ (or r a , r s in the COVID-19 case) per unit time, which should carry over from model to model. Similarly, we will also reuse the parameter from literature [8] [31] which represents the rate that an infected individual develops symptoms.

In the SIS case, we choose the constant ζ (from definition of matrix D 2 in Eq. (2.4)) to match the initial growth of the models in the literature rate; the "initial growth rate" is λ max (A) − γ. We do likewise in the SIR case.

In the case of COVID-19, we also need to choose the transmission rates β a , β s , which we do as follows. Our first step is to let β a =αβ s and assume we can reuseα from the existing literature, as this scalar measures the transmission rate difference of symptomatic individuals and asymptomatic individuals. Thus we only need to decide how to choose β s . Our second step is to choose β s to match the initial growth rates of the models in the literature. This is λ max (M(t 0 )), where M(t 0 ) is defined in (2.2). Note that as a consequence of this procedure, the matrix A used for in COVID-19 case may not match the matrix A used in the SIS case.

Finally, we remark that the initial growth rate of the models in literature [8] [31] [7] can also be computed with the similar method. The resulting parameter values offered by literature [8] [31] [7] are summarized in SI Table 2 . Our procedure is essentially the same as matching the R 0 of the models (see e.g., [25, 40] ).

One caveat is the in the model from [7] , the parameters r a , r s , andα are missing. In our experiments using data from [7] , we simply choose r a = r s = γ, and letα = 0.6754 as in [31] .

To define the populations of each node in the network, we adopt the 2010 Census Bureau data [15] at the level of the counties in the New York state. These populations (N i ) will be used to construct the matrix A as in (2.4) in the three models.

Economic coefficients (c i ). In the cost function, c i represents the economic coefficients, which captures the cost of closing down site or location i. As the economic activity closely related to the number of employees, we let c i be proportional to the number of employees in location i. Specifically, suppose e i is the number of employees of node i, e max is the maximum of the number of the employees in the network, we define c i = e i /e max . We use LEHD Origin-Destination Employment Statistics (LODES) dataset [61] to obtain the number of employees in each county of NYS. Particularly, we use the Residence Area Characteristics (RAC) type of this dataset from [61] . This dataset is based on census at the block geographic level; hence we aggregate the data to the county level.

Travel rate (τ i j ). To construct matrix A for the three models, we need travel rate matrix τ, where τ i j represents the rate at which an individual travels from location i to location j. We use the Social Distancing Metrics dataset [74] from SafeGraph to generate τ. This dataset was collected using a panel of GPS pings from anonymous mobile devices, and it is based on Census Block Group levels. For each device/individual, the dataset identifies a "home" CBG, and the median daily home-dwell-time is provided for each CBG. Additionally, this dataset provides the daily number of trips that the people go from their home CBG to various destination CBGs.

In our empirical simulations, we only consider the network of New York State (i.e., we do not consider the trips to places outside the New York State). For each node, we aggregate the number of trips to the county level and obtain the number of trips from one node to another. We can also obtain the home-dwell-time of each node as the median of the home-dwell-time among all the CBGs (daily median home-dwell-time) in this county. Then, we define

This manuscript is for review purposes only.

where h i is the home-dwell-time of node i (measured in minutes), k ia is the number of trips from node i to node a. We divide h i by 1440 because the latter is the total number of minutes in a day.

Initial susceptible rate s(t 0 ). For the network SIR model and COVID-19 model, we need the initial susceptible rate s(t 0 ) to derive the optimal lockdown rate. For each node of our network, we can get the cumulative confirmed cases on each day from the dataset provide in [34] . However, recent studies [58] suggest that many people are infected with COVID-19 but not showing symptoms. To account for this, we divide the cumulative confirmed cases by a reporting rate to estimate the actual total infections. For the reporting rate, we use the estimate 0.14 provided in recent literature [8] [41] . Then we obtain s i (t 0 ) = 1 − I i (t 0 ) 0.14N i , where I i (t 0 ) represents the number of the cumulative confirmed cases of node i at time point t 0 and, as before, N i is the number of the population of node i.

Initial recovery rate and initial symptom rate. To estimate the rate of active cases and the rate of cumulative cases for SIS model and SIR model, we need initialize the recovery rate. For COVID-19 model, besides the recovery rate, we also need initialize the rate of symptomatic individuals and the rate of asymptomatic individuals. Since we do not consider individuals who died for the pandemic in all three models, in our simulation, we simply regard these people as recovered patients. We get the number of cumulative death cases in each county of New York State on April 1st from New York Times [78] . For the truly recovered people of the pandemic, unfortunately, we can not find specific numbers for each county in New York State. However, we learn from [23] that the total number of recovered cases in USA on April 1st, 2020 is 8878, and the total number of cumulative cases in USA on April 1st, 2020 is 215215. Therefore, we initialize the recovery rate of county i as:

where D i (t 0 ) represents the number of cumulative death cases in county i at time point t 0 , we divide 0.14 as we also need consider the reporting rate. For the initial symptom rate, we use the estimation in literature [8] [41] again, we assume 86% of active cases are asymptomatic individuals, and the remaining active cases are symptomatic individuals. The number of active cases are obtained as the difference of the cumulative cases and the recovered cases.

7. Numerical Calculations. We now describe the details of the synthetic experiments we have performed. Because in the empirical data all the variables of interest co-vary together, synthetic experiments are necessary to dis-entangle the effect of aspects of graph variation. Figure 1e . To clearly demonstrate the optimal lockdown issues we consider and compare our non-uniform lockdown policy and the uniform lockdown policy, we implemented a synthetic experiment on a simple three-nodes network as shown in Figure 1b . This specific network consists of a city (A) with large population, and two suburbs (B, C) with small population (the population as well as other data we used are given in Supplementary Table 1 ). We will assume employment is proportional to the population so that the economic cost coefficient c i is also proportional to the population. Moreover, we assume there exists no direct trips between location B and location C. For the choice of the travel rate matrix τ, we choose a matrix that is similar to the New York State data, but with rounder numbers; specifically, we define

where h i is the home-dwell-time of node i, k ia is the number of trips from node i to node a, and we let This manuscript is for review purposes only. from [7] where considers a SEIR model; we set the missing parameters asα = 0.6754 (where β a =αβ s ) and r a = r s = γ. We still choose the target decay rate α = 0.0231 that corresponds to halving every 30 days. It can be seen from Figure 1e that our policy outperforms the uniform lockdown policy in all of the three locations. Besides, we can see that our optimal lockdown policy tends to shutdown the city more stringently than the suburbs even if the epidemic mainly happens at the city. We will discuss the details about this counter-intuitive phenomenon in Section 2.2.

To understand the impact of various disease parameters, we implement sensitivity experiments. In each experiment, we vary the value of one parameter while fix the values of the others. The normal values of γ, , and initial growth rate are chosen as in [7] , and the normal valueα is chosen as in [31] (α is not provide in [7] ). When analyze the impact of γ, we let the decay rate α = 0.02. In other cases, we let the decay rate α = 0.2r s = 0.04, since it is impossible to achieve a decay rate better than γ, so we make it our goal to reach halfway there. For simplicity, we assume γ = r a , and r s =γr a , where the normal value ofγ = 1. To compare the economic cost of the obtained optimal lockdown policy and the uniform lockdown (decay matching) policy for each scenario. We define efficiency as the economic cost ratio of the optimal lockdown policy to the best uniform lockdown policy, i.e.,

where z uni * represents the uniform lockdown rate. Here, the uniform lockdown (decay matching) implies the uniform policy which decays at a rate greater than equal to α. SI Fig 6 shows the corresponding experimental results.

Observe each column of SI Fig 6, it can be seen that the value of the optimal lockdown rate z * i and the value of the efficiency is inversely proportional to each other. In other words, the higher the level of the allowed economic activities, the more effective the optimal lockdown policy (compare to the best uniform lockdown). Besides, it can be observed from SI Fig 6(a) , 6(f) that the recovery rate γ is the parameter that has the most significant impact on the value of z * i . In particular, when γ ranges from 0.03 to 0.1, the value of z * i for all the three models increases more than 0.5. Therefore, if we can make the patients recover faster, we can send the number of infections to 0 quickly while maintaining a high level of economic activity.γ decides the value of the recovery rate of symptomatic individuals in the COVID-19 model, hence we allow a higher level of economic activities when the value ofγ increases, this can be seen in SI Fig 6(e) . Moreover, the value of the initial growth rate is closely related to the value of the transmission rate β s , the larger the initial growth rate is, the larger β s will be for each model. As a consequence, we have to maintain a lower level of economic activities to make sure the infections go to 0 quickly (SI Fig 6(b) ).

It can be seen from SI Fig 6(d) and 6(c) that the value of z * i is not quite sensitive to the symptom rate andα. When the value ofα increases, the corresponding value of z * i slightly increases. The reason may be that the increasing ofα can lead to the increasing of λ max (M(t)). As we fix the value of the initial growth rate λ max (M(t 0 )), the obtained value of transmission rate β s can be smaller, which means we can allow a higher level of economic activities.

7.3. Other Parameter Analysis.. Except the disease parameters, the structure of the data may also have a great impact on the value of the generated optimal lockdown rate z * i . To clearly investigate these relationships, we depict the optimal z * i with respect to the centrality, the home-stay rate (i.e., daily homedwell-time in minutes/total number of minutes in a day), the population and the employment for each county 38 This manuscript is for review purposes only.

in SI Fig 7. It can be observed that the value of z * i for all three models increases when the home-stay rate increases. However, we do not observe any direct relationship between centrality, population, employment, and the value of z * i . To further verify this observation, we implemented random permutation experiments.

7.4. Random Permutation Analysis.. To fully understand the relationship of the value of z * i and the data structure, we implemented random permutation experiments with respect to degree, home-stay rate, population, employment and initial susceptible rate. In a random permutation experiment, we randomly permute one parameter while keep all the others fixed, then we fit the models to the randomly permuted data and apply the proposed algorithms get the optimal lockdown rate. We repeat this process 100 times for each parameter, and compare the average values of the lockdown rate on permuted data with the original data. The disease parameters are chosen as in [7] . In order to do a random permutation of the degree of the nodes, we fix the network, and randomly permute all the other parameters. We implemented such experiments based on data on different time points. The experimental results are reported in SI Fig 8 and SI Fig 9. It can be observed that the shapes of the histogram after randomly permuting the employment and initial susceptible rate are similar to the shape of the histogram of the original data. This implies the value of z * i does not have a close relationship with employment and susceptible rate. However, the shapes of the histogram after randomly permuting the degree, home-stay rate, and population are quite different from the original data. In other words, the centrality of the node, home-stay rate, and population have a close relationship with the value of z * i . To find out what kind of relationship between them, we implemented further experiments on synthetic data. 7.5. Impact of Centrality.. To study the impact of centrality, we generated some geometric random graphs and then added additional "hotspots" with high degree. This was done by randomly choosing several nodes in the initially generated geometric random graphs, and letting the edges leading from these nodes to other nodes be present with probability 0.9.

We let the remaining parameters be identical for all the nodes. Specifically, we let the population of all the nodes be 4000, the home-stay rate of all the nodes be 0.8, and we assume the number of employees of each node is proportional to the population. Since the SIS model and the SIR model will be identical if we let the initial susceptible rate of each node be the same constant (we choose β s to make the initial growth rate of each model be equal). Besides, according to our previous analysis, when the values of the initial susceptible rate of the nodes are close, it will not impact the value of the optimal lockdown rate of each node too much. As a consequence, we choose initial susceptible rate from interval [0.8, 0.9] uniformly at random. Moreover, according to the data [74] provided by SafeGraph, people tend to spend more time in their own counties rather than travel to other counties when they are not stay at home. To make the synthetic data similar to the real data, we update the definition of the travelling matrix τ as following:

where h i represents the home-stay rate of node i. In addition, the disease parameters of this experiment are set as in [7] , and the missing parameters are set as:α = 0.6754 (β a =αβ s ), r a = r s = γ. Similar to the experiments on real data, we let the decay rate α = 0.2γ = 0.04 to satisfy the assumption that α ≤ γ.

The results of this experiment are presented in Figure 4a -b. It can be observed that centrality only matters for the value of z * i when there exist hotspots (i.e., highly central nodes) in all the three models. Surprisingly, beyond such hotspots, the effect of centrality is essentially nonexistent.

To further study the impact of the centrality, we generated random graphs based on Barabási-Albert model, where there exist few nodes with unusually high degree compared to other nodes in the network. By 39 This manuscript is for review purposes only.

varying the the number of the nodes, the number of links and the initial seed for the B-A model, we can generate various random graphs. Besides the adjacency matrix A, all the other parameters about the data (population, τ, employment, initial susceptible rate, home-stay-rate), the disease parameters and the decay rate α are all chosen as in the experiment based on geometric graph. The experimental results are presented in SI Fig 10(a) . From SI Fig 10(a) , we can observe the similar phenomenon: the hotspots are assigned with smaller values of z * i by the optimal lockdown policy in all the three models, however, such effect of the centrality is nonexistent for all the other nodes.

Moreover, we also define a kind of random graphs to provide further evidence for the relationship of the centrality and the value of z * i . To generate graphs where the degree of each node can be different and can be decided by us, we define the adjacency matrix A of the random graph as following: the uppertriangular element A i j will be 1 with probability p i , and be 0 with probability 1 − p i , where p is the given probability vector; the values of the elements in the lower-triangular matrix are equal to the symmetry of the upper-triangular matrix. In other words, the edges leading from different nodes can be present with different probabilities. Similarly, the disease parameters, the decay rate α, and all the other parameters about the data are chosen as in the experiment based on geometric random graphs. The experimental results are presented in SI Fig 10(b) . From SI Fig 10(b) , we can also observe that centrality only matters for the value of z * i when there exist hotspots. Except the hotspots, the effect of centrality is centrality is nonexistent. 7.6. Impact of Population.. In the experiments about the population, we fix all other parameters and vary the population of each node from 100 to 20000. Since the employment is proportional to the population, and the economic cost c i of each node is decided by the corresponding number of employees, the economic cost of different nodes will also be different. To make sure the centrality of all the nodes are close, we experimented with random regular graph, Erdos-Renyi random graph, and the geometric random graph. All the other data parameters, disease parameters, and the decay rate α are set as in the experiments about the centrality. The experimental results are presented in Figure 4c -d, it can be observed that nodes with small populations are assigned with smaller values of z * i , surprisingly once the population is not very small, the effect is almost nonexistent. 7.7. Impact of Home-Stay Rate.. To study the impact of home-stay rate, we fixed all other model parameters and tuned the home-stay rate of the nodes. We chose random regular graph, Erdos-Renyi random graph, geometric random graph, and the 2d grid as our network as the centrality of the nodes on these graphs are similar. The disease parameters, the dacay rate α, and all the other data parameters are set as in the experiments about the centrality. The experimental results are presented in SI Fig 11. We found that z * l increases with increasing home-stay rate, which agrees well with our intuition.

8. City-Suburb Model. We now revisit the phenomenon we have observed in our analysis of NY, which is that the optimal lockdown tends to shutdown the outside of NYC harder than the NYC itself. To isolate this phenomenon in the simplest possible setting, we implement a simple synthetic experiment of a network with two nodes: node 1 will be referred to as "the city" while node 2 will be referred to as "the suburb." The city will have a large number of population while the suburb will have a smaller population. We will assume employment is proportional to the population so that the economic cost coefficient c i is also proportional to the population. Then we apply the proposed algorithms to design the optimal shutdown policy to this city-suburb model.

The disease parameters are taken from [7] where considers a SEIR model; we set the missing parameters asα = 0.6754 (where β a =αβ s ) and r a = r s = γ. We still choose the target decay rate α = 0.2r s = 0.04 as before. For the choice of the travel rate matrix τ, we choose a matrix that is similar to the New York State data, but with rounder numbers; specifically, we define τ i j = (1− h i 1440 )· k i j a k ia , where h i is the home-dwell-time 40 This manuscript is for review purposes only.

of node i, k ia is the number of trips from node i to node a, and we let Table 3 . We see the same phonemon as in our New York State experiments: the optimal lockdown policy shutdown the suburb more stringently than the city even though, in cases 1 and 2, the epidemic is mainly localized in the city. Comparing the results of Case 1 and Case 2, we see that this trend gets stronger when the population difference between the city and the suburb increases.

9. Beyond Eigenvalue Bounds. In this section, we further test our finding that the optimal stabilizing shutdown using New York State data is more stringent outside NYC. Our goal is to directly compare lockdowns by comparing the total number of infections. Unfortunately, we do not know of any way to optimize lockdowns efficiently to minimize the total number of infections; indeed, this difficulty is what motivates optimization of eigenvalue bounds in the first place.

We will compare the optimal stabilizing lockdown computed by our method, which we will say has cost c * , with lockdowns of the following structure:

• Two-parameters lockdown: shutdown the counties in NYC by z 1 , and shutdown the other counties in New York State by z 2 , where z 2 > z 1 , such that the economic cost of the shutdown is at most c * . • Uniform lockdown (cost matching): shutdown all the counties in New York State by z, such that the economic cost of the shutdown is at most c * . Because these lockdowns are characterized by, respectively, two and one parameters, they can be found via direct search (i.e., discretizing the parameters and trying every possibility). To summarize, we will perform an exhaustive search, not only over uniform lockdowns, but overall lockdowns defined by two z i , one inside NYC and one outside, with the former being smaller (corresponding to a more stringent lockdown).

We will apply the disease parameters from literature [8] , [7] , and [31] respectively. Similarly, we let the decay rate α = 0.2r s as before. The experimental results are presented in SI Fig 12. These findings further substantiate our results: our policy outperforms the best two-parameters lockdown as well as the best uniform lockdown in all scenarios.

We conclude with a cautionary tale about what happens when the eigenvalues are pushed too far into the left-hand plane. We have already remarked, in the main body of the paper, that this could result in poor performance as the improvement in asymptotic rate starts to come at the expense of transient performance; we next demonstrate how this phenomenon occurs in a series of charts. We experimented with decay rate α = 0.5r s , the experimental results are presented in SI Fig 13. Note that, in two of the three charts, this corresponds to an asymptotic rate with roughly 10% decrease in the number of cases per day. This extremely aggressive rate results in an optimal shutdown which, while attaining this rate, does not perform well. In particular, in the figures with a 10% decrease rate, we see that our optimal stabilizing shutdown underperforms the uniform shutdown and the two-parameters shutdown with the same cost.

Sensitivity analysis: To further study the impact of decay rate as well as various disease parameters to the total number of infections, we implement sensitivity experiments. In each experiment, we vary the value of one parameter while fix the values of the others. The normal values of γ, , and the initial growth rate are 41 This manuscript is for review purposes only. chosen as in [7] , and the normal value ofα (β a =αβ s ) is set as in [31] . We let the normal value of the decay rate α = 0.2r s = 0.04, since it is impossible to achieve a decay rate better than γ. For simplicity, we assume γ = r a , and r s =γr a , where the normal value ofγ = 1. The experimental results are presented in SI Fig 14  and 15 .

As expected, it can be observed from SI Fig 14(a), 14(d) , and 15(a) that the total number of the infections of our policy is much better for small α; but if we increase α too much, this starts to change and we no longer outperform. Additionally, it can be seen from SI Fig 14(c) , 14(f), and 15(b) the total number of infections of all the polices increases as the initial growth rate increases, again as expected; however, somewhat surprisingly, we can see from SI Fig 14(b), 14(e) , and 15(c) that the total number of infections increases when the recovery rate γ increases. The reason for this counter-intuitive phenomenon is that we assume the initial growth rate is fixed in this experiment, so that an increase in the recovery rate γ means an increase in the transmission rate β s to obtain the same initial growth rate, and the latter effect increases the total number of infections.

For COVID-19 model, besides α, γ, and the initial growth rate, we also studied the impact of parameters ,α, andγ. We can observe from SI Fig 15(d) that the total number of the infections is not sensitive to parameter for all the three polices. As the transmission rate of asymptomatic individuals β a =αβ s , then the total number of infections increases asα grows (SI Fig 15(e) ). Parameterγ is closely related to the recovery rate of the symptomatic individuals, its effect is similar to γ, when the value ofγ increases, the transmission rate β s increases and further lead to the increasing of the total number of the infections (SI Fig 15(f) ).

Outside NYC Harder? An Intuitive Explanation. Our main finding in the empirical experiments -that an optimal stabilizing lockdown will shut down the outside of NYC more stringently than NYC -is counter-intuitive. At least part of the reason why this happens is that our costs are taken to be proportional to employment: we do not attempt to take into account either the larger salaries of workers in NYC or the importance of NYC to the national economy. Nevertheless, the finding remains counter-intuitive, and we are not aware of any previous literature pointing out that something like this can occur.

We now provide an explanation which makes this finding more intuitive. Let us consider an idealized city-suburb model with two nodes. Imagine now that the coefficients c i are proportional to population, and imagine that, initially, there is no travel between city and suburb. For clarity, let us consider a related model where, instead of minimizing cost subject to a constraint on how fast infections decay, we instead put a price on each infection; this is closely related to what we do and slightly simpler to reason about. In that case, the optimal shutdown should be invariant to scaling in population since doubling the population of a city will both double the cost of a shutdown as well as double the benefit in terms of number of infections averted. Now let us consider further what happens when we change the system by stipulating that: (1) 1% of trips now happen between city and suburb (2) the population N of the city goes to infinity while the population of the suburb remains fixed at, say, 1000.

In this case, whenever we shut down the suburb, the cost scales with the population 1000; while the benefit also scales with N, the population of the city, since any shutdown in the suburb decreases the total number of infections in the city (since the epidemic can spread in the city through interactions with the suburb). Thus, even though only 1% of the trips go between city and suburb in each direction, the benefit of shutting down the suburb will overwhelm the cost as N → +∞. The same argument does not apply to the city: the increased benefits from reducing the number of infections in the suburb is small relative to the cost which scales with N. This provides an explanation why an optimal shutdown might choose to be more stringent in the suburb.

Of course, this is an idealized example. We note that in our empirical results, costs were taken to be proportional to employment, not population; so that we will never encounter a situation where the costs of 42 This manuscript is for review purposes only.

shutting down a suburb do not scale with the size of the city (because if 1% of the trips go between city and suburb, the suburb will also have employment that will scale with N to provide services to visitors). Nevertheless, our empirical results suggest that this idealized example is not too far from what happens once we use real data on number of trips and employment in NYC and suburbs, as the optimal stabilizing shutdown does choose to be more stringent outside NYC.

11. Extended Lockdown Cost Functions. In this section, we explore the optimal lockdown with additional cost functions. Our goal is to show the counter-intuitive phenomenon we observed earlier -that the optimal lockdown tends to shutdown the outside of NYC harder than the NYC itself -is not due to the particular choice of the economic cost functions in Eq. (2.6) in the main text. The choice of the new economic cost functions follow the similar rules in (2.6).

We thus consider other cost functions. First, we consider

By using the similar method as in SI Sec. 4.4, we can write the optimal lockdown issue with economic cost (11.1) as a constrained convex optimization problem as following

The constraints of (11.2) is convex SDP constraints, and the cost function is also convex. Hence, we can use the projected gradient descent (PGD) algorithm to obtain the optimal solution of this problem. Each step of the projection involves solving a semi-definite program.

Besides the cost (11.1), we also consider the following economic cost functions:

Similarly, we can also write the optimal lockdown issues with these costs as the convex optimization problem with convex SDP constraints, which can be addressed by the PGD algorithm.

To check if we can still observe the similar counter-intuitive phenomenon, we experimented with the data about COVID-19 break in NYS with these extended cost functions. In our simulations, we employed 43 This manuscript is for review purposes only. the toolbox from [54] to solve the projection step in each step of PGD problem. The disease parameters are set as in [31] . The decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). The experimental results are presented in SI Fig 16 and SI Fig 17. It can still observed that the value of z * i for counties in NYC are relatively higher than other counties in NYS in each scenario, which implies we should shutdown the outside of NYC harder than itself. Additionally, the median of z * i is greater than or close to the value of the uniform lockdown in all the cases. In other words, compared to the uniform lockdown, the optimal lockdown policy not only leads to less economic losses but also allows the majority of the counties to have more economic activity.

12. Robustness Check of Travel Rate Matrix τ . In this section, we check the robustness of one of our main findings, that the optimal lockdown policy should shutdown the outside of NYC harder than itself, with respect to the travel rate matrix τ. In our simulations of NY, we used the data from SafeGraph [74] to construct the travel rate matrix τ. The data was collected by using a panel of GPS pings from anonymous mobile devices, and might be inaccurate due to various practical reasons. For instance, we do not know from how many of the minutes recorded as travel outside might have been spent alone in a vehicle, with no possibility for transmission to other individuals. To test the robustness of our finding against the uncertainty of the travel rates, we introduce two types of perturbation to the travel rate matrix τ by either adding noise to its entries or removing a fraction of its entries. Then we check if the optimal lockdown policy still suggests that we should shutdown the outside of NYC harder. In both scenarios, we perturb the travel rates of each county while keep the home-stay-rate of each county be the same.

In the first scenario, we add Gaussian noises to each entry of the matrix τ. To do this, for each τ i j , we generate Gaussian noise g i j with mean 0 and variance θτ 2 i j , then we set

Next, we let

where β i is chosen to make j τ new i j = j τ i j . In this case, the size of the noises will be proportional to τ i j , and we can tune the size of the noises by varying the parameter θ. Meanwhile, recall the definition of the travel rate matrix (6.1), we still have j τ new i j = j τ i j = 1 − h i 1440 , which means the home-stay-rate of county i keeps the same.

In our second scenario, we remove a p-fraction of entries from each row of matrix k uniformly at random, and replace these values by 0, where k i j is the daily number of trips from county i to county j. This is to mimic missing data. We call the updated travel natrix as k new , then we generate τ with our original method,

Similarly, we can tune the size of the missing data by simply varying p. We implemented the two different perturbation scenarios with the data about COVID-19 break in NY. The disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. The simulation results are presented in SI Fig 18. It can be observed that the average of z * i for NYC cities is greater than the average of other counties for a wide range of the noise level or the fraction of missing data in the travel rate matrix. Thus our finding that the optimal lockdown should shutdown NYC harder is robust against considerable uncertainty in the travel rate matrix τ. 44 This manuscript is for review purposes only.

13. Nonuniform Transmission Rate β: Potential Urban-Rural Differentials. We now consider the possibility that the spread of an epidemic, as measured by the parameters β s and β a depends on the location. Specifically, we take a simple model where transmission is proportional to a power h of the population density of each county.

There are a number of reasons why that might be the case. In an epidemic where the primary mode of transmission is through outside contact, the transmission speed will naturally scale with population density. Additionally, higher density counties will on average have more extensive public transport, which could also facilitate transmission even in the presence of masking guidelines. However, for COVID-19 in New York State, it is not clear that COVID-19 transmission substantially increases with density.

The data can be seen in SI Fig 20. The bottom two graphs of that figure show population density vs R t (average and range). We do not in general see an upward curve. The first graph of that Figure presents a plot of R t for New York Country (approximately 70, 000 people per square mile) and Hamilton County (approximately 2.6 per square mile); these are, respectively, the most and least dense counties in New York. Although R t fluctuates over time for both counties, SI Fig 20 does not show any persistent differences, in spite of a ×30, 000 density differential. This rules out the possibility of a strong density dependence, e.g., β being proportional to density, though it leaves open the possibility that population density could have some slighter influence which is counteracted by the differences in other characteristics between these counties.

We thus consider what happens if we introduce some scaling with a power h of population density into our models; we stress that h here could be quite close to zero. It is clear that, as we increase h, at some point our finding that NYC should be shut down less harshly than the rest of NYS will reverse: as infections comparatively spread faster and faster in NYC, the NYC shutdown should get more stringent, at some point becoming more stringent that for the remainder of the state. We want to see how big h has to be for this reversal to occur. As before, will find it helpful to frame this h in terms of the differential between the most and least populated counties. Thus we seek to answer the following question: how much faster does the infection need to spread in New York County relative to Hamilton County for this reversal to occur?

More formally, suppose p l is the population density of location l. For COVID-19 model, we let the transmission rates associated with location l be

where, recall,α is a fixed scalar which measures the transmission rate ratio between symptomatic individuals and asymptomatic individuals. The parameter k is chosen to match the initial growth rate of the available COVID-19 models, just as before. In this case, we can write A 0 in Eq. (4.10) as

We can then apply our approach to get the optimal stabilizing lockdown. We do this, getting the population density data from [35] . The simulation results are presented in SI Fig 19. It can be observed that the value z * i of NYC counties still larger than the other counties when h ≤ 0.1; however, this phenomenon becomes unintelligible by h ≥ 0.15. Coming back to the New York County/Hamilton County divide, we find that our results still hold if β corresponding to the former is roughly 2.8 times the β corresponding to the latter.

14. Impact of The Cases Decline Speed. In this section, we consider the cases where the lockdown sends the confirmed cases to zero with different rates. Our goal is to show the counter-intuitive phenomenon we observed earlier-that the optimal lockdown tends to shutdown the outside of NYC harder than NYC 45 This manuscript is for review purposes only.

itself -is not due to the particular choice of the decay rate α. In particular, we consider the cases when the confirmed cases are reduced very fast, with 10 % (or even 20%) daily decline speed.

Suppose we aim to reduce the confirmed cases with a speed at least a * through lockdown, we will compare the optimal stabilizing lockdown computed by our method with the following two simple lockdown policies:

• Best two-parameters (NYC) lockdown: shutdown the counties in NYC by z 1 , and shutdown the other counties in New York State by z 2 , where z 2 > z 1 , such that the daily decline speed of the COVID-19 model is at least a * and the corresponding economic cost is minimal. • Best two-parameters (outside) lockdown: shutdown the counties in NYC by z 1 , and shutdown the other counties in New York State by z 2 , where z 1 < z 2 , such that the daily decline speed of the COVID-19 model is at least a * and the corresponding economic cost is minimal. Since both of the two lockdown policies are characterized by two parameters, we will find them via grid search (i.e., discretizing the parameters and trying every possibility). To summarize, we will fix the decline speed of the epidemics and compare the economic cost of different lockdown policies.

Similar to before, we will apply the disease parameters from literature [7, 31] and [8] respectively. Since we need to ensure α < min(r a , r s ), and the values of r a , r s are different in these literatures, the range of a * we consider are different for different groups of parameters. The experimental results are presented in SI Fig  21. It can be observed, when the decline speed is fixed, the best two-parameters (NYC) lockdown policy leads to more economic losses than the best two-parameterss (outside) lockdown and our method, even when the cases decline very fast. This finding provides further support for the counter-intuitive phenomenon we observed that the optimal lockdown tends to shutdown the outside of NYC harder.

15. Robustness Check: Activity Differential between the Symptomatic and Asymptomatic Individuals. In this section, we consider what happens if symptomatic and asymptomatic individuals have different activity levels. We consider a simple model where the travel rate of symptomatic individuals is a constant fraction of the asymptomatic individuals. This could happen because some symptomatic individuals may choose to isolate or quarantine themselves.

If the activity level of symptomatic individuals is zero, the COVID-19 model we have proposed here reduces to the standard SIR model, as people in the symptomatic class will not infect anyone. In general, however, the activity level of symptomatic individuals may not equal zero: some people may ignore regulations, others may ignore their symptoms, and still others may believe they are suffering from something other than COVID. Because we are not aware of any research allowing us to choose a specific activity reduction, we will consider a number of possibilities for how much symptomatic individuals reduce activity relative to asymptomatic individuals, ranging from 10% to 90%.

Specifically, we set τ s i j = κτ a i j . Equivalently, we may multiply our coefficientα with 1 κ and obtain the same model, because, recall, α = β a β s , and we choose the value of transmission rate β s to match the given initial growth rate. In our simulations, we experimented with k = 0.1, 0.5, 0.9, 1.0 respectively. The simulation results are presented in SI Fig  22. We find that the optimal stabilizing lockdowns z * i vary when κ takes on different values, however, the patterns of z * i are quite similar. Specifically, we still can observe that the optimal stabilizing lockdown tends to shutdown the outside of NYC harder than itself. 46 This manuscript is for review purposes only. 

. . . . . . . . . . . .

:-, : : . d-f, the estimated cumulative cases in NY from Apr. 1st, 2020 to Aug. 14th, 2021 (or May 10th, 2024). In a, d, the disease parameters are set as in [7] , the decay rate α is chosen as 0.0231 which corresponds to halving every 30 days. In b, e, the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In c, f, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. Uniform lockdown, random lockdown, and uniformly-bouded-decline lockdown are defined as in Fig 1. It can be observed that our policy outperforms all the other lockdown polices.

This manuscript is for review purposes only. . . d-f, the estimated cumulative cases in NY from Apr. 1st, 2020 to Aug. 14th, 2021 (or May 10th, 2024). In a, d, the disease parameters are set as in [7] , the decay rate α is chosen as 0.0231 which corresponds to halving every 30 days. In b, e, the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In c, f, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. Uniform lockdown, random lockdown, and uniformly-bouded-decline lockdown are defined as in Fig 1. It can be observed that our policy outperforms all the other lockdown polices. Experimental results of real data: lockdown rate of each county given by different polices for SIR model based on available data about COVID-19 outbreak in NY on April 1st, 2020. a-c, optimal lockdown rate z * i given by our method . d-f, uniform lockdown rate z i . g-i, random lockdown rate z i . j-l, uniformly-bounded-decline lockdown rate z i . Uniform lockdown, random lockdown, and uniformly-bouded-decline lockdown are defined as in Fig 1. In a, d, g, j, the disease parameters are set as in [7] , the decay rate α is chosen as 0.0231 which corresponds to halving every 30 days. In b, e, h, k the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ).

In c, f, i, l, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. It can observed from a-c that the value of z * i for counties in NYC are relatively higher than other counties in New York State, which implies we should shutdown the outside of NYC harder than itself. Besides, it can be seen that such counter-intuitive phenomenon does not appear in any other lockdown polices. 50 This manuscript is for review purposes only. [7] . The decay rate is chosen as α = 0.2r s = 0.04 so that α < min(r a , r s ). It can be seen that the value of z * i tends to increase from March to June, as people travel less frequently due the impact of COVID-19. 51 This manuscript is for review purposes only. Fig. 6 Disease parameter sensitivity analysis with respect to z * i and efficiency: optimal lockdown rate z * i and corresponding efficiency given by Theorem 4.6 for SIS, SIR, and COVID-19 model based on available data about COVID-19 outbreak in New York State on April 1st. In a-e, the shaded region represents the standard deviation of z * i for all 62 counties, the solid line is the mean value of z * i for all 62 counties. Parameter ,α,γ only appears in COVID-19 model, hence c-e and h-j only show the results of COVID-19 model. It can be seen that the value of z * i and the corresponding economic cost are sensitive to recovery rate and the initial growth rate but not to other parameters. Fig. 7 Other parameter analysis: the relationship of the optimal lockdown rate z i and the centrality, home-stay rate, population and employment. z * i are produced by Theorem 4.6 based on available data about COVID-19 outbreak about COVID-19 outbreak in New York State on April 1st, 2020. The disease parameters are set as in [7] . In a, the disease parameters are set as in [7] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. In b, the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In c, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. Degree closely related to centrality, and the number of employees decides the economic cost coefficients c i . Each point represents a county in New York State. It can be observed that the value of z * i increases as the home-stay rate of the node increases, the impacts of the degree, population and the employment to the value of z * i are not obvious. 53 This manuscript is for review purposes only. Fig. 8 Experimental results of random permutation: optimal lockdown rate z * i after doing random permutation with respect to degree, home-stay rate, population, employment, and initial susceptible rate. a, shows the results based on available data about COVID-19 outbreak in NYS on March 2nd, 2020. b, the results based on available data about COVID-19 outbreak in NYS on April 1st, 2020. The decay rate is chosen as α = 0.2r s = 0.04 so that α < min(r a , r s ). The results are average of 100 repeat. It can be observed that the distribution of z * i can be very strongly affected by permutations of centrality, population, and the home stay rate. On the other hand, the distribution of z * i is not altered much by permuting employment and the initial susceptible rate. 54 This manuscript is for review purposes only. (b) Fig. 9 Experimental results of random permutation: optimal lockdown rate z * i after doing random permutation with respect to degree, home-stay rate, population, employment, and initial susceptible rate. a, the results based on available data about COVID-19 outbreak in NYS on May 1st, 2020. b, the results based on available data about COVID-19 outbreak in NYS on June 2nd, 2020. The disease parameters are set as in [7] . The decay rate is chosen as α = 0.2r s = 0.04 so that α < min(r a , r s ).The results are average of 100 repeat. It can be observed that the distribution of z * i can be very strongly affected by permutations of centrality, population, and the home stay rate. On the other hand, the distribution of z * i is not altered much by permuting employment and the initial susceptible rate. 55 This manuscript is for review purposes only. Fig. 10 Experimental results of synthetic data: the impact of centrality of the network to the value of optimal lockdown rate z * i . a, the results based on Barabási-Albert model. b, the results based on a kind of generated random graph. Each point represents a node in the network. It can be observed that centrality only matters for the value of z * i when there exist hotspots in all the three models. Surprisingly, beyond such hotspots, the effect of centrality is essentially nonexistent. 56 This manuscript is for review purposes only. Fig. 11 Experimental results of synthetic data: the impact of home-stay rate to the value of optimal lockdown rate z * i . Each point represents a node in the network. It can be observed that z * i increases as the home-stay rate increases, and in fact the home-stay rate has by far the biggest influence on z * i compared to the other parameters we consider. 57 This manuscript is for review purposes only. (c) Fig. 12 Experimental results of real data: the estimated rate of accumulative cases in the total population by applying different lockdown policies based on available data about COVID-19 outbreak in New York State on April 1st, 2020. It can be observed our policy outperforms the best two-parameters lockdown as well as the uniform lockdown in all the three models. In a, the disease parameters are set as in [7] , and the decay rate α = 0.2r s = 0.04. In b, the disease parameters set as in [31] , and the decay rate α = 0.2r s = 0.0034. In c, the disease parameters set as in [8] , and the decay rate α = 0.2r s = 0.058. 58 This manuscript is for review purposes only. (c) Fig. 13 Experimental results of real data: the estimated rate of accumulative cases in the total population by applying different lockdown policies based on available data about COVID-19 outbreak in New York State on April 1st, 2020. In a, the disease parameters set as in [7] , and the decay rate α = 0.5r s = 0.1. In b, the disease parameters set as in [31] , and the decay rate α = 0.5r s = 0.0085. In c, the disease parameters set as in [8] , and the decay rate α = 0.5r s = 0.145. In b, it can be seen that our policy outperforms the best two-parameters lockdown as well as the uniform lockdown in SIS model and SIR model. In a, c, and the COVID-19 model in b, it can be seen that our policy underperforms the best two-parameters lockdown as well as the uniform lockdown, the reason for such phenomenon is provided in SI Sec. 9. 59 This manuscript is for review purposes only. 14 Sensitivity analysis of disease parameters with respect to final cumulative cases: the estimated rate of final accumulative cases (at day 500) in the total population for SIS model and SIR model, respectively based on available data about COVID-19 outbreak in New York State on April 1st, 2020. It can be observed from a and d that the total number of the infections of our policy is much better for small α; but if we increase α too much, this starts to change and we no longer outperform. Besides, from c and f we can see that the total number of infections of all the polices increases as the initial growth rate increases. In addition, it can be observed from b and e that the total number of infections increases when the recovery rate γ increases, the reason for this counter-intuitive phenomenon is provided in SI Sec. 9. 66 This manuscript is for review purposes only. Economic cost best two-parameters (NYC) best two-parameters (outside) ours a b c Fig. 21 The economic cost of different lockdown policies when the daily decline speed of the confirmed cases is fixed. In a, the disease parameters are set as in [7] , the daily cases decline speed ranges between [1%, 15%]. In b, the disease parameters set as in [31] , the daily cases decline speed ranges between [0.01%, 0.3%]. In c, the disease parameters set as in [8] , the daily cases decline speed ranges between [2%, 24%]. The best two-parameters (NYC) lockdown and the best wo-parameters (outside) lockdown are defined in SI Sec 14. [31] 0.034 0.034 0.017 0.125 0.6754 [7] 0.20 --0.32 - Table 3 Optimal lockdown rate z * i for the City-suburb model. The first node is the city and the second one is the suburb. It can be observed that the optimal lockdown policy shutdown the suburb more stringently than the city even though, in cases 1 and 2, the epidemic is mainly localized in the city.

Case 1 Case 2 Case 3

Optimal targeted lockdowns in a multi-group SIR model

A simple planning problem for COVID-19 lockdown

A multi-city epidemic model

A deterministic model in biomathematics. asymptotic behavior and threshold conditions

Measured voluntary avoidance behaviour during the 2009 A/H1N1 epidemic

Nonnegative matrices in the mathematical sciences

The challenges of modeling and forecasting the spread of COVID-19

Controlling epidemic spread: Reducing economic losses with targeted closures. University of Chicago

Optimal control and basic reproduction numbers for a compartmental spatial multipatch dengue model

Time-optimal control strategies in sir epidemic models

A multi-layer network model to assess school opening policies during the covid-19 vaccination campaign

Optimal design of lock-down and reopening policies for early-stage epidemics through sir-d models

A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2

Lectures on network systems

Population -NYS & Counties

This manuscript is for review purposes only

Applying optimal control theory to complex epidemiological models to inform real-world disease management

Model predictive control to mitigate the covid-19 outbreak in a multi-region scenario

More than nine-in-ten people worldwide live in countries with travel restrictions amid COVID-19

Mobility network models of covid-19 explain inequities and inform reopening

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

COVID Reproduction Rate

Matrix scaling and balancing via box constrained Newton's method and interior point methods

Covid-19 dashboard

A network model of italy shows that intermittent regional strategies can alleviate the covid-19 epidemic

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

Optimal lockdown in a commuting network

Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe

Matrix theory. Courier Corporation

Applications of the Theory of Matrices. Courier Corporation

Mitigation strategies for pandemic influenza in the United States

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Jue insight: How much does covid-19 increase with mobility? evidence from new york and four other us cities

Control systems: principles and design. Tata McGraw-Hill Education

New York State Statewide COVID-19 Testing

Population, Land Area, and Population Density by County

Epidemic spreading and control strategies in spatial modular network

Global stability of the endemic equilibrium of multigroup SIR epidemic models

Association of state-issued mask mandates and allowing on-premises restaurant dining with countylevel covid-19 case and death growth rates?united states

Differentiability of the hartman-grobman linearization

The type-reproduction number T in models for infectious disease control

Estimating the fraction of unreported infections in epidemics with a known epicenter: an application to COVID-19

A review of matrix scaling and sinkhorn's normal form for matrices and positive maps

Positive semidefinite programming: mixed, parallel

Covid-19 dashboard

Population flow drives spatio-temporal distribution of COVID-19 in China

On the complexity of matrix balancing

On the optimal control of virus spread in networks

Stability of epidemic models over directed graphs: A positive systems approach

Viral dynamics of sars-cov-2 infection and the predictive value of repeat testing. medRxiv

An optimal control problem for a spatiotemporal sir model

A deterministic model for gonorrhea in a nonhomogeneous population

The control handbook

Optimal control for information diffusion over heterogeneous networks

Yalmip : A toolbox for modeling and optimization in matlab

A mathematical model for predicting the geographic spread of new infectious agents

Convergence of the forward-backward sweep method in optimal control

On the dynamics of deterministic epidemic propagation over networks

Estimation of the asymptomatic ratio of novel coronavirus infections (COVID-19)

Analysis and control of epidemics: A survey of spreading processes on complex networks

Optimal resource allocation for control of networked epidemic models

Employment situation summary

Adaptive susceptibility and heterogeneity in contagion models on networks

Association of public health interventions with the epidemiology of the COVID-19 outbreak in Wuhan

The complexity of the matrix eigenproblem

Virus spread over networks: Modeling, analysis, and control

Modelling and predicting the effect of social distancing and travel restrictions on covid-19 spreading

Random geometric graphs

Optimal resource allocation for network protection against spreading processes

Optimal resource allocation for network protection against spreading processes

Distributed control of positive systems

The optimal control of infectious diseases via prevention and treatment. CEPR Discussion Paper

Optimal control of epidemics in metapopulations

This manuscript is for review purposes only

A structured epidemic model incorporating geographic mobility among regions

Concerning nonnegative matrices and doubly stochastic matrices

Convex optimization of the basic reproduction number

Coronavirus (Covid-19) Data in the United States

COVID-19 live updates: U.s. hospitalizations top 61

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Synchrony, waves, and spatial hierarchies in the spread of influenza

Analysis, prediction, and control of epidemics: A survey from scalar to dynamic network models

On assessing control actions for epidemic models on temporal networks

This manuscript is for review purposes only. Experimental results of real data: lockdown rate of each county given by different polices for SIS model based on available data about COVID-19 outbreak in NY on April 1st, 2020. a-c, optimal lockdown rate z * i given by our method . d-f, uniform lockdown rate z i . g-i, random lockdown rate z i . j-l, uniformly-bounded-decline lockdown rate z i . Uniform lockdown, random lockdown, and uniformly-bouded-decline lockdown are defined as in Fig 1. In a, d, g, j, the disease parameters are set as in [7] , the decay rate α is chosen as 0.0231 which corresponds to halving every 30 days. In b, e, h, k the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ).In c, f, i, l, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. It can observed from a-c that the value of z * i for counties in NYC are relatively higher than other counties in New York State, which implies we should shutdown the outside of NYC harder than itself. Besides, it can be seen that such counter-intuitive phenomenon does not appear in any other lockdown polices. 49 This manuscript is for review purposes only. This manuscript is for review purposes only. Fig. 15 Sensitivity analysis of disease parameters with respect to final cumulative cases: the estimated rate of final accumulative cases (at day 500) in the total population for COVID-19 model based on available data about COVID-19 outbreak in New York State on April 1st, 2020. It can be observed from a that the total number of the infections of our policy is much better for small α; but if we increase α too much, this starts to change and we no longer outperform. Besides, from b, we can see that the total number of infections of all the polices increases as the initial growth rate increases. In addition, it can be observed from c that the total number of infections increases when the recovery rate γ increases, the reason for this counter-intuitive phenomenon is provided in SI Sec. 9. The analysis for d-f can also be found in SI Sec. 9. 61 This manuscript is for review purposes only. Experimental results with extended lockdown cost function: optimal lockdown rate z * i given by Theorem 4.6 for SIS, SIR, and COVID-19 model based on available data about COVID-19 outbreak in New York State on April 1st, 2020. The disease parameters are set as in [31] . The decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In a, the cost function is chosen as i c i ( 1

. It can observed that the value of z * i for counties in NYC are relatively higher than other counties in New York State, which implies we should shutdown the outside of NYC harder than itself. Besides, the median of z * i is greater than the value of the uniform lockdown in all the cases. 62 This manuscript is for review purposes only. Fig. 17 Experimental results with extended lockdown cost function: optimal lockdown rate z * i given by Theorem 4.6 for SIS, SIR, and COVID-19 model based on available data about COVID-19 outbreak in New York State on April 1st, 2020. The disease parameters are set as in [31] . The decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In a, the cost function is chosen as i c i (min( 1 z i , 10) − 1). In b, the cost function is chosen as i c i (min( 1 z i , 20) − 1). In c, the cost function is chosen as i c i (min( 1 z i , 100) − 1). It can observed that the value of z * i for counties in NYC are relatively higher than other counties in New York State, which implies we should shutdown the outside of NYC harder than itself. Besides, the median of z * i is greater than or close to the value of the uniform lockdown in all the cases. 63 This manuscript is for review purposes only. [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. Solid lines show the average of 50 runs, the shaded regions show the min-max interval of the 50 runs. It can be observed that the average of z * i for NYC cities is greater than the average of other counties for a wide range of the noise level or the fraction of missing data. This suggests that our finding that the optimal lockdown should shutdown NYC harder is robust against the uncertainty of the travel rate matrix τ. 64 This manuscript is for review purposes only. Fig. 19 Optimal lockdown rate of each county calculated by our method for COVID-19 model when the transmission rate β is nonuniform. In the first row, we assume the transmission rate β at each county are the same. In the second row, we assume the transmission rate β l is proportional to p 0.05 l , where p l represents the population density at county l. In the third row, we assume the transmission rate β l is proportional to p 0.1 l . In the last row, we assume the transmission rate β l is proportional to p 0.15 l . In the first column, the disease parameters are set as in [7] , the decay rate α is chosen as 0.0231 which corresponds to halving every 30 days. In the second column, the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In the third column, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. The data used is about COVID-19 outbreak in NY on April 1st, 2020. 65 67 This manuscript is for review purposes only. Fig. 22 Optimal lockdown rate of each county calculated by our method for COVID-19 model when the activity level of symptomatic and asmptomatic individuals are different. In a-c, the travel rate τ i j of symptomatic individuals is 10% of the asymptomatic individuals. In d-f, the travel rate τ i j of symptomatic individuals is 50% of the asymptomatic individuals. In g-i, the travel rate τ i j of symptomatic individuals is 90% of the asymptomatic individuals. In j-l, the travel rate τ i j of symptomatic individuals are equal to the asymptomatic individuals. In the first column, the disease parameters are set as in [7] , the decay rate α is chosen as 0.0231 which corresponds to halving every 30 days. In the second column, the disease parameters are set as in [31] , the decay rate is chosen as α = 0.2r s = 0.0034 so that α < min(r a , r s ). In the third column, the disease parameters are set as in [8] , the decay rate α is chosen 0.0231 that corresponds to halving every 30 days. The data used is about COVID-19 outbreak in NY on April 1st, 2020. 68