Vaccination and three non-pharmaceutical interventions determine the dynamics of COVID-19 in the US

The rapid rollout of the COVID-19 vaccine 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 model that maps human social behaviors into 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 accurately in 381 MSAs and avoid over-fitting. Our numerical results demonstrate our approach’s excellent predictive power with R2 > 0.9 for 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 Re to a manageable level and how the variants of concern diminish the likelihood for disease elimination at each location. Our analytical results provide insights into the debates surrounding the elimination of COVID-19. NPIs, if tailored to the MSAs, can drive the pandemic to an easily containable level and suppress future recurrences of epidemic cycles.


Introduction
T he ongoing global pandemic of coronavirus disease 2019  has caused devastating loss of human lives and inflicted severe economic burden in the US (WHO, 2020). By 20 March 2021, more than 510,000 people were killed by COVID-19, the unemployment reaches 11.5% (CBO, 2021), and fiscal shortfall reaches over 200 billion (The Council of State Governments, 2020). It is known that findings from scientific research on non-pharmaceutical interventions (NPIs) and pharmaceutical interventions (e.g., available vaccinations (Le et al., 2020;Bubar et al., 2021;Markoviv et al., 2021) and drugs (Rabby, 2020)) have successfully guided policymakers on implementing public health strategies during a seasonal or pandemic flu. To control the COVID-19 pandemic, countries of the world have deployed a wide range of (NPIs) (Haushofer and Metcalf, 2020;Priesemann et al., 2021), including personal NPIs (e.g., home isolation, hand hygiene, face mask), community NPIs (e.g., school closure, social distancing, limiting for mass gatherings), and environmental hygiene (e.g., surface cleaning) (CDC, 2021). Methodologies, from randomized controlled trials (Haushofer and Metcalf, 2020), econometric methods (Hsiang et al., 2020;Cho, 2020) to mathematical models (Dehning et al., 2020;Chen et al., 2020;Ruktanonchai et al., 2020;Brauner et al., 2021;Zhong et al., 2021;Yang et al., 2021;Wang et al., 2022), have been developed to measure the effects of these NPIs. Table 1 provides a brief overview of some existing methodologies for estimating NPIs' effectiveness at reducing COVID-19 transmission. But many of them are simulation studies, which often cannot be used without predicting the pandemic when the NPIs are adjusted. The other data-driven studies, which mostly incorporate binary data of whether the NPIs are implemented instead of human behaviors/opinions toward interventions, would fail to capture the impact of human and social development. Moreover, 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 (Yan et al., 2015), 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, manufacturing, aircraft, satellite, and robots, has 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 (Liu et al., 2011) or the lowest control costs (Yan et al., 2015(Yan et al., , 2012 for node control, edge control (Nepusz and Vicsek, 2012), target control (Gao et al., 2014), multilayer control (Pósfai et al., 2016;Menichetti et al., 2016), temporal control (Li et al., 2017), and data-driven control (Baggio et al., 2021). 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, for example. to control a vehicle's speed, 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 (Liu and Barabási, 2016;L'´ober, 2016). 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 (L'´ober, 2016). We validate this hypothesis by developing a parsimonious model that predicts how the interventions influence the spreads in 381 MSAs in the US and ultimately estimates the control of COVID-19.

