key: cord-0695438-j1huk63f
authors: Gulbudak, H.; Browne, C.; Macdonald, J. C.
title: Modeling COVID-19 outbreaks in United States with distinct testing, lockdown speed and fatigue rates
date: 2021-01-06
journal: nan
DOI: 10.1101/2021.01.04.21249231
sha: 6d9a556af9fe6ef32729874f20eed4cc489082df
doc_id: 695438
cord_uid: j1huk63f

Each state in the United States exhibited a unique response to the COVID-19 outbreak, along with variable levels of testing, leading to different actual case burdens in the country. In this study, via per-capita testing dependent ascertainment rates, along with case and death data, we fit a minimal epidemic model for each state. We estimate infection-level responsive lockdown entry and exit rates (representing government and behavioral reaction), along with the true number of cases as of May 31, 2020. Ultimately we provide error corrected estimates for commonly used metrics such as infection fatality ratio and overall case ascertainment for all 55 states and territories considered, along with the United States in aggregate, in order to correlate outbreak severity with first wave intervention attributes and suggest potential management strategies for future outbreaks. We observe a theoretically predicted inverse proportionality relation between outbreak size and lockdown rate, with scale dependent on the underlying reproduction number and simulations suggesting a critical population quarantine "half-life" of 30 days independent of other model parameters.

The COVID-19 pandemic began in Wuhan, China and spread rapidly throughout the world in early 2020 provoking a wide range of interventions. China gained rapid control of its epidemic, likely due to it's rapid and strict lockdown as discussed in several studies [4] , [7] , [15] , in contrast to several other countries with slower and less unified efforts. In particular, the United States employed a heterogeneous response, with different scales of quarantine measures (stay at home orders and business closures), despite the first reported case of COVID-19 occurring on January 21, 2020 in Washington state and cases being reported in all 50 states by mid-March [13] . Furthermore, each state displayed distinct exit strategies and lockdown fatigue even though cases never were brought down to low levels. Another feature of the US response was variable testing through time and by state, growing from sparsely to more widely available tests, leading to distinct increasing trends of case detection.

Each state represents distinct realizations of how different response characteristics and case ascertainment correlated to true case burden in the United States, which provides fertile ground for testing outbreak containment strategies. However, successfully fitting an epidemic model to this heterogeneous response exhibited in the first wave COVID-19 outbreak in the US presents a number of challenges. The multitude of epidemic curve shapes induced by the distinct dynamic behavioral and government interventions can lead to issues of over-fitting or a large number of parameters that blur the most important factors [8] . In addition, the unknown true number of cases calls for methods to incorporate testing, mortality or other sources of data (e.g. seroprevalence studies), which is complicated by changing quantities, such as detection [1] , [12] , [11] , or antibody levels [12] . Here we focus on a unified model both simple and flexible enough to fit the wide range of COVID-19 outcomes in the US (through May 31, 2020), which incorporates the critical factors in epidemic trajectory. Furthermore, we construct an appropriate per-capita testing dependent ascertainment rate calibrated with mortality data to both allow more accurate model fits and provide a means of estimating the actual number of cases.

The wide range of quarantine responses in the different states and resultant outcomes allow us to analyze potential responses to future outbreaks utilizing correlation and sensitivity analyses as well as counterfactual simulation studies. In this way, our minimal number of identifiable parameter values can be compared among the states representing crucial control quantities, such as lockdown speed and fatigue, and via analytically derived relations, we link outbreak size as inversely proportional to population quarantine rate. We also provide estimates of commonly used quantifiers of outbreak severity, epidemic trajectory and effectiveness of control measures, such as Infection Fatality Ratio (IFR) and the ratio of (estimated) true cases to reported cases, accounting for common sources of error in their prediction [14] , [6] . By estimating these quantities for all states and territories under one modeling framework we provide a comprehensive overview of the first wave outbreak in the United States. Ultimately we synthesize these results into a range of potential future responses which highlight the critical nature of widespread, swift quarantine measures, with a model suggested critical duration, for a rapidly spreading outbreak. These results can be used to inform management strategies both as COVID-19 vaccine is not yet widely available, and for future outbreaks of emergent pathogens.

