## Abstract

The rapid rollout of the COVID-19 vaccine global raises the question of whether and when the ongoing pandemic could be eliminated with vaccination and non-pharmaceutical interventions (NPIs). Despite advances in the impact of NPIs and the conceptual belief that NPIs and vaccination control COVID-19 infections, we lack evidence to employ control theory in real-world social human dynamics in the context of disease spreading. We bridge the gap by developing a new analytical framework that treats COVID-19 as a feedback control system with the NPIs and vaccination as the controllers and a computational and mathematical model that maps human social behaviors to input signals. This approach enables us to effectively predict the epidemic spreading in 381 Metropolitan statistical areas (MSAs) in the US by learning our model parameters utilizing the time series NPIs (i.e., the stay-at-home order, face-mask wearing, and testing) data. This model allows us to optimally identify three NPIs to predict infections actually in 381 MSAs and avoid overfitting. Our numerical results universally demonstrate our approach’s excellent predictive power with *R*^{2} > 0.9 of all the MSAs regardless of their sizes, locations, and demographic status. Our methodology allows us to estimate the needed vaccine coverage and NPIs for achieving *R*_{e} to the manageable level and the required days for disease elimination at each location. Our analytical results provide insights into the debates on the aims for eliminating COVID-19. NPIs, if tailored to the MSAs, can drive the pandemic to an easily containable level and suppress future recurrences of epidemic cycles.

The ongoing global pandemic of coronavirus disease 2019 (COVID-19) has caused devastating loss of human lives and inflicted severe economic burden in the US^{1}. By 20 March 2021, more than 510,000 people were killed by COVID-19, the unemployment reaches 11.5%^{2}, and fiscal shortfall reaches over 200 billion^{3}. In response to the disease, federal, state, and local governments have implemented non-pharmaceutical interventions (NPIs) from encouragement and recommendations to full-on regulation and sanctions^{4–6} and pharmaceutical interventions (PIs) with available vaccinations^{7, 8} and drugs^{9}. However, these NPIs (e.g., social distancing, face mask-wearing, hand hygiene, testing, contact tracing, isolation, etc.) are often loosened and re-tightened without rigorous empirical evidence. Methodologies, from randomized controlled trials^{5}, econometric methods^{10, 11} to mathematical models^{12–18}, have been developed to measure the effects of these NPIs. But many of the methodologies cannot be used without predicting the pandemic when the NPIs are adjusted. Besides, different sets of NPIs enforced in different places at different times often come with different effects^{10, 16}. Without considering NPIs’ varieties on space, timing, and duration, we cannot understand whether these NPIs have had the desired effect of controlling the epidemic. In this study, we aim to tackle the challenge by modeling the COVID-19 spreading as a feedback control system, where the NPIs and vaccination note the controllers, and then to help policymakers determine the magnitude and timing of interventions’ deployment as circumstances change.

Engineering perspectives are useful in epidemic modeling^{19}, including the principle of control theory that provides a theoretical basis for NPIs’ functioning and transmission management. Control theory, originally developed for engineered systems with applications to power grids, manufacture, aircraft, satellite, and robots, has recently been adapted to understand the controllability of complex systems emerging in ecology, biology, and society. The recent work about network control enables us to identify the minimal driver nodes^{20} or lowest control costs^{19, 21} for node control, edge control^{22}, target control^{23}, multilayer control^{24, 25}, temporal control^{26}, and data-driven control^{27}. However, we continue to lack general answers to practically applying control theory to human and natural systems, like the dynamical system for disease spreading. The difficulty is rooted in the fact that we know the pedals and the steering wheel are the drivers prompting a car to move with the desired speed and in the desired direction, but the practical drivers are unknown for complex human and natural systems^{28, 29}. Specifically, when focusing on the COVID-19 pandemic, we hypothesize that non-pharmaceutical and pharmaceutical interventions are the “drivers” to determine the dynamics of the COVID-19 pandemic in each location through controlling the infection, recovery, and death rates^{29}. We validate this hypothesis by developing a parsimonious model that excellently predicts how the interventions influence the spreads in 381 MSAs in the US and ultimately estimates the end of COVID-19.