Model
Model COVID-19 spreading as a feedback control system. Increasing evidence shows that the spread of COVID-19 follows compartmental models (Hsiang et al., 2020;Lai et al., 2020;López and Rodó, 2020;COVID et al., 2020), such as the SIRD (Susceptible-Infectious-Recovered-Death) model, which is mathematically described by the nonlinear equations expressing a population balance as follows (Keeling and Rohani, 2011): 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. Here, V is the count of people fully vaccinated with an initial efficacy rate 90% (Thompson et al., 2021) (see that V(t) changes as a function of the cross variant immunity and waning immunity in Eq. (14)). The parameter μ is the crude birth and death rate, and the epidemiological parameters β, γ, and δ are the infection, recovery, and death rates (Θ = [β, γ, δ] T ). Studies show that epidemiological parameters are time-dependent, which adapt accordingly to the change in interventions (Dehning et al., 2020;Chen et al., 2020). By defining β 0 , γ 0 , and δ 0 as the infection, recovery and death rates before interventions (Θ 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 (Nepusz and Vicsek, 2012) and the vaccination V works as the node control (Liu et al., 2011). Based on nonlinear feedback control law, we develop the controllers U Θ (t) using the feedback linearization approach. Subject to our designed control actions, 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 three-dimensional data when the feedback effect is included, while it fails to fit when it is 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. Specifically, the human behavior toward NPIs θ I ¼ ½θ s ; θ f ; θ g T , which are designated as stay-at-home order θ s , face-mask wearing θ f , and testing θ g , are measured according to Google Trend and survey data from YouGov (see Dataset for the details). The designed controller, U Θ ¼ ½U β ; U γ ; U δ T , are the changes in preintervention infection, recovery and death rates. Then, this model could measure the effect of NPIs θ I on the controllers U Θ by comparing the changes in infection dynamics (designed controller) before and after the same region's NPIs deployment. Beyond existing models, which assume the effects on policies are approximately linear (Hsiang et al., 2020), we also identify the interactions between NPI policies. Thus, we are able to compile the compound evolution of human behavior toward the NPIs, to the evolution of the control signals with their respective effects, In short, the feedback control system characterizes how the input of human behavior toward NPIs determines the output of designed controllers (changes in pre-intervention infection, recovery, and death rates) in the first step. And it exhibits how the designed controllers and vaccine rollout determine the final output of disease dynamics in the second step. Thus, the twosteps approach captures how the NPIs and vaccine rollout govern the disease dynamics, through combining the principles of control theory (L'´ober, 2016;Stewart et al., 2020) and statistical estimation (Hsiang et al., 2020;Gertler et al., 2016).
Deriving the feedback controllers. The feedback controllers measure the output of the SIRD model and then manipulate the inputs on infection rate, recovery rate, and death rate as needed to drive the model output toward the desired COVID-19 pandemic trajectory. 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: are the differences/errors between real-word data and simulated data. 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 IðtÞ ¼ e Àk 1 tĨ ð0Þ;R ðtÞ ¼ e Àk 2 tR ð0ÞDðtÞ ¼ e Àk 3 tD ð0Þ; ð4Þ and the errors ðĨðtÞ;R ðtÞ;DðtÞÞ ! ð0; 0; 0Þ. 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 The control is subject to a positivity constraint. To prevent nonphysical transmission rate U β (t) + β 0 , recovery rate U γ (t) + γ 0 , and death rate U δ (t) + δ 0 , the controllers are designed as The SIRD-model based feedback control system reveals how the interventions' magnitudes and timing govern the dynamics of disease. a Based on nonlinear feedback control law, we develop the controllers U Θ ¼ ½U β ; U γ ; U δ T working as the edge control on infection rate, recovery rate, and death rate. On the other side, the vaccination V works as the node control on the susceptible population. b With these controllers, the output trajectory of disease X = [S, I, R, D] T of the feedback control system fit with the real-world infection data of the disease. c Through the model, which includes the nonlinearity or interactions between NPIs, the human behavior toward NPIs (e.g., stay-at-home order θ s , face-mask wearing θ f , and testing θ g ) are linked to the controllers U Θ through their effects (e.g., w s Θ , w f Θ , w g Θ ). Icon credits for c: stay-at-home logo, Tinypolkadoz/Dreamstime.com; face mask, Elenabsl/Dreamstime.com; testing, Zubada/Dreamstime.com; vaccination, www.creazilla.com. 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 (Gao and Han, 2012). For the controllers, irrespective to the values of initial condition Θ 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 the 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 factor 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 0 to time t ð10Þ 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, whereV is the saturation of vaccination coverage. The cumulative vaccination adoption (vaccination coverage) is It means, the cumulative vaccination coverage will reach saturation ofV in T s days. For testing the case of a well-controlled COVID 19, we tune the u and σ to fit the real-world vaccination coverage in the United States from 12 January 2021 to 12 January 2022 in Fig. 6.
To consider the waning nature of vaccine (Mallapaty et al., 2021), we define aðt À t 0 Þ as the decreasing vaccine efficacy since t 0 with initial efficacy 0.9, where b is the decreasing rate. When 0:9 < bðt À t 0 Þ, the vaccine totally loses its efficacy and the vaccinated people again become susceptible. Then, the effectively vaccinated population at t is where hðt 0 Þ is the count of new vaccinated population at t 0 .