. . . population, and infected (I) individuals whom ultimately progress to reported (R), unreported (U ), or dead (D) cases.

Our model assumes susceptible individuals become infected relative to force of infection λ(t) = βI(t) N , where β is the transmission rate and N is the population size. In addition susceptible quarantine acts proportionally to force of infection λ(t) with lockdown (rate) factor ψ, as a proxy for the population (or government mandated) reaction to incidence by self-quarantine. While individuals are in quarantine it is assumed that they do not contact with infected population, and that individuals exit quarantine with rate α, which we also call lockdown fatigue. Infected individuals either have their cases reported (model compartment R), or fail to have their cases reported (model compartment U), with mean infectious period T = 7.5 days, with the proportion of individuals in each compartment determined by ascertainment rate ρ(t), or the proportion of true cases at a given time captured by testing. Regardless of report status death is expected to occur with Infection Fatality Ratio (IFR), represented by model parameter ξ after a mean delay of µ = 21 days after infection. The model parameters and descriptions are also listed in table 1. In addition to the fixed parameter values for mean infectious period T and mean time until death µ, we fit the remaining parameters of the model as described in next section, obtaining good fits for the data from all 50 states, Washington D.C. Guam, Puerto Rico, the Northern Mariana Islands, and the US Virgin Islands (see Figure 7 and the SI figs. S1-S14). 

Lockdown (rate) Factor α Self-quarantine Exit rate (Lockdown fatigue) T mean infectious period fixed at T = 7.5

A/N +τ (t)/N + k 0 time variable ascertainment rate testing data, τ (t), from [3] k maximum ascertainment rate A half saturation point for testing k 0 minimum ascertainment rate k 0 = 0.03 for states, K 0 = 0.005 whole US fit µ mean time to death fixed at µ = 21 η gamma distribution shape parameter ξ Infection Fatality Ratio (IFR) All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint

We now describe the methods of our fitting procedure and analysis of the model. All relevant data and code are available at [9] . Where shown, indicated stay at home orders are based upon [2] . Population levels are from US Census Bureau 2019 projections [5] . To obtain parameter estimates for our model in each state and territory, we utilized cumulative death and case totals [10] as well as daily testing data, inferred from [3] . The data used for the whole US fit was obtained by combining data from all subsidiary states and territories due to concerns about data reliability at the federal level [16] . Even at the state level the raw testing data contained a number of reporting irregularities, as evidenced by several days with negative daily testing instance. To account for this we smoothed the data by averaging with nearest neighbors in such a way as to not change the cumulative testing totals prior to finding the moving average. We define the ascertainment rate function ρ(t) as