## Results

### Model COVID-19 spreading as a feedback control system

Increasing evidence shows that the spread of COVID-19 follows compartmental models^{10, 12, 30–32}, such as the SIRD (Susceptible-Infectious-Recovered-Death) model, which is mathematically described by the nonlinear equations expressing a population balance as follows^{33}:
where *S, I, R*, and *D* are the susceptible, infected, recovered and death numbers, respectively. *S* + *I* + *R* + *D* = Ω and Ω is the population at the given place. *V* is the ratio of people full vaccinated with efficacy rate 90%^{34}. The parameter *µ* is the crude birth and death rate, and epidemiological parameters *β, γ*, and *δ* are the infection, recovery and death rates (Θ = [*β, γ, δ*]^{T}). Studies show that epidemiological parameter are time-dependent, which adapt accordingly to the change of interventions^{12, 13}. By defining *β*_{0}, *γ*_{0}, and *δ*_{0} as the infection, recovery and death rates before interventions are imposed (Θ_{0} = [*β*_{0}, *γ*_{0}, *δ*_{0}]^{T}), we define Θ(*t*) = Θ_{0} + *U*_{Θ}(*t*), where *U*_{Θ}(*t*) = [*U*_{β}(*t*), *U*_{γ}(*t*), *U*_{δ}(*t*)]^{T} is a vector of the control input signals. As shown in Fig. 1a, the controllers *U*_{Θ}(*t*) work as the edge control^{22} and the vaccination *V* works as the node control^{20}. Based on nonlinear feedback control law, we develop the controllers *U*_{Θ}(*t*) using the feedback linearization approach. Then the SIRD model could perfectly track the reference trajectories which is generated from reported real-world pandemic data. As illustrated in Fig. 1b, the output model trajectory fits the real-world 3-dimensional data when the time-varying controllers are included, while it fails to fit when the controllers are not considered. Note that our approach is general and can be extended to other models that consider *n*-dimensional data.

Next, we propose the parsimonious model by employing the *difference-in-difference* estimation, which maps human behavior toward NPIs into the designed controller, as shown in Fig. 1c. This method enables us to measure the effect of NPIs on the controllers *U*_{Θ} by comparing the infection dynamics before and after the same region’s NPIs deployment. Beyond existing models, which assume the effects on policies are approximately linear^{10}, we also identify the interactions between policies. Thus, we are able to compile the evolution of human behaviour toward the NPIs, *θ*_{I} = [*θ*_{s}, *θ*_{f}, *θ*_{g}]^{T}, (e.g., stay-at-home order *θ*_{s}, face-mask wearing *θ*_{f}, and testing *θ*_{g}) to the evolution of the control signals with their respective effects, .

In short, our two-steps approach captures how the NPIs and vaccinations govern the disease dynamics, through combining the principles of control theory^{29, 35} and statistical estimation^{10, 36}.

### Tracking the pandemic’s trajectory with designed controllers

To validate its accuracy in predicting the disease evolution in the future and its applicability in determining the needed interventions to control the infection process, we test the NPIs data and infection data at 381 Metropolitan statistical areas (MSAs) from 1 April 2020 to 20 February 2021. The MSAs, defined as the core areas integrating social and economic adjacent counties, is contiguous areas of relatively high population and traffic density. We consider each MSA as a “closed population” the disease evolution could be modeled by the nonlinear SIRD dynamical model, Eq. (1), in a compact form, as
with *X* = [*S, I, R, D*]^{T} is the state/output vector and *U*_{Θ} = [*U*_{β}, *U*_{γ}, *U*_{δ}]^{T} is the input vector. Using feedback linearization control design, we construct a nonlinear feedback control law , to track the real-world pandemic trajectories *X*_{d} = [*S*_{d}, *I*_{d}, *R*_{d}, *D*_{d}]^{T} (see Eqs. (6-11) in Methods). The feedback law *ϕ*(.) relies on the measurements of the model full state *X* and the reference trajectory (*X*_{d}) and their time derivative ( and ).