Results
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, are 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 (2)- (7)) in Model). The feedback law ϕ(.) relies on the measurements of the model full state X and the reference trajectory (X d ) and their time derivative ( _ X and _ X d ). 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 (Pei et al., 2020), "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 R e ðtÞ ¼ ½β 0 þU β ðtÞSðtÞ ½γ 0 þU γ ðtÞþ½δ 0 þU δ ðtÞþμ , 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 threshold 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 control 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 behaviors 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 HUMANITIES AND SOCIAL SCIENCES COMMUNICATIONS | https://doi.org/10.1057/s41599-022-01142-3 ARTICLE HUMANITIES AND SOCIAL SCIENCES COMMUNICATIONS | (2022) 9:149 | https://doi.org/10.1057/s41599-022-01142-3 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 loosened 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Û Θ ðtÞ ¼ ½U β ; U γ ; U δ T ¼ f ðθ I ðtÞ; W I Θ Þ as (see [Eqs. (8)-(11)] in Model) whereÛ Θ ðtÞ is the estimate of the control action based on NPIs with their magnitudes θ j (t) and their impact value w j Θ . θ I is the vector representing different NPIs θ I ¼ ½θ 1 ; θ 2 ; :::; θ i ; :::; θ n T . W Θ is a vector representing the effects of NPIs at θ I and W I Θ ¼ ½w 1 Θ ; w 2 Θ ; :::; w i Θ ; :: For each specific NPI j, large w j Θ value demonstrates that NPI j has strong impact. If w j Θ ¼ ½0; 0; 0 T , 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 effect on the controller. For either community NPIs ϑ c or personal NPIs ϑ p , 1 À ∑ θ j 2ϑ i w j Θ θ j ðtÞ 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 (see Table 2), we evaluate the goodness of fit versus different combinations of NPIs to find an effective parsimonious model for Eq. (16). A parsimonious model, which has great explanatory predictive power, accomplishes a better prediction of controllers with 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. Figure 3 shows that Fig. 2 The designed feedback controllers U Θ drive the output trajectories ( _ X) in the SIRD model ingeniously fits the reported trajectories ( _ X d ) for 381 MSAs. a The reported data (dots) of infected/recovered/dead cases are fitted with the output trajectory (solid line) for example MSAs, i.e., "New York" MSA and "Houston" MSA. b Comparison between the reported data and output data at all MSAs. Each MSA represent a dot. c All MSAs' temporal infection rate, [β 0 + U β (t)], recovery rate, [γ 0 + U γ (t)], and death rate, [δ 0 + U δ (t)] with marking out the examples MSAs' rates. d shows the MSA's temporal effective reproductive ratio (R e ) and e shows each MSA's average effective reproductive ratio from 14 February 2020 to 20 February 2021 in the cartogram map, in which the geometry of regions is distorted according to their population. The effective reproductive ratio reaches the lowest as of 20 February 2021.
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. (16) 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 facemask 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 fewer infected cases (Chiu et al., 2020). Representing the three NPIs as θ I ¼ ½θ s ; θ f ; θ g T , 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. (17), we learn the parameter-byintervention specific marginal effects w s Θ , w f Θ , and w t Θ , reflecting the variations ofÛ Θ as a function of the evolution of the timedependent 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. (17) 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. S1), 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. (17) 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. (17), we trained and validated the model iteratively on different MSAs'  , stay-at-home order, school closure, quarantine, working from home, face-mask wearing, testing, frequent hand wash, and avoid crowding, we test the accuracy of the predictive model with different combinations of NPIs for U β (b), U γ (c), and U δ (d). The model achieves high parsimony (with fewer NPIs) and a high level of goodness of fit (with the lowest mean absolute error) using three NPIs, which are stay-at-home order, face-mask wearing, and testing.
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 Eq. (17) in capturing NPIs' effects despite the huge regional heterogeneity of human behaviors. 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 U γ ðtÞ ¼Û β ðtÞ=2:264 to the SIRD model, we find that the new output trajectoriesX fit the reported infection X d with R 2 ≥ 0.9, as shown in Fig. S2. 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 a manageable level. The two-steps approach, integrating the designed controllers in Eqs. (1)-(15)) and the effect of NPIs on controllers in Eqs. (17), 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 an efficacy rate 90% (Thompson et al., 2021)), that is, From the equation above, we can determine the needed magnitude of NPIs for achieving R e (t) < 1, the criterion for the disease die out, under a given vaccination coverage. Here, the vaccination Fig. 4 The learned effects enable us to use the magnitude of NPIs to predict infection at each MSA. Using the 70% NPIs and infection data, we learned the NPIs' effects in the model of Eq. (17). Based on the learned effects, a, b illustrates how the magnitude of NPIs determine the controllers (Û β ðtÞ and U δ ðtÞ) in the example MSAs. To test the accuracy of model, we compare [Û β ðtÞ andÛ δ ðtÞ] with analytical controllers [U β (t) and U δ (t)] with the left 30% NPIs and infection data (vertical shaded area) in the right-side plots of a, b. c All MSAs' magnitudes of NPIs. The solid line represents the average, and the shaded area represents the standard variance for all MSAs. Though the large difference between MSAs' magnitudes of interventions, in d, the marginal effects for stay-at-home order and face-mask wearing are mainly negative. Testing, which may help to find more infection and death, have either negative or positive effect. Thus, by taking the estimated controllersÛ Θ ¼ ½Û β ;Û γ ;Û δ as the input of the SIRD model, we find the estimated infection/death assemble with reported infection/death with R 2 > 0.9, see the two example dates in e and f. Figure S2 shows the estimation accuracy for all others dates. It should be noted that as theÛ γ keeps constantly proportional toÛ β (see Fig. S1), there is no need for further exploringÛ γ .
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 is below the threshold, R e (t) < 1 (Dowdle, 1998), 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, facemask wearing, and testing. The horizontal slices of the "triangular prism" are the visualization for the needed magnitude of stay-athome 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-athome order and face-mask wearing, the testing intervention have 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. (18). 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% of 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 caution.
When the variants are considered in the model. The world has endured multiple variants of the SARS-CoV-2, and there is a possibility that COVID-19 will become endemic. Although it means large-scale NPIs enforced by governments are likely to lessen, it is worth noting that it also means the world may bear long-term suffering. A COVID endemic could claim millions more lives each year and cause devastating economic burdens on immunization, treatment, and prevention (Lee et al., 2020;Heywood and Macintyre, 2020). Additionally, it is possible that repeated tightening and loosening interventions are needed in the future due to recurrent outbreaks (Stewart et al., 2020;Metcalf et al., 2020;Thompson et al., 2020). Therefore, it is valuable to explore the "zero COVID," which was achieved temporarily in New Zealand (Baker et al., 2020), Vietnam, Brunei, and Island states in the Carbbean (Lee et al., 2020). We utilize the vaccination data (CDC, 2021) and investigate the possibility of theoretically reaching zero infections, i.e., no new cases at least for three months in a given MSA (Diekmann et al., 2012) when within/ without variants. Furthermore, the day when the following 3 months having zero new cases is defined as the curbing day.
We assume that the ratio of fully vaccinated people (also called vaccination coverage), follows the innovation adoption model according to Rogers (Oldenburg and Glanz, 2008), which is a normal distribution, see Eq. (12)). Then, the total number of immunized people is defined as V(t) in Eq. (14) with the waning vaccine effectiveness. As shown in Fig. 6a, the reported data reveal that the ratio of fully vaccinated people from 12 January 2021 to 18 February 2022 in the US increases to 64%. We obtained the vaccination coverage as a curve that gradually reaches saturation show the needed magnitude for NPIs in order to keep R e ≤ 1 with zero vaccine coverage. 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. c When people are fully vaccinated at different levels (vaccine coverage), how much extents of stay-at-home order and face-mask wearing could be relaxed while keeping R e = 1. The magnitude of testing is fixed, the same as each MSA's recent testing capacity. These results demonstrate relaxing interventions need thorough and careful consideration when the ratios of full vaccinated people <0.6. within example periods of 270 and 360 days (T s ) by fitting the real-world ratio. Considering a full vaccination of 60% of people after 360 days since 12 January 2021, the average curbing 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 the curbing stage after 197 days while 'Houston' MSA achieves a total curbing stage after 187 days. These results reveal that broader vaccination coverage further reduces both the time until the curbing day and the death toll. However, when we consider variants of concern: Alpha, Beta, Gamma, Delta, and Omicron in the model, we find that curbing day is hard to achieve in three years as new cases will re-emerge even after a temporary period of "zero new cases." We simply assume that infection rates increase and death rates decrease according to the emergence of variants (see Table S2). Though highly contagious variants would potentially exacerbate the infections, their less-deadly characteristics would help reduce the deaths. (see Fig. 6c, d).