where τ (t) is the three day moving average of the raw daily testing totals, with parameters k as maximum ascertainment rate, A as half saturation point in terms of tests, and k 0 as minimum ascertainment rate. Observe that ρ(t) is chosen to be a saturating function of daily tests as a percent of population (see Figs 2 (b)-(d) and 4 (e). As the number of tests increased, so did the ascertainment rate because the increase in tests was not driven by more demand solely based on rising case numbers, but wider availability and convenience leading more individuals overall to take the tests. The negative concavity of ρ as a function number of tests reflects the principle of diminishing returns associated with wait times, and also a larger influence of demand due to actual rising cases once tests become widely available, which together lead to saturation of ρ to a maximum level. It is assumed even absent or at low testing, there is a minimal ascertainment rate calibrated based upon low ascertainment levels during the early stages of an outbreak with bounds of .005 and .03 for the whole of the United States and its territories, and for each state or territory respectively. This difference in lower bound was chosen because, as indicated by Figure 5a , many states' outbreak began significantly before their first reported case, thus when considering the United States as a whole, few cases were ascertained in the initial phase of the outbreak. For each state and territory fit and the whole United States fit, the maximum ascertainment rate was fixed based upon the maximum daily test count as a percent of the All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint . . . states population. This bound was set based in part on the estimated cumulative case total for Connecticut and New York City provided in [12] . The half saturation constants were fit to each state, calibrated by the simultaneous fit of mortality and reported cases, to reflect state specific relationships between testing and ascertainment (see Fig. 2 ). The obtained ratios of total to reported cases (see table 2 ) are in line with CDC estimates [1] .

We numerically solved the non-linear least squares problem, using the interior-reflexive Newton method as implemented by MatLab's lsqcurvefit function, minimizing objective function

where the ratio of cumulative deaths and reported cases at end-times, ω = D(t D f )/R(t R f ), was the chosen positive weight for each state. For each state or territory the time interval considered was from the first confirmed case up until 31 May 2020 for cases, by which time most state's mandated stay at home orders had ended [1] , and June 21, 2020 for deaths. For our fittings the mean infectious period was fixed at, T = 7.5 and the mean to death at µ = 21 to improve model identifiability. These values are in line with the weighted (by proportion of total cases) mean parameter values obtained from fitting the model without fixing these parameters (see program files at [4] ). All other model parameters, as well as I 0 , the initial case total, were fit for each state or territory. The fit values (see table 2 and the SI for all parameter values) were then used to provide approximate cumulative unreported case totals (model compartment U ) and true cumulative case totals (sum of model compartments U + R). We additionally checked if further model parameters could be fixed while retaining model robustness (see section 1 of the SI), but found that fitting both case ascertainment, ρ, and infection fatality ratio, ξ (bounded within a feasible range from .001 to .05), were necessary to obtain good fits across all states and territories.

In order to compare % model predicted positives (MPP) with % positives inferred from the data (PID) these quantities were plotted together with τ (t) for each state and territory,

where d(t) is the three day rolling average of case totals inferred from reported cumulative case totals after applying the same smoothing to daily case totals as was applied to daily testing incidence, and δ(t) is the daily reported case incidence inferred from the fit cumulative case totals (see Figures 4, 7) . For all states and territories, with the exception of the Northern Marianas Islands (due to that this territory's fit parameters indicating that all cases were imported), the relationship between model parameters and outbreak trajectory were assessed. The reproduction number R 0 can be utilized as a single parameter. Indeed, for model (2.1),

(noting that the susceptible quarantine does not affect R 0 because it is proportional to force of infection). Then since T is fixed (at 7.5 days), the quantity R 0 can be identified with parameter β. We computed time and test variable ascertainment rate, ρ(τ (t)) as well as time variable reproduction number, Re = (S/N )R 0 (see Figure 2 ). We also analyzed the association of estimated parameter values with outbreak measures, along with correlation between parameters, using Spearman's rank correlation (see Figure 5 and SI figure S19 for statistically significant results). In order to assess identifiability and confidence in our fitting procedure, we conducted uncertainty analysis. For the whole United States and its territories, the uncertainty quantification was carried out in the following manner:

(i) Simultaneously fit the model (2.1) to cumulative case and death totals, and national testing data as described above. (ii) Obtain the inferred fit daily case total curve, and under the assumption that the reporting error is normally distributed and relative in magnitude to the reported total at each data All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint 6 . . . point:

where g is the true number of daily cases andξ the set of true parameter values. We generated 10,000 cumulative case data sets with noise level of 60%, we then summed these daily case results and refit the model to each of them and the original death data simultaneously. For the results in table 2,3 and Figure 4 a value of s = .6 was used because this value causes the synthetic data-sets to cover the original data (see the SI figure 18 ). (iii) Arrange the generated values in increasing order and remove the top and bottom 2.5% in order to obtain the desired approximate 95% confidence intervals.

The resulting fit, confidence intervals, and average relative errors (see Figure 4 and table 3) indicate that the key model parameters are all practically identifiable assuming the 60% noise level in daily case reporting.

In addition, to further examine the relationship between ψ, R 0 , α and outbreak trajectories, we also highlight several states with different relative relationships between ψ and α representing the range of outbreak responses in the United States, and conduct both sensitivity analysis and counterfactual simulation studies with these examples (see Fig. 7 ). We fixed model parameters as those obtained from our fitting process, and varied R 0 ,ψ and α to examine their differential impacts on epidemic trajectory. The fit plots for all states and territories are included in the SI.

Our unified model simultaneously fits each state's cumulative mortality and reported case data, along with daily percent positive tests through our calibrated ascertainment rate, starting at first reported case to May 31, 2020 (as described in Methods). The estimated parameters, particularly quarantine (rate) factor, ψ, lockdown exit rate (fatigue), α, reproduction number, R 0 , and Infection fatality ratio, ξ, provide key information about the variability of COVID-19 spread and response in the United States. The resulting parameter values (table 2) and associated fits provide true (reported and unreported) cumulative case estimates for all 55 states and territories considered. These range in value from as low as 0.42% of population (Hawaii) to as high as 10.41% of population (New York State) over the time period, with an estimated 3.51% of the entire US population having been infected as of May 31, 2020 (Confidence Interval 2.72-4.43% population). These values represent both the wide range of observed outcomes and highlight the importance of incorporating testing data to obtain sufficient model flexibility (Fig. 2 , see also SI figures S15,S16 and S4-S16). Indeed, the ratio of true cumulative cases to reported cases across states ranged from as low as low as 2.99 (Rhode Island) to as high as 17.88 (the US Virgin Islands) with a value of 6.62 (CI 5.54-8.20), in line with CDC estimates [3] . Our fit IFR, both for every state and for the entire United states, works to counter both a common source of overestimation and a common source of underestimation [14] , [6] . Indeed, our gamma distributed delay (mean µ = 21) of days after infection until death and fitted true cumulative case totals with ascertainment rate ρ(t), together estimate IFR ranging from as low as 0.001 (Guam) to as high as .0156 (Rhode Island). There were outbreaks from mild to severe toll; with the fit IFR for the United states being 0.009 (CI .007 -.011) as of June 21 (cut-off date for deaths associated with infections on or before May 31), similar to estimates from a different retrospective study [14] .

To further dissect the relationships between model parameters and estimated quantities we next turn to correlation analysis. The following notable pairs of epidemiologically and statistically significant variables were found: I 0 vs date of first reported case (after Jan 21) with Spearman correlation = .48 (p = 2.3 × 10 −4 ); mortality (% population) vs. ψ with = −.53 (p = 3.7 × 10 −5 ; true Case Estimate vs. ψ with = −.68 (p = 2.1 × 10 −8 ); ψ vs. α with = .52 (p = 5.8 × 10 −5 ; peak daily cases vs. ψ with = −.72 (p = 1 × 10 −9 ; time to peak daily cases vs. R 0 with = −.34 (p = 1.6 × 10 −2 ; cumulative case estimate vs. IFR, = .33 (p = 1.6 × 10 −2 ); cumulative reported cases vs. IFR with = .71 (p = 4.7 × 10 −5 ); cumulative case estimate/reported cases vs. IFR with All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted January 6, 2021. = −.63 (p = 2.8 × 10 −7 . This analysis (see Figure 5 and SI figure S19 for all correlation plots) together with our three variable scatter plots in Figure 3 suggest that ψ dominates α and R 0 in terms of statistically significant correlation between model parameters and quantities of interest such as cumulative and peak daily case total as well as mortality and lockdown fatigue.

Furthermore we observe that ψ and estimated true case totals (shown in 5 (c)) as well as peak daily case totals ((d) of the same figure) follow an inverse proportionality relationship with respect to ψ. Indeed, as derived in [4] , for the case α = 0 (indefinite quarantine period), the final cumulative infected, C = 

where s∞ = lim t→∞ S(t) N and I 0 ≈ 0 (at start of outbreak). In each formula, the factor 1 1+ψ is multiplied by the corresponding classical relations for final and peak outbreak size. In view of our model fitting, the true cumulative and peak case estimates of each state fall roughly into All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in Figure 4 : Model fit together with 95% confidence intervals for United States Fit an inverse proportionality with ψ, although states with high lockdown fatigue (α much larger than zero) stray from this pattern, as shown in Figure 6 . In particular, by simply calibrating the inverse proportionality relation, C = K N Y 1+ψ , to the estimated parameters and cumulative cases of All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted January 6, 2021. New York (which has small ψ and α ≈ 0), we approximately fit the rest of the state values, with the analogous result for peak cases (see Fig. 6 ). The different values of R 0 and α in each state, however, shift the respective case numbers from the expected values predicted by inverse relation with ψ (Fig. 6 ). Thus while ψ has strongest impact, with rapidly escalating costs as the lockdown (rate) factor decreases below a critical value, R 0 and α still have an influence on outcome. We will examine these relationships further with counterfactual simulation studies. Also of note is the suggested relationships between R 0 and I 0 as well as R 0 and time to peak daily case incidence. Figure 5 panel a, which shows the aforementioned positive correlation between estimated I 0 and date of first reported case (in terms of days after Jan. 21), and figure in the SI, which shows the negative correlation between R 0 and time to peak daily cases. This is further supported by the plots of time variable Re in Figure 2 panel (a) where it can be seen that in general states which have a later date for their first reported case tend to have higher reproduction number (for specific values see table 2). Other relationships of note are the positive relationship between IFR and cumulative case estimate as well as IFR and reported case estimate as well as the negative relationship between IFR and ratio of true to reported cases. This first relationship is likely explained by demographics and the law of large numbers: with higher cumulative case totals the likelihood both that hospital systems are overwhelmed and more individuals and clusters in high risk categories will become infected increases (for all fit IFR values see table 2). The last two relationships likely explain each other: as more cases and deaths were recorded greater All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint . . . Figure 6 : Both Cumulative cases and ψ (top right) and peak cases (bottom right) and ψ follow the inverse proportionality relationship described in [12] . The plotted lines indicate change in cumulative (peak) case estimate as ψ is varied from 10% it's fit value up to the fit value itself. the color gradient of the lines from dark blue to dark red is indicative of increasing reproduction number R 0 This suggests that while ψ dominates these relationships, as our correlation analysis also indicates, R still has an impact on outcome with regard to both peak daily cases and cumulative case total. Top left shows the inverse proportionality relationship for the New York parameter set going from fit ψ to the maximum observed value of 6000. Bottom left is the same except for peak daily cases where the red lines follow equations (4.1) and (4.2) respectively effort was put into case reporting. And even in locales with successful outbreak responses such as Finland [11] , or Hawaii in our own model, we note that the proportion of cases undetected in the early states of the outbreak is very high. For our United States fit we estimated an IFR of 0.009 (95% confidence interval .007 to .011), which is in line with the observed IFR for countries with very high levels of testing such as South Korea [16] and Finland [11] .

Correlation analysis alone does not give the full picture. To further examine and differentiate the effects of α, ψ and R 0 on epidemic trajectories and suggest potential alternative strategies for management of subsequent waves of COVID-19 we carried out both counterfactual simulation studies and sensitivity analysis on these key model parameters. This was done using the fit model parameters for selected states which broadly speaking represent the range of outbreak responses which occurred in the United States. These are, in order of increasing cumulative outbreak size as of May 31 (see Figures 6, 7, 8) : As can be seen in these figures the states all reported their first confirmed case within a week of one another, but had significantly different outcomes. Taken as percent of population Hawaii had both the lowest cumulative case total across all states and territories and one of the briefest outbreaks in terms of increasing daily case totals. All other considerations aside response RLLF represents the optimal strategy in terms of both peak and cumulative daily case totals.

Tennessee highlights the importance of lockdown fatigue. In the baseline fit epidemic trajectory, despite having a lockdown factor higher than that of Hawaii, Tennessee's cumulative All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint Figure 7 : Representative examples of different relative value combinations of ψ and α, quality of strategy decreases from top to bottom row case totals are an more than 4 times as high as a result of significant linear growth in daily cases not far off the peak value. While this response will likely still keep hospitals from being overwhelmed provided sufficiently high ψ it both prolongs the outbreak, and leads to more cases than response RLLF. Georgia represents a number of states for which both the fit ψ and α values falls into a middle ground between the maximal and minimal fit values across all states and territories. As can be seen in the fit trajectory and table 2 despite have a ψ value an order of magnitude lower than that of Tennessee and similar dates of first reported case, the two states have relatively similar All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint . . . cumulative case totals, suggesting that there may be a critical threshold for α which is examined further below.

New York is representative of the final observed response, slow lockdown with low lockdown fatigue. Despite having both the highest reported cumulative and estimated true case totals across all states and territories New York's outbreak was essentially over at a time when many states were experiencing significant daily case totals despite similar dates of first reported case. All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted January 6, 2021. Our results suggest that the optimal response strategy in terms of case totals, and therefore also fatalities, remains rapid lockdown with low lockdown fatigue (response RLLF). Here the speed of lockdown is the critical factor as our model considers that individuals and/or government will eventually react to escalating cases, reaching similar numbers of self-quarantined (distinct for each state parameter set) in the counter-factual simulations of varying ψ with the difference being delayed lockdown resulting in order of magnitude more cases. While no substitute for the rapid lockdown strategy, the next best intervention may be sustained public social distancing and mask wearing, targeting transmission reduction rather than removing susceptibles all together, to reduce R 0 . This has the added benefit of reducing the number of individuals that need to be quarantined, with those that are able to self quarantine doing so preferentially throughout the entire outbreak, but with at least half the quarantined population not returning to normalcy for a period of 30 days (see figure insert) . Indeed, across all sample states, and regardless of model parameters, self-quarantine periods lasting longer than this threshold greatly reduce the daily case total on May 31, suggesting that if this is achieved other methods of managing the outbreak, such as contact tracing, could then be employed to manage subsequent cases.

Sufficiently high lockdown fatigue can lead to sustained linear growth in cumulative cases (or sustained steady state of daily cases) regardless of ψ value. Observe that despite Tennessee and Hawaii having fairly similar fit ψ values, both the daily case load on May 31 and the peak All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint daily case incidence for Tennessee are an order of magnitude higher than those of Hawaii when ψ is reduced by 90% (see Figure 8 ). Comparatively this difference is closer to a factor of five in baseline simulations. Figure 9 suggests that this is attributable to their different α values, which fall on opposite ends of the considered range. Thus our model suggests that lockdown fatigue, α, primarily affects whether or not daily case incidence on May 31 is close to zero relative to peak cases, where lockdown factor (or speed), ψ, modulates the peak level. A similar result is suggested with respect to R 0 and peak daily case total as well as date of peak, importantly with high α significantly dampening the impact of lowered R 0 (see Figures 8, 9) . While caution must be exercised with short duration or high-turnover lockdowns (responsive to accumulating cases), governments and individuals are also reluctant to return to strict lockdown. Given this reality, while response RLLF has the best outcomes, it is important to consider the viability of alternative strategies which might be implemented for future outbreaks.

In this study we fit a unified model to case, death, and testing data from all 50 states as well as Washington DC and four outlying territories in order to quantify impacts of a wide variety of lockdown entry/exit responses, which we correlate with true case estimates for insight into the dynamic COVID-19 outbreak control in the United States. Despite the heterogeneous nature of exhibited responses to the outbreak within each state, we were able to obtain quality fits for all considered states and territories. A crucial step in doing so was incorporating testing data, τ (t), by calibrating the ascertainment rate, ρ(τ (t)), as a saturating function of per-capita tests, which provided extra confidence in our fitting because daily percent positive tests were mimicked in each state. The inferred true cumulative case totals demonstrate high varying levels of under-reporting, with an estimated 6.62 cases per reported cases (cumulative) overall in the US up to May 31, 2020. Furthermore, our modeling framework allows us to infer Infection Fatality Ratio, found to be 0.9% for the US over time period. These are likely only estimates for (reported and unreported) symptomatic cases which do not account for a significant extent of asymptomatic individuals, since testing, reported case, and mortality data are all derived mostly from symptomatic infections.

Our model captured increasing case ascertainment rates through time (as testing increased), with a significant positive correlation to the amount of cases reported in a state, in contrast to a recent study in France showing a negative relationship between ascertainment ratio and reported cases as larger case counts hampered their substantial contact tracing program [11] . In the United States, comparatively little effort was put into contact tracing, and the particular geographical or political features of COVID-19 spread, along with higher case burden focusing community attention on the disease, may have contributed to the reverse correlation. Faced with a spreading outbreak, high numbers of undetected infection, and no nationwide control policy, the US adopted heterogeneous responses centered around some form of lockdown, which we characterize according to parameter fits of quarantine entry/exit rate for each state.

Our model fitting results suggest, broadly speaking, that the distinct responses can be divided into one of four categories (arranged in order of generally increasing cumulative case totals): Rapid Lockdown with Low Fatigue (RLLF), Rapid Lockdown with High Fatigue (RLHF), Intermediate Lockdown with Intermediate Fatigue (ILIF), and Slow Lockdown with Low Fatigue (SLLF). Using correlation analysis of the key parameters and outcomes over all US states, along with counterfactual simulation studies for representative states of these range of responses (Hawaii, Tennessee, Georgia, and New York), we translate the unified model fits into assessing the variable aspects of reactive lockdowns or self-quarantine (speed, scale, duration/fatigue) in combination with sustained interventions aimed at reducing R 0 , such as public face mask wearing and social distancing, contact tracing.

Lockdown (rate) factor ψ played the largest role in outbreak severity, in particular both cumulative and peak cases roughly followed an inverse proportionality relationship with ψ, All rights reserved. No reuse allowed without permission. perpetuity. preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint which was theoretically derived and applied to the epidemic in China [12] in the instance selfquarantine (lockdown) exit rate, α, is zero. Differences in R 0 also influenced outbreak shape, where its reduction mainly results in both a delayed peak in daily case totals and decreased size of peak daily case incidence, along with lessening the scale of lockdown needed to control the diseases. High lockdown fatigue (α significantly larger than zero) in some US states did impact their outbreaks though, the daily case load at the end of the time-frame considered (May 31) being large, with magnitude determined primarily by ψ, but also R 0 . Importantly, independent of all other model parameters our model suggests a critical threshold of 30-60 days for half-life return to normalcy, that is for 50% of the quarantined population to exit lockdown (see Figure 9 ). With this threshold indicating the approximate point at which there is significant reduction in daily case load within approximately 90 days of first reported case, perhaps allowing for other less invasive control methods, such as contact tracing, to be implemented.

Our results provide insights for management of future outbreaks. The optimal strategy purely in terms of cumulative, and peak daily cases, and therefore also fatalities remains response RLLF (rapid lockdown with low fatigue). However, in order to balance economic/social concerns in addition to the purely epidemiological, we suggest a new response to an outbreak compared to those observed in the United states first wave; Rapid measured lockdown with intermediate fatigue. This strategy describes implementation of rapid reactive lockdown as soon as possible in conjunction with subsequent wave being detected, lasting at least 30 days before 50% return to normalcy, with the underlying R 0 and public adherence to proper social distancing and mask wearing determining the scale of closures during lockdown (where at-risk or sectors responsible for large-scale spread are preferentially included in quarantine). The quick response and critical duration of quarantined sectors will allow case numbers to be sufficiently reduced after the 30 day period for contact tracing to be feasible and in combination with broader measures aimed at lowering R 0 (e.g. face masks) can prevent any substantial subsequent wave until effective vaccines are widely taken by the population.

Figs. S1-S19 Table S1 , S2 All rights reserved. No reuse allowed without permission.

perpetuity.

preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in

The copyright holder for this this version posted January 6, 2021. ; https://doi.org/10.1101/2021.01.04.21249231 doi: medRxiv preprint . . . 

Cdc covid-19 pandemic planning scenarios

Timing of state and territorial covid-19 stay-at-home orders and changes in population movement -united states

The COVID Tracking Project

Differential impacts of contact tracing and lockdowns on outbreak size in covid-19 model applied to china

Ecommunicating the risk of death from novel coronavirus disease (covid-19)

The effect of human mobility and control measures on the covid-19 epidemic in china

Prediction models for diagnosis and prognosis of covid-19: systematic review and critical appraisal

COVID-19 In the United States Project

An Ongoing Repository of Data on Coronavirus cases and deaths in the US

Underdetection of covid-19 cases in france threatens epidemic control, pre-print

Estimating the cumulative incidence of sars-cov-2 infection % and the infection fatality ratio in light of waning antibodies, pre-print

Estimating the case fatality ratio for covid-19 using a time-shifted distribution analysis, pre-print

Report of the who-china joint mission on coronavirus disease

How the united states flunked the covid-19 test: Some observations and

Acknowledgment JCM, CJB, and HG are supported by a U.S. National Science Foundation RAPID grant (DMS-2028728). HG was also supported by an NSF grant (DMS-1951759) and a grant from the Simons Foundation/SFARI(638193). CJB is partially supported by an NSF grant (DMS-1815095).The authors would like to thank Fadoua Yahia (Department of Mathematics, University of Louisiana at Lafayette) for her contributions to an earlier version of this paper.