Theoretically, for each MSA, the output trajectory *X* perfectly fits the reported infection trajectory *X*_{d} when the control action *U*_{Θ} is applied. This fact is illustrated in Fig. 2, which shows the evolution of the predicted and reported trajectories from 1 April 2020 to 20 February 2021 for two MSAs, namely the “New York” MSA (with New York as the core) and ‘Houston’ MSA (with Houston as the core). Moreover, the excellent alignments between the real-world and model data of all 381 MSAs indicate our approach’s predictive power, as shown in Fig. 2b. On the other side, Fig. 2c shows all MSAs’ infection rate [*β*_{0} + *U*_{β}(*t*)], recovery rate [*γ*_{0} + *U*_{γ}(*t*)], and death rate [*δ* + *U*_{δ}(*t*)], respectively, and mark out the respective rates for the two examples of MSA. Overall speaking, the infection rate and recovery rate decrease till May, rebound in June, decrease again to the lowest in September and start to fluctuate till February. The death rate continues to decrease in October and stays relatively constant beyond. As validated in literature^{4}, “New York” MSA enforces interventions more effectively and earlier, leading to a relatively lower infection rate, recovery rate, and death rate than most MSAs till October. One can observe that “Houston” MSA follows medium patterns.

Having the daily transmission rate, recovery rate, and death rate, we compute the effective reproductive ratio *;*, representing the expected number of secondary infected cases at time *t* when no vaccinations are rolled out. The function *R*_{e}(*t*) is a critical thresh-old in understanding whether the outbreak is under control. More precisely, if *R*_{e}(*t*) *<* 1, then the ongoing outbreak will eventually fade out, whereas *R*_{e}(*t*) > 1 means an acceleration of the infection dynamics leading to a substantial growth of infected cases and deaths. Fig. 2d-e respectively visualize the temporal and spatial distribution of the MSAs’ effective reproductive ratio. Our results reveal that just a few MSAs’ effective reproductive ratios ever reach *R*_{e}(*t*) *<* 1, implying the need to implement more rigorous interventions to achieve an effective controlling of the pandemic. For example, from 14 February 2021 to 20 February 2021, the average effective reproductive ratios in “New York” and “Houston” MSAs are evaluated as 1.468 and 1.393, respectively, revealing a critical need for stronger interventions.

### Mapping human behaviours as actuators to steer NPIs as controllers

Although our controller, *U*_{Θ}, precisely predicts the infectious, death toll, and recovery, it is thus far unknown how the measurable interventions change the controllers. Compared with the driving car, it is similar that we know the desired speed and direction in order to move the car to the target location, but we do not know what the pedals and the steering wheel for epidemic control are and how the changes in the pedals and the steering wheel determine the speed and direction. In another word, as depicted in Fig. 1, we assume that the changes of interventions, directly steering the control signals, will shape the disease dynamics when they are tightened and loosen to different levels. For multiple NPIs *θ*_{I}, we divide them into two sets, i.e., the set of community NPIs *ϑ*^{c} (e.g., social distancing and quarantine) and the set personal NPIs *ϑ*^{p} (e.g., face covering, test, and frequent hand wash). Then we develop the following mathematical model as (see [Eqs. (12-15)] in Methods)
where *Û*_{Θ}(*t*) is the estimate of the control action based on NPIs with their magnitudes *θ*_{j}(*t*) and their impact value . For a the specific NPI *j*, large value demonstrates that NPI *j* has strong impact. If , the control *Û*_{Θ}(*t*) is independent with the NPI *j*. The term Π _{i∈{c,p}} indicates that community NPIs *ϑ*^{c} and personal NPIs *ϑ*^{p} have a joint affect on the controller. For either community NPIs *ϑ*^{c} or personal NPIs term denotes the combined impact of the NPI set. Usually, *Û*_{Θ} is a non-positive value, showing the reductions in each controller.