Discussion
We present the COVID-19 pandemic as the feedback control system to show how the evolution of disease adapts to human behaviors to NPIs and the vaccine coverage. This approach enables us to model the linear and nonlinear effects of multiple NPIs on the SIRD model-based feedback controllers on the epidemiological parameters. By reducing eight NPIs to three representative ones (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 (Lai et al., 2020;Flaxman et al., 2020;Group et al., 2021), we could directly link the NPIs' marginal effects to the effective reproductive number R e . 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 manage COVID-19 locally. The approach is also universal and thus can be used by policymakers elsewhere.
By 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.  6 Disease severity with/without variants of concern. By (1) defining the day reaching zero new infection for 90 days as the curbing day, (2) assuming that interventions (i.e., stay-at-home order, face-mask wearing, and testing) keep the same as that from 5 January 2021 to 11 February 2021, (3) consider the ratio of people fully vaccinated as the vaccination coverage, we test the curbing days for MSAs since 12 January 2021. a The vaccination adoption model H(t) illustrates how the saturation of people fully vaccinated is reached in period T s . Specifically, the orange line is the reported cumulative vaccination coverage. b MSAs' estimated curbing days if 60% are full vaccinated with 360 days. c, d Box plots for MSAS' additional death in the following 1000 days for different saturation of vaccination coverage with/without the presence of variants.
It is currently being debated whether COVID-19 ("Zero COVID") can be eliminated and some view it as the only means of preventing future crises and insist that mass vaccination is necessary (Lee et al., 2020;Heywood and Macintyre, 2020). However, the goal is undermined by the new variants. The continuous emergence of new, highly contagious, and less deadly variants gives rise to a long process with fluctuations in new infections. The best outcome would be that future mutation of the COVID-19 virus keeps less lethal and becomes endemic. Thus, to persistently fight against COVID-19, our experiments reinforce the argument that maintaining NPIs and encouraging people to get vaccines are necessary and key strategies to lower the new deaths 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 the elimination test. The third limitation is assuming vaccination adoption follows the curve of diffusion of innovations without considering people's attitudes to vaccinations. As the attitudes towards vaccination vary by age, race, ethnicity, and education, it is hard to capture the full complexity. The fourth limitation is that our model outcomes are determined by the people's responses to NPIs, which are time-dependent. When it's being used for prediction, it would be more suitable for predicting short-term rather than long-term dynamics especially when the NPIs deployment changes. Nevertheless, our study provides practical insights into tightening or relaxing NPIs for the aim of living with COVID-19.

Data 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 at https://github.com/lucinezhong/ controllers_ode_COVID.git.