Abstract
Countries across the world observed dramatic rises in COVID-19 cases and deaths in March 2020. In the United States, delays in the availability of diagnostic testing have prompted questions about the extent of unobserved community transmission. Using a simulation model informed by reported cases and deaths, we estimated that tens of thousands of people (median: 22,876, 95% posterior predictive interval: 7,451 - 53,044) were infected by the time a national emergency was declared. Our results also indicate that fewer than 10% of locally acquired, symptomatic infections in the US were detected throughout much of late February. These results point to the need for immediate, large-scale efforts to mitigate the impacts of SARS-CoV-2 on the US.
Main text
SARS-CoV-2 is a newly emerged coronavirus that is causing a global pandemic (1). The unprecedented spread of SARS-CoV-2 owes to its high transmissibility (2), pre-symptomatic transmission (3), and transmission by asymptomatic infections (4). An appreciable fraction of infections are asymptomatic (5), and many others result in mild symptoms that could be mistaken for other respiratory illnesses (6). These factors point to a potentially large reservoir of unobserved infections (7), especially in settings where capacity to test for SARS-CoV-2 has been limited (8). The United States is one such country in which limited testing has been a major concern, particularly as imported cases, and now local cases, have increased over time (9). Until February 27, testing criteria in the US were limited to close contacts of confirmed cases and those with recent travel to China (9). This means that any local infections resulting from an unobserved imported infection would have gone unnoticed. Community transmission occurred without notice while testing was still being rolled out (10, 11), albeit to an unknown extent.
Our goal was to estimate the extent of community transmission of SARS-CoV-2 in the US that occurred prior to its widespread recognition. Unlike other countries where testing and containment measures were pursued aggressively (12, 13), rollout of testing in the US was slow (9) and widespread social-distancing measures did not go into effect until several weeks after the first reported case (14, 15). Understanding the extent of community transmission has major implications for the effectiveness of different options for control (16) and for anticipating the trajectory and impact of the pandemic (17).
To estimate the extent of community transmission of SARS-CoV-2 in the US, we used a stochastic simulation model that combined importation and local transmission processes. We informed model parameters with estimates from other countries, where available (Table S1), and estimated values of two unknown parameters by fitting the model to data on local reported deaths in the US (18). To model importation, we simulated observed and unobserved imported infections based on the number and timing of imported cases reported in the US (19) and assumptions about the proportion of different infection outcomes (5, 20). To model local transmission, we used a branching process model informed by estimates of the serial interval and reproduction number of SARS-CoV-2 from Singapore (3). Due to aggressive containment efforts there (12), we considered our model to be a conservative representation of community transmission in the US. To relate our model’s predictions to US data on reported cases and deaths, we also simulated the timing of symptom onset (3), case reporting (18), and death (21), for simulated infections for which those outcomes occurred.
By March 12, there were a total of 1,514 reported cases and 39 reported deaths that resulted from local transmission of SARS-CoV-2 in the US. We used this information to estimate the probability of detecting imported symptomatic infections, ρtravel, by seeding our model with imported infections, simulating local transmission, and comparing simulated and reported local deaths. Under our baseline scenario, this resulted in a median estimate of ρtravel = 0.39 (95% posterior predictive interval: 0.15 - 0.90). Simulating from January 1, we obtained 22,876 (95% PPI: 7,451 - 53,044) local infections cumulatively in the US by March 12 (Fig. 1A). Due to the exponential growth posited by our model, 2,958 (95% PPI: 956 - 7,249) local infections were predicted to have occurred on March 12 alone (Fig. 1B). Had we performed a simple extrapolation of reported cases and deaths based on ρtravel, our estimate of cumulative local infections by March 12 would have been only 5,018 (95% PPI: 2,350 - 12,445). This suggests that detection of local infections was less sensitive than detection of imported infections.
These results derive from our baseline scenario and show A) cumulative and B) daily incidence of all local infections, including observed and unobserved. In B, black shows the median and gray shading shows the 95% PPI.
We estimated the probability of detecting local symptomatic infections, ρlocal, by comparing our model’s predictions of symptomatic infections to local case reports on a daily basis. Over the course of February, daily estimates of ρlocal decreased from our uniform prior down to a low of 6.4×10−3 (95% PPI: 2.4×10−4 - 4.8×10−2) on March 1, as increases in simulated local infections outpaced newly reported local cases (Fig. 2B, black). As testing increased in March (Fig. 2B, red), so too did reported cases (Fig. 2A, red) and daily estimates of ρlocal (Fig. 2B, black). By March 12, we estimated ρlocal to be 0.80 (95% PPI: 0.34 - 1.00). Between February 23 (low estimate of ρlocal) and March 8 (last day of verified testing numbers), our daily estimates of ρlocal were well correlated with daily numbers of tests administered (Pearson’s correlation, median: 0.57, 95% PPI: 0.48 - 0.65). Although these results are consistent with the possibility that testing might have improved case detection in March, they also indicate that case detection was likely very low at times in February when containment might have been feasible.
A) Local symptomatic infections predicted under the baseline scenario increased exponentially, whereas reported cases increased more sharply in March. B) Based on this, we estimated how the probability of detecting local symptomatic infections changed daily in the US. Black lines show the median and gray shading shows the 95% PPI.
Successful fitting of our model was demonstrated by its predictions of local deaths by March 12 (median: 33, 95% PPI: 9 - 74), which were consistent with the 39 reported (Fig. 3). Although we did not fit our model to deaths on a daily basis, 85.5% of the deaths predicted by our model occurred within the same range of days over which local deaths were reported (February 29 - March 12). This indicates that, collectively, our model’s assumptions about the timing of importation, local transmission, and delay between exposure and death are plausible. Deaths caused by COVID-19 often occur several weeks after exposure (22). Thus, our baseline model predicts that there will be a median of 395 (95% PPI: 125 - 948) additional deaths as a result of infections that occurred by March 12. Relative to deaths reported by then, this represents an increase by a factor of 12.2 (95% PPI: 7.03 - 21.3).
Our baseline model’s predictions were A) consistent with reported deaths through March 12 (dashed line) and B) indicate that many more deaths should be expected after then based solely on infections that occurred by March 12. Results to the right of the dashed line do not reflect additional deaths that would result from new infections occurring March 13 or after. The black line shows the median and gray shading shows the 95% PPI.
There are several limitations of our analysis that should be acknowledged. First, our results were, in some cases, sensitive to deviations from baseline assumptions (Supplementary Text). Although most parameter scenarios we explored resulted in similar cumulative infections, higher values of R0 and earlier importation resulted in estimates in excess of 100,000 (Fig. S4). Second, our branching process model assumes exponential growth, which could be affected by social distancing (23) or the buildup of immunity (24). Neither of those factors were likely to have had much influence on local transmission in the US before March 13, however. Third, our parameter assumptions were based on analyses of data collected outside the US. Similar information has proven useful for other pathogens, such as Zika and Ebola (25, 26), in past public health emergencies. Fourth, we did not make use of airline data to model importation (27), but future applications of our method could incorporate such data.
These limitations mean that results from our baseline scenario should be interpreted cautiously. Nonetheless, based on our sensitivity analysis, we conclude that unobserved SARS-CoV-2 infections in the US by March 12 likely numbered in the tens of thousands, and quite possibly in excess of 100,000. This result, considered together with extensive pre-symptomatic and asymptomatic transmission of SARS-CoV-2 (3, 4), suggests that the US was well past the possibility of containment by March 12. Other modeling work (16) suggests that the feasibility of containing SARS-CoV-2 is highly sensitive to the number of infections that occur prior to scale-up of containment efforts. Our estimate that fewer than 10% of local symptomatic infections were detected by surveillance for much of February suggests that a crucial opportunity to limit the impact of SARS-CoV-2 on the US may have been missed. Although the number of tests administered increased in March (9), so too did the number of infections and, consequently, the demand for testing.
Coincident with the March 13 declaration of a national emergency (14), social distancing measures went into effect across the US (15). Our estimate of several thousand active SARS-CoV-2 infections at that time suggests that large-scale mitigation efforts, rather than reactionary measures (28), are necessary. Even if these efforts begin to reverse increases in SARS-CoV-2 transmission in the US, our results show that a downturn in COVID-19 deaths may not appear until several weeks later. Analyses of the impact of large-scale mitigation efforts in China (29) provide reason for optimism that they can be effective.
Materials and Methods
We calibrated a stochastic model, including separate importation and local transmission steps, to two publicly available datasets on cases of COVID-19 internationally and in the United States. All code and data used are available at http://github.com/TAlexPerkins/sarscov2_unobserved.
Data
We obtained data on the number of imported cases and deaths from line list data compiled by the Models of Infectious Disease Agent Spread (MIDAS) Network (1). These data informed the number and timing of imported infections predicted by our importation model. We obtained data on the total number of US cases and deaths and total number of cases and deaths globally from time series compiled by the Johns Hopkins University Center for Systems Science and Engineering (2). These data informed our estimates of the proportion of local infections detected. We also used these data in an alternative importation scenario in which the timing of imported infections was sampled proportional to daily global incidence.
Importation model
We considered cases associated with international travel in the MIDAS dataset to be imported. We removed SARS-CoV-2-positive individuals who were repatriated from the Diamond Princess cruise ship from our analysis, due to the fact that they were quarantined (3), leaving 153 imported cases (including one death). We first estimated the number of imported infections based on the probability that an infection would be symptomatic, the probability of an imported symptomatic infection being detected, and the probability of death among symptomatic infections (case fatality risk, CFR). The CFR and the probability that an infection is symptomatic were drawn from beta distributions with parameters given in Table S1, with means of 2.29% and 17.9%, respectively. We jointly estimated the probability of detection of imported symptomatic infections, ρtravel, and the relative offspring number of asymptomatic infections, α, by running the importation and branching process models across a range of values of those parameters and calculating the probability of observing the number of reported deaths through March 12; this approach is described in more detail in the parameter calibration section below. The probability of the number of unobserved imported infections being between 0 and 20,000, along with the 152 observed cases and 1 observed death, was calculated using a multinomial distribution; the number of imported infections was then sampled from that distribution. We then smoothed the date of known imported infections with a Gaussian kernel and sampled dates of all imported infections from that distribution. As an alternative scenario, we distributed the timing of imported infections based on the timing of international incidence, with cases in China excluded after February 3, due to a ban on entrance by non-resident foreign nationals who had been to China within the past 14 days enacted on February 2. For each scenario and parameter combination, we generated 1,000 sets of imported infections.
Transmission model
We simulated local transmission in the United States from January 1 to March 12 using a branching process model, seeded by the aforementioned importation model. Each replicate draw of the number and timing of imported infections seeded one simulation of the branching process model, to maximally represent uncertainty in both importation and transmission processes. The number of secondary infections generated by each infection in the branching process model was drawn from a negative binomial offspring distribution with mean R and dispersion parameter k. Under our baseline scenario, we used a dispersion parameter of k = 1,000, approximating a Poisson distribution, due to a lack of estimates of k for SARS-CoV-2. Under alternative scenarios for k, we considered values of 0.15 and 0.30 to account for superspreading observed in outbreaks of SARS and MERS (4, 5). The number of secondary infections generated by asymptomatic individuals was also drawn from a negative binomial distribution, but with mean αR0, where α ∈ [0, 1]. Whether an individual was symptomatic was determined by a Bernoulli trial with probability equal to the proportion of infections that were asymptomatic in that replicate. Each secondary infection’s exposure time was drawn from a log-normal generation interval distribution with mean 4.56 days. In doing so, we assumed that the generation interval followed the same distribution as the serial interval.
In addition to exposure, we simulated three additional outcomes, and the timing thereof, in a subset of infections.
Symptom onset: The number of new symptomatic infections on day t was drawn from a binomial distribution with the number of trials equal to the number of infections with time of potential symptom onset on day t, and the probability of success equal to the proportion of infections that are symptomatic. For infections that were simulated to result in symptoms, the time of symptom onset was drawn from a Weibull incubation period distribution with mean 7.07 (6) and added to each individual’s exposure time.
Case reporting: The number of cases reported on day t was drawn from a binomial distribution with the number of trials equal to the number of infections with time of potential case reporting on day t, and the probability of success equal to the proportion of infections that are symptomatic. This accounts for the delay in reporting, but not underreporting, which is addressed below when we calculate the probability that a symptomatic infection is detected, ρlocal. The time of potential case reporting was drawn from a gamma distribution of the period between symptom onset and case reporting with mean 6 days, and added to each infection’s time of symptom onset.
Death: The number of deaths on day t was drawn from a binomial distribution with the number of trials equal to the number of infections that could have experienced death on day t, and the probability of success equal to the case fatality risk. The time of death was drawn from a log-normal distribution of time from symptom onset to death with mean 14 days (7), and added to each individual’s time of symptom onset.
All parameter values, and their associated distributions, are described in Table S1. Where parameter distributions were described in the literature using medians and interval measures of spread, we used the optim function in R to estimate parameters of those distributions that matched distribution moments reported by those studies. In that sense, all parameters in our analysis were treated as random variables, with associated uncertainty accounted for throughout our analysis. For the delay between symptom onset and case notification, we fitted a gamma distribution to data on the delay between symptoms and reporting for 26 US cases in the MIDAS line list data; the gamma distribution fitted the data better than negative binomial or log-normal distributions according to AIC (133.5, 134.6, and 134.0, respectively) (Fig. S1). Our mean estimate of 6.0 for this delay is in line with previous estimates from China of 5.8 by Li et al. (8) and 5.5 by Bi et al. (9). Three key parameters – R, the serial interval, and the incubation period – were taken from a single reference (6) to ensure that those estimates were consistent with each other. That is important because R and the serial interval jointly control the epidemic growth rate (10), so taking estimates of R and the serial interval from different studies could have led to unrealistic projections of epidemic growth rate.
The curve shows the maximum-likelihood fit of a gamma distribution (shape = 3.43, rate = 0.572) to those data.
We estimated how the probability of detecting locally acquired, symptomatic infections, ρlocal, changed over time. These estimates were based on the number of symptomatic cases reported each day, C(t), and our model’s predictions for the number of symptomatic infections that could have been reported each day, S(t), after accounting for a delay between symptom onset and reporting. We assumed a uniform prior for ρlocal, and on each day estimated a posterior equal to ρlocal (t) ∼ Beta(1 + C (t), 1 + S(t) − C (t)).
To understand how many deaths may occur after the time period of our analysis based on infections occurring through then, we set R0 = 0 from March 13 onwards and simulated our model forward to May 31. This allowed any infections occurring by March 12 enough time to result in death, for the proportion expected to result in that outcome.
Parameter calibration
Due to a lack of prior estimates for two parameters, we jointly estimated the proportion of imported symptomatic infections that were detected, ρtravel, and the relative infectiousness of asymptomatic infections, α. We fitted these parameters to the total number of deaths resulting from locally acquired SARS-CoV-2 infections in the US by March 12. To approximate a likelihood for given values of ρtravel and α, we simulated 200 replicate time series of imported infections, each based on the same value of ρtravel, and then simulated local transmission using the same value of α for each of the 200 replicates. For each of these 200 replicate simulations, we calculated the cumulative number of infections, ID, that, based on their timing, could have resulted in death by March 12. We then calculated the likelihood of the reported number of deaths, D, according to a binomial distribution in which D is the number of successes among ID trials that each have probability of success IFR, where IFR is equal to the CFR times one minus the probability of being asymptomatic. Each of the 200 replicates used independent draws from the uncertainty distributions of other model parameters, so we took the average of the 200 likelihoods to obtain a single marginal likelihood for a given value of ρtravel and α. After calculating this marginal likelihood across a grid of values between 0 (or 0.01 for ρtravel) and 1 in increments of 0.05 for each parameter, we smoothed this marginal likelihood surface using the bicubic.grid function in the akima package in R to create a gridded marginal likelihood surface with a 0.001 x 0.001 mesh. Finally, we drew samples from the posterior probability distribution of these parameters by resampling from this smoothed marginal likelihood surface, which implicitly assumed a uniform prior on the two parameters. We repeated this calibration procedure for each scenario that we explored, obtaining different estimates for ρtravel and α for each of our sensitivity analyses.
Sensitivity analysis
In addition to the alternative importation models, we also undertook a one-at-a-time sensitivity analysis for each parameter shown in Table S1, with the exception of the calibrated parameters (the last two rows). These last two parameters were re-calibrated as described in the previous section for each new parameter set and importation timing combination. Including the baseline scenario, there were a total of 18 scenarios (i.e., the baseline plus two explored values for each of seven parameters plus one additional scenario with different importation timing). For some parameter values explored in sensitivity analyses, we did not directly use literature estimates, but instead chose values which were plausible minima or maxima for that parameter; these are indicated by “lower” or “higher” in Table S1. For the dispersion parameter, we wanted to explore a value that allowed for superspreading but that generated less overdispersion than was observed for SARS; this formed our intermediate value in the sensitivity analysis. All baseline values were taken directly from literature estimates, with the exception of reporting delay, which was calibrated as described in the branching process model section. For that parameter, we obtained the low and high scenarios by multiplying the shape parameter by 0.5 and 1.5, respectively, while keeping the rate parameter the same. In this way, the reporting delay is the sum of one, two, or three identically distributed gamma random variables in the low, baseline, and high scenarios, respectively.
All time periods are given in days.
Supplementary Text
Imported infection predictions and estimates of imported case detection probability, ρtravel
By March 12, there were a total of 152 reported cases and one reported death in the US that were classified as imported on the basis of international travel to areas with known SARS-CoV-2 transmission (21). By jointly estimating ρtravel and the relative infectiousness of asymptomatic infections, α, we obtained a median estimate of 0.39 (95% PPI: 0.15 - 0.90) for ρtravel under our baseline scenario. This resulted in a median of 452 (95% PPI: 206 - 1068) imported infections. Under the alternative importation timing scenario, where importation timing was based on international case reports, we estimated ρtravel = 1.00 (95% PPI: 0.98 - 1.00) and 187 (95% PPI: 174 - 202) imported infections. An estimate of ρtravel = 1.00 implies that all symptomatic imported infections were detected, but it still means that asymptomatic infections would have been undetected.
Posterior predictive check against reported local cases
Using our estimate of ρlocal (t), we simulated the number of reported cases through time and compared this with the actual number of reported cases. By March 12, our model predicted that there should have been 1,530 (95% PPI: 475 - 3,496) reported cases, commensurate with the actual number of 1,514 reported cases (Fig. S2). As expected, this confirms that our estimates of ρlocal (t) were consistent with the model and the data.
The number of cases reported in the US compared to the number our model predicts were reported.
Sensitivity analysis of unknown parameters
Estimates of the proportion of imported symptomatic infections that were detected, ρtravel, and the infectiousness of asymptomatic infections relative to symptomatic infections, α, varied based on the values of the other parameters. In general, higher values for parameters expected to increase transmission (e.g., R) were associated with higher estimates of ρtravel (Table S2). Compared to a baseline estimate of ρtravel = 0.39 (95% PPI: 0.15 - 0.90) with R = 1.97, the estimate of ρtravel was 0.83 (95% PPI: 0.47 - 0.99) with R = 2.7 and 0.08 (95% PPI: 0.04 - 0.19) with R = 1.5. For a shorter serial interval with a mean of 4.7 days, the estimate was ρtravel = 0.52 (95% PPI: 0.19 - 0.96), and with a longer mean serial interval of 7.5 days, the estimate was 0.06 (95% PPI: 0.03 - 0.14). The estimated value of ρtravel was also lower if the CFR was low (ρtravel =, 95% PPI: 0.08 - 0.53), compared to the scenario with a higher CFR (ρtravel = 0.54, 95% PPI: - 0.96). Higher ρtravel estimates correspond to fewer undetected imported infections; therefore, fewer undetected importations are required to account for the observed number of local deaths through March 12 if the CFR is high, R is high, or the serial interval is short. In addition, when we based the timing of importations on international incidence (excluding China after travel restrictions were implemented on February 3) the estimate of ρtravel was 1.00 (95% PPI: 0.98 - 1.00) due to the increased probability of early importations – and more time for local infections to increase – under this scenario. There was greater uncertainty in our α estimates under most sensitivity scenarios, and in most scenarios the estimates of ρtravel and α were positively correlated (Fig. S3).
Samples (104) from the joint posterior distribution of the proportion of imported symptomatic infections detected (ρtravel) and the relative infectiousness of asymptomatic infections (α) under different parameter-sensitivity scenarios.
Median estimates and 95% posterior predictive intervals of the marginal distributions of proportion of imported symptomatic infections detected (ρtravel) and the relative infectiousness of asymptomatic infections (α) under different parameter-sensitivity scenarios.
Sensitivity analysis of cumulative infections
Because ρtravel and α were estimated for each parameter-sensitivity scenario, cumulative infections were relatively similar under the low, baseline, and high scenarios for many parameters. Cumulative infections were most sensitive to assumptions about R, the serial interval, and the timing of imported infections (Fig. S4, Table S3). The former two affect how quickly local infections increase, and the latter affects how much time they have to increase. Cumulative infections were also somewhat sensitive to assumptions about case fatality risk and the delay between exposure and death, because assumptions about those parameters influenced estimates of ρtravel and α, which were based on reported deaths.
Unlike other parameters, importation timing was not described in terms of simple numerical values; in that case, “mid” refers to our baseline assumption that the timing of unobserved imported infections followed the timing of observed imported cases, and “high” refers to the alternative scenario that their timing followed international incidence patterns.
Median estimates and 95% posterior predictive intervals of cumulative infections under different parameter sensitivity scenarios.
Sensitivity analysis of imported case detection probability, ρlocal
The proportion of symptomatic infections detected over time followed a similar pattern under all parameter sensitivity scenarios, with low values of ρlocal throughout late February followed by increases in March (Fig. S5). Long delays in case detection (9 days) were associated with the lowest proportion of symptomatic infections detected; in that scenario, ρlocal exceed 10%. mostly did not
Each panel represents a different parameter-sensitivity scenario.
Sensitivity analysis of the ratio of deaths after and before March 12
The ratio of deaths expected March 13 and after, relative to before then, was higher with changes in parameters that resulted in faster growth in local infections and later arrival of imported infections (Fig. S6, Table S4). The proportion of deaths expected to occur after March 12 also increased with increases in the delay between symptom onset and death (Table S4). Overdispersion (lower k) did not drastically alter our estimates of ρtravel or α (Table S2) or the number of cumulative infections (Table S3), but it did extend the lower and upper bounds on the range of the ratio of deaths after and before March 12 (Table S4).
Unlike other parameters, importation timing was not described in terms of simple numerical values; in that case, “mid” refers to our baseline assumption that the timing of unobserved imported infections followed the timing of observed imported cases, and “high” refers to the alternative scenario that their timing followed international incidence patterns.
Median estimates and 95% posterior predictive intervals of the ratio of deaths after and before March 12 under different parameter sensitivity scenarios.
Acknowledgements
Thanks to Jason Rohr and Moritz Kraemer for feedback on the manuscript.
References
References
Subject Area
Reviews and Context
0
Comment
0
TRIP Peer Reviews
0
Community Reviews
0
Automated Services
4
Blogs/Media
Author Videos