Given the collected data sets for eight NPIs, we evaluate the goodness of fit versus different combinations of NPIs to find an effective parsimonious model for Eq. 3. A parsimonious model, which has great explanatory predictive power, accomplishes a better prediction of controllers as few NPIs as possible. With 70% of the dataset of NPIs as the training data, we use the mean absolute error to test the predictive accuracy for the rest 30% of the dataset. Fig. 3 shows that including three NPIs in the model gains the highest predictive accuracy for the designed controllers (see supplementary text for the details). We find that the most representative NPIs for Eq. 3 are: (1) stay-at-home order, represented by the normalized ratio of excessive time of staying at home, *θ*_{s}; (2) face-mask wearing, represented by the fraction of people wearing face masks, *θ*_{f}; (3) testing, represented by the normalized fraction of tested population, *θ*_{g}. Commonly, stay-at-home order and face-mask wearing have a positive impact on decreasing the number of reported cases. Differently, testing *θ*_{g}(*t*) may positively or negatively impact the reported infected cases. The reason is that more testing could allow identifying more cases when the number of testing is not sufficient, and yet, more testing may also lead to less infected cases^{6}. Representing the three NPIs as *θ*_{I} = [*θ*_{s}, *θ*_{f}, *θ*_{g}], then,

In the following studies, we only use these three selected NPIs for predictions.

### Effects of NPIs on shaping the disease dynamics

Based of the parsimonious model of Eq. (4), we learn the parameter-by-intervention specific marginal effects , and , reflecting the variations of *Û*_{Θ} as a function of the evolution of the time-dependent local NPIs *θ*_{s}(*t*), *θ*_{f} (*t*), and *θ*_{t}(*t*) at MSAs. With the objective to directly infer the *Û*_{Θ} with the NPIs data, we first use 70% dataset of NPIs as the training data to learn the marginal effects. Then, we estimate the counterparts *Û*_{Θ} given as Eq. (4) using the 30% left testing datasets and evaluate the fit between *U*_{Θ} and *Û*_{Θ}. Taking the “New York” MSA and “Houston” MSA shown in Fig. 4a-b as illustrative examples; given the evolution of NPIs, the estimated control signals *Û*_{β} and *Û*_{δ} with the learned marginal effects, fit the predicted model-based control law defined in terms of *U*_{β} and *U*_{δ}. Based on the proportionality between *U*_{γ}(*t*) and *U*_{β}(*t*) (see Fig. S2), i.e., *U*_{β}(*t*)*/U*_{γ}(*t*) = 2.664, here, we only consider *U*_{β}(*t*) and *U*_{δ}(*t*).

It is remarkable that for all MSAs, Eq. (4) stays robust to the huge heterogeneity of locality. As depicted in Fig. 4c, the NPIs have a great high standard deviation (*σ*) across MSAs, especially for stay-at-home order and face-mask wearing. One can notice that the average magnitude of stay-at-home order decreases from 80% to 40% recently with *σ* = 0.140. Besides, nearly 53.56% of people are wearing face masks when being outdoors with *σ* = 0.20 while the average magnitude of testing increases to 9.19% with *σ* = 0.032. To test the robustness of the model of Eq. (4), we trained and validated the model iteratively on different MSAs’ datasets for NPIs. The distributions of marginal effects across the three NPIs are shown in Fig. 4d. Most MSAs’ marginal effects for stay-at-home order and face-mask wearing are negative. The statistics imply the generic feature of the Eq. (4) in capturing NPIs’ effects despite the huge regional heterogeneity of human behaviours. However, the average marginal effect of testing has two opposite outcomes. For *Û*_{β} (*t*), 73.5% of MSAs’ testing’ marginal effects are positive, meaning 26.5% of them are negative. For *Û*_{δ} (*t*), 46.0% of MSAs’ testing’ marginal effects are positive, meaning 54.0% of them are negative. The positive effects suggest that testing is favorable to finding more infections and deaths, while the negative effects show that testing reduces infection or death.

Applying the estimated control signals *Û*_{β} (*t*), *Û*_{δ} (*t*), and *Û*_{γ} (*t*) = *Û*_{β} (*t*)*/*2.264 to the SIRD model, we find that the new output trajectories fit the reported infection *X*_{d} with *R*^{2} ≥ 0.9, as show in Fig. S3. As well, Fig. 4e-f depict the predicted infected/dead cases with reported infected/dead cases on 20 October 2020 and 20 February 2021. All the results further validate the model in assessing NPIs’ effects on disease dynamics.

### Needed NPIs and vaccinations for achieving *R*_{e} to manageable level

The two-steps approach, integrating the designed controllers in Eqs. (1-2) and the effect of NPIs on controllers in Eqs. 4, has successfully mapped human behaviors to NPIs to infection rate, recovery rate, and death rate of the SIRD model. Then, the effective reproductive number *R*_{e} can be equivalently translated in terms of magnitudes of NPIs (*θ*_{I}) and *V* (the ratio of people full vaccinated with efficacy rate 90% ^{34}), that is,
From the equation above, we can determine the needed magnitude of NPIs for achieving *R*_{e}(*t*) *<* 1, criteria for the disease die out, under a given vaccination coverage. Here, the vaccination coverage is viewed as an open-loop control action or a given parameter whose assigned value affects the speed of propagation of the disease.

When there is no vaccination administered in the US (like before 12 January 2021 *V* = 0), the needed magnitude of NPIs to achieve an effective reproductive number below the threshold, *R*_{e}(*t*) *<* 1^{37}, for ‘New York’ and ‘Houston’ MSAs are shown in Fig. 5a-b. The “triangular prism” in three dimensions represents the needed magnitude of the stay-at-home order, face-mask wearing, and testing. The horizontal slices of the “triangular prism” are the visualization for the needed magnitude of stay-at-home order and face-mask wearing if the magnitude of testing is fixed. If stay-at-home order and face-mask wearing interventions are enforced at a level greater than about 80%, *R*_{e}(*t*) *<* 1, which forms a “triangle” as shown in Fig. 5a-b. As opposed to stay-at-home order and face-mask wearing, the testing intervention has a positive marginal effect on ‘New York’ MSA and ‘Houston’ MSA. Hence, in both cases, the “triangle” becomes smaller for larger testing. Some MSAs (‘Log Angles’ MSA and ‘Miami’ MSA in Fig. S4) have a negative marginal effect for testing intervention, and their “triangle” becomes bigger for larger testing. According to their testing interventions ‘ marginal effects, the shapes of “triangular prism” for more MSAs are shown in Fig. S4.

Since 12 January 2021, the US began the first jab of vaccine, and undoubtedly, *R*_{e}(*t*) could be smaller with mass vaccinations according to the Eq. (5). For different ratios of fully vaccinated people with 90% efficiency, we could find at least how much NPIs could be relaxed (reduced) to achieve *R*_{e}(*t*) = 1 in Fig. 5c. Setting aside testing intervention, more magnitude of stay-at-home order and face-mask wearing interventions could be eased if more people are fully vaccinated. However, when 60% of people are fully vaccinated, only 55.3% of stay-at-home order and 64.8% face-mask wearing could be eased. When 70% people are fully vaccinated, the stay-at-home order and face-mask wearing could be eased entirely. However, to achieve *R*_{e}(*t*) *<* 1, the policymakers should relax NPIs with more cautions.

### Needed days of eliminating the pandemic locally with vaccinations

Moreover, due to the unknowns around mutations in escalating the disease’s infectiousness^{38} and human’s long-term immunity to the disease^{39}, suppressing the disease to an acceptable level with NPIs but not aiming to eliminate it could make COVID-19 an endemic^{40, 41}. A COVID endemic could claim millions of more lives each year and cause devastating economic burdens on immunization, treatment, and prevention. Also, repeated tightening and loosening interventions due to recurrent outbreaks, without a doubt, increase the difficulty to decide the right moment of enforcing exit strategies^{14, 35, 42, 43}. Currently, three vaccines are authorized and recommended in the United States from Pfizer-BioNTech, Moderna, and Johnson & Johnson. By 10 March, more than 37.4 million adults in the U.S. have been fully vaccinated^{44}. Despite continued efforts of NPIs and vaccine roll-out, the plateaued numbers of infection cases and the fast-spreading variants still pose significant uncertainty about whether “zero COVID” could be achieved locally in the US. Given the examples of New Zealand^{45}, Vietnam, Brunei, and Island states in the Carbbean^{40}, which ever reached a “Zero COVID” stage, we investigate the possibility of reaching zero contamination in the MSAs. We refer to “Zero COVID” to describe a situation that leads to no new cases at least for three months in a given location, referring to literature^{46}. Furthermore, the day when the following three months having zero new cases is defined as the elimination day.

We assume that the ratio of fully vaccinated people (also call vaccination coverage), follows the innovation adoption model according to Rogers^{47}, which is a normal distribution, *H*(*t*), see Eq. (16)-(17). Then, the total number of immunized people is defined as *V* (*t*) = *aH*(*t*) where *a* is the vaccine effectiveness (*a* = 90%). As shown in Fig. 6a, the reported data reveal that the ratio of fully vaccinated people from 12 January 2021 to 20 February 2021 in the US increases to 5.40%. By fitting the real-world ratio, we obtained the vaccination coverage as a curve that gradually reaches saturation within different periods (*T*_{s}). Considering a full vaccination of 50% of people after 360 days since 12 January 2021, the average elimination day for all MSAs corresponds to the hundred seventieth day. Under this configuration, from in Fig. 6b, one can notice that ‘New York’ MSA reaches a total elimination stage after 197 days while ‘Houston’ MSA achieves a total elimination stage after 187 days. Of course, one could define a vaccination strategy rollout accounting for nonuniform saturation coverage in different periods. Fig. 6c comprehensively presents all MSAs’ elimination days for vaccination coverage ranging from 0% to 90% with vaccination period (*T*_{s}) equals to 90 days, 180 days, 270 days, and 360 days. Our results reveal that broader vaccination coverage further reduces both the remoteness of the elimination day and the death toll. On the other side, a larger vaccination period expands the time range needed to reach the elimination day and leads to a more significant death toll. For the case in Fig. 6b, a total of 23,920,837 more people would be infected, and 1,199,202 more people would be dead because of COVID-19 if 50% of people are fully vaccinated in 360 days. On the contrary, if 90% of the population is fully vaccinated in 360 days, the average elimination day for all MSAs is about 129 days after implementing the vaccination strategy, with 20,428,296 more infected cases and 9,127,85 more dead cases overall. As shown in Fig. S6-S7, “New York” MSA, ‘Los angels’ MSA, “Chicago” MSA, “Dallas” MSA, “Phoenix” MSA, “Boston” MSA and “Philadelphia” MSA have far more death than other MSAs, meaning they need stricter interventions in reducing death because of COVID-19. In summary, all these results emphasize that the earlier, aggressive vaccination strategies are deployed in each MSA, the earlier the pandemic could be eliminated, leaving fewer deaths from the disease.

## Discussion

We present the COVID-19 pandemic as the feedback control system to show how the evolution of disease adapt to human behaviours to NPIs and the vaccine coverage. This approach enable us to model the linear and nonlinear effects of multiple NPIs on the SIRD model-based feedback controllers on the epidemiological parameters. Through reducing eight NPIs to three representative three (i.e., stay-at-home order, face-mask wearing, and testing), we obtain the model having great explanatory predictive power toward disease dynamics. By studying the NPIs at 381 MSAs from 1 April 2020 to 20 February 2021, we find that the two-steps approach is robust and efficient to predict daily infection and death accurately. Beyond in line with existing studies that NPIs are effective in suppressing the disease outbreaks^{30, 48, 49}, we could directly link the NPIs’ marginal effects to the effective reproductive number *R*_{e}. It implies that, without the need for up-to-date knowledge of current infections and ‘nowcasting’. Besides accurate forecast on both case counts and deaths, this approach could provide practical information for policymakers regarding the extent of safely relaxing NPIs and the needed vaccination coverage to end COVID-19 locally. The approach is also universal and thus can be used by policymakers elsewhere.

Through analyzing the needed NPIs for keeping *R*_{e} *<* 1, we find that MSAs should continue to enforce their stay-at-home order and face-mask wearing. Loosening the degree below 80% could lead to a resurgence of COVID-19. Our results show that MSAs should continue to require wearing masks and unless 60% of people are fully vaccinated, the MSA should not relax the stay-at-home order and face-mask wearing intervention.

There is an ongoing debate on if the elimination of COVID-19 (”Zero COVID”) is possible^{40, 41}. It is considered by some as the only way to prevent future crises and cannot be achieved without mass vaccination. Using the “Zero COVID” as a guiding criterion, i.e., no new infection at least for three months, we test how many days it can be achieved if 50% MSA populations are vaccinated in a 360-day period. We find that, on average, it takes 170 days to reach the “Zero COVID” milestone. It is worth noting that the promising result can be undermined by many external factors, like importation cases from other areas, making the total elimination a long process with fluctuations on new infections. In addition, some scientists argue pursuing the “Zero COVID” comes at a huge cost to normal human life^{40, 41} and can thus out-weigh the benefits. Nevertheless, our experiments reinforce the argument that maintaining NPIs and encouraging people to get vaccines are necessary and key strategies to lower the new infections as much as possible.

Our study has some limitations. The first limitation is rooted in the dataset for COVID-19 and the dataset for NPIs. Given the mass mild or asymptomatic infections, the inaccuracy of reported infections and deaths would increase the uncertainty of the SIRD model-based feedback controllers. The second limitation is assuming each MSA as the closed population, which ignores the fluctuation of infections caused by case importation and exportation. This simplification would influence the results from needed NPIs and PIs for suppressing COVID-19 to elimination test. The third limitation is assuming vaccination adoption follows the curve of diffusion of innovations without considering people’s attitude to vaccinations. As the attitudes towards vaccination vary by age, race, ethnicity, and education, it is hard to capture the full complexity. In our current work, we use feedback linearization to design the control signals, and we may improve the accuracy by considering other controller design strategies, for example, adaptive controller^{50, 51}, model predictive control^{52}, or intelligent control^{53}. Nevertheless, our study provides practical insights into tightening or relaxing NPIs for the aim of living with COVID-19. Also, we provide the possibility of achieving “Zero COVID” in metropolitan areas if vaccination is stable and efficient enough.

## Methods

### Dataset

#### 2019-COVID pandemic

We acquire the second-administrative units’ COVID-19 cases from the 2019 Novel Coronavirus Visual Dashboard operated by the Johns Hopkins University Center for Systems Science and Engineering^{54}. It covers the counties’ reported infected and dead cases at the US and the whole countries’ infected and dead cases from 21 January 2020 to 20 February 2021. We integrate counties’ COVID-19 cases for 381 MSAs defined by the United States Office of Management and Budget (OMB), which serves as a high degree of social and economic core areas.

#### Testing Capacity

The testing data at each county is also collected by the Johns Hopkins University Center for Systems Science and Engineering^{54}. By integrating the counties’ testing data, we get each MSA’s daily ratio of total testing over MSA’s population.

#### Stay-at-home Data

Daily data about the observed minutes at home and observed minutes outside of home for all devices are counted at each county by the Safegraph^{55}, which is a platform of collecting the points of interests (POIs) from anonymous mobile devices. Integrating the data, we could get each MSA’s ratio of time at home, i.e, the fractions of the observed minutes at home over the sum of observed minutes at home and outside home, which is reflective of people mobility pattern for local governments’ anti-contagious policies. By subtracting the benchmark ratio (the average ratio of time at home before 13 March 2020) and then normalizing the ratios from 0 to 1, we get each MSA’s ratio of excessive time at home.

#### People Wearing Face Mask

We assume that the search of “face mask” collected by Google Trend equals to the fraction of people having the awareness of wearing face masks in each MSA. Combined the fractions and the weekly percentages of people wearing face mask in the US provided by YouGov^{56}, we could obtain each MSA’s daily percentage of people wearing face mask. We use the similar approach to collect the data of people support school closure, quarantine, working from home, frequent hand wash, and avoid crowding.

#### Vaccination Data

CDC provides the overall US COVID-19 Vaccine Distribution and Administration, mainly for Pfizer-BioNTech and Moderna.

### Deriving the feedback controllers

The feedback controllers measure the output of the SIRD model and then manipulates the inputs on infection rate, recovery rate, and death rate of as needed to drive the model output toward the desired COVID-19 pandemic trajectory. In another word, given the controllablility of the SIRD model (see supplementary text), then we could rewrite the Eq. (1) as
with the controller set *U*_{Θ}(*t*) = [*U*_{β}, *U*_{γ}, *U*_{δ}]^{T}. *U*_{β}, *U*_{γ}, *U*_{δ} are the controllers on infection rate, recovery rate, and death rate. For the reference COVID-19 pandemic trajectory *X*_{d} = [*S*_{d}, *I*_{d}, *R*_{d}, *D*_{d}]^{T} and its corresponding output trajectory *X*, the error dynamics *X*_{d} −*X* is governed by the following equations
where, , and . Here, *k*_{1}, *k*_{2}, and *k*_{3} are positive gains of the designed feedback controller. The solution to the linear ordinary differential equations above are expressed as
and . Therefore, the error system is exponentially stable and *X*(*t*) → *X*_{d}(*t*). The simulation of the closed-loop system is performed selecting *k*_{1} = *k*_{2} = *k*_{3} =0.1. The controllers are
To prevent un-physical negative transmission rate *U*_{β}(*t*), recovery rate *U*_{γ}(*t*), and death rate *U*_{δ}(*t*), the controllers are clipped such that the time-dependent epidemiological parameters are always positive. Here Θ_{0} = [*β*_{0}, *γ*_{0}, *δ*_{0}]^{T} is the pre-intervention epidemiological parameters learned by infection of ‘New York’ MSA from 1 March 2020 to 13 March 2020 with Nelder-Mead simplex algorithm^{57}. For the controllers, whatever the values for Θ_{0}, the error system is stable and *X*(*t*) → *X*_{d}(*t*).

### Using difference-in-difference method to learn the NPIs’ marginal effects to the controllers

For NPI set *θ*_{I}, we use the linear regression model for difference-in-difference to calculate their effects on the vector of infection rate, recovery rate, and death rate *U*_{Θ}(*t*) + Θ_{0}. Consider the model,
where *λ*(*t*) is the facotr for time trend and *ϵ*_{Θ}(*t*) is the residual term. Then, the difference of outcome controllers from time *t* − 1 to time *t* is
Adding up the all difference from time *t* = 0 to time *t*
As when *t* = 0 no NPIs are implemented, thus, *U*_{Θ}(*t*) = 0 and *f* (*θ*_{I}(0)) = 0,

### Vaccination adoption model

Like the innovation adoption model, the daily newly vaccination adoption (vaccination coverage) will follow the bell curve, the normal distribution,
where is the saturation of vaccination coverage. The the cumulative vaccination adoption (vaccination coverage) is
with *t* ∈ [1, *T*_{s}]. It means, the cumulative vaccination coverage will reach saturation of in *T*_{s} days. In the test of “Zero COVID”, we tune the *u* and *s* to fit the real-world vaccination coverage in the United States from 12 January 2021 to 20 Februray 2021 for ∈ [0.1, 0.9] and *T*_{s} ∈ {120, 180, 250, 360} in Fig. 6.

## Data Availability

Derived Data and codes that support the findings of this article are available in https://github.com/lucinezhong/controllers_ode_COVID.git.

## Funding

This work was supported in part by the Department of Mechanical Aerospace and Nuclear Engineering Department at Rensselaer Polytechnic Institute, Troy, NY.

## Author Contributions

M.D. and J.G. conceived the project and designed the experiments; Q.W. and L.Z collected the data set and analyzed the data; L.Z., M.D., and J.G. carried out theoretical calculations and performed the experiments; all authors wrote the manuscript.

## Competing Interests

The authors declare that they have no competing financial interests.

## Data and materials availability

All data needed to evaluate the paper’s conclusions are presented in this paper. Derived Data and codes that support the findings of this article are available in https://github.com/lucinezhong/controllers ode COVID.git.

## Additional information

Supplementary text including Fig. S1-S7, Tab. S1, and Eqs. (S1)-(S12).

## Acknowledgements

We thank Bhathiya Rathnayake for fruitful discussions on developing the feedback system and the paper.