Chopping the tail: How preventing superspreading can help to maintain COVID-19 control.

Disease transmission is notoriously heterogeneous, and SARS-CoV-2 is no exception. A skewed distribution where few individuals or events are responsible for the majority of transmission can result in explosive, superspreading events, which produce rapid and volatile epidemic dynamics, especially early or late in epidemics. Anticipating and preventing superspreading events can produce large reductions in overall transmission rates. Here, we present a stochastic compartmental (SEIR) epidemiological model framework for estimating transmission parameters from multiple imperfectly observed data streams, including reported cases, deaths, and mobile phone-based mobility that incorporates individual-level heterogeneity in transmission using previous estimates for SARS-CoV-1 and SARS-CoV-2. We parameterize the model for COVID-19 epidemic dynamics by estimating a time-varying transmission rate that incorporates the impact of non-pharmaceutical intervention strategies that change over time, in five epidemiologically distinct settings-Los Angeles and Santa Clara Counties, California; Seattle (King County), Washington; Atlanta (Dekalb and Fulton Counties), Georgia; and Miami (Miami-Dade County), Florida. We find that the effective reproduction number (RE) dropped below 1 rapidly in all five locations following social distancing orders in mid-March, 2020, but that gradually increasing mobility starting around mid-April led to an RE once again above 1 in late May (Los Angeles, Miami, and Atlanta) or early June (Santa Clara County and Seattle). However, we find that increased social distancing starting in mid-July in response to epidemic resurgence once again dropped RE below 1 in all locations by August 14. We next used the fitted model to ask: how does truncating the individual-level transmission rate distribution (which removes periods of time with especially high individual transmission rates and thus models superspreading events) affect epidemic dynamics and control? We find that interventions that truncate the transmission rate distribution while partially relaxing social distancing are broadly effective, with impacts on epidemic growth on par with the strongest population-wide social distancing observed in April, 2020. Given that social distancing interventions will be needed to maintain epidemic control until a vaccine becomes widely available, "chopping off the tail" to reduce the probability of superspreading events presents a promising option to alleviate the need for extreme general social distancing.

Introduction

In the face of emerging epidemics with limited pharmaceutical options for treatment and prevention of infection, non-pharmaceutical interventions such as social distancing are critical for slowing epidemic growth. Shelter-in-place and other social distancing orders have helped to slow the pace of the COVID-19 pandemic, reducing the effective reproduction number – or the number of secondary infections produced by each infected person – to one or below in most places. In doing so, social distancing has effectively kept most regional healthcare systems operating under maximum capacity. However, after only a few weeks of declining numbers of daily cases due to an at or below one in the spring of 2020, most state and county governments in the United States began relaxing social distancing orders, citing their major economic impacts. In order to avoid epidemic resurgence, it is vitally important that governments employ long-term strategies that maintain epidemic control while economic reopening commences.

One obstacle to designing effective long-term strategies is a notoriously heterogeneous transmission process. It is now widely recognized that the minority of infections generate the majority of secondary cases, leading to the so-called 20–80 rule in epidemiology (the rule-of-thumb that 20% of people generate 80% of cases) (Woolhouse et al., 1997). Work on SARS-CoV-1, measles, and other respiratory viruses suggests that this skew in secondary cases is even larger (Lloyd-Smith et al., 2005). This heterogeneity gives rise to events in which a single infected person transmits a disease to dozens or hundreds of people – called superspreading events – which have played an important role in the COVID-19 pandemic (Adam et al., 2020, Kumar et al., 2020, Liu et al., 2020, Althouse et al., 2020, Hébert-Dufresne et al., 2020). Indeed, the frequency of asymptomatic and presymptomatic transmission, potential disconnect between infection and clinical presentation (Davies et al., 2020), and potential transmission via direct contact, aerosols, and surfaces (World Health Organization, 2020, Guo et al., 2020) are all features of SARS-CoV-2 that tend to promote superspreading. As local and national governments search for viable alternatives to shelter-in-place, a critical question is how effective curtailing superspreading events could be in controlling epidemic spread.

Practically, one strategy to help prevent superspreading is to prohibit medium to large indoor gatherings such as exercise classes, sporting events, concerts, and weddings for an extended period after allowing smaller and lower-risk activities to resume. From a modeling standpoint, predicting the effects of this straightforward intervention is difficult for two reasons: (1) local epidemiological dynamics are changing with evolving intervention strategies; and (2) information may not be available to parameterize detailed models of disease spread through heterogeneous populations. Despite these difficulties, it is important to consider some individual-level heterogeneity in transmission because model analyses of mean transmission rates alone may over-estimate the effectiveness of interventions, overlook potentially effective interventions that act on the heterogeneity within populations, miss potentially explosive resurgences, and poorly predict the final epidemic size (Lloyd-Smith et al., 2005, Hébert-Dufresne et al., 2020).

Studies of superspreading often empirically estimate secondary case distributions from recorded transmission chains and/or using branching process models (Lloyd-Smith et al., 2005, Liu et al., 2020, Althouse et al., 2020, Zhang et al., 2020b). These studies estimate a dispersion parameter, , that describes the variance in secondary cases based on a Negative Binomial distribution, where smaller values indicate more heterogeneity and skew and large values approach a Poisson distribution. Estimated values for SARS-CoV-2 remain uncertain, but are thought to range from 0.04–0.3 (Althouse et al., 2020, Hébert-Dufresne et al., 2020, Zhang et al., 2020b, Endo et al., 2020a), similar to the estimate of 0.16 for SARS-CoV-1 (Lloyd-Smith et al., 2005), which we use for this analysis. These empirical and branching process approaches are ideal for characterizing heterogeneity in secondary cases, but not for projecting epidemic trajectories through time, without being further embedded in a compartmental or network modeling framework.

Here, we present a mechanistic susceptible, exposed, infectious, removed (SEIR) model that uses data on cases, deaths, and mobility for parameter estimation, incorporates heterogeneity in transmission rates, and is realistic enough to be useful for scenario exploration but simple enough to be adapted to a wide range of settings. The key innovation in our model is in using the average of Gamma-distributed individual transmission rates at each time step, as supported by previous work on secondary case distributions, to generate the distribution of population-average transmission rates. This formulation allows us to both generate more realistic variation in trajectories than SEIR models that assume a single average transmission rate, and explore and quantify the impact of altering individual-level transmission distributions on population-level dynamics without more detailed information on contact networks, age structure, or other social information.

The model, with accompanying open-access code, can be used to fit to any county in the U.S. using publicly available data; here we focus on five contrasting epidemiological settings—Seattle (King County), Washington; Los Angeles (Los Angeles County), California; Santa Clara County, California; Atlanta (Dekalb and Fulton Counties), Georgia; and Miami (Miami-Dade County), Florida. For each location we estimate a time-varying effective reproduction number, , which represents the average number of secondary infections produced by each infected person, and is an important (though imperfect (Hébert-Dufresne et al., 2020)) metric of epidemic control. Using each fitted model, we truncate the individual-level transmission rate distribution and stochastically simulate epidemic dynamics into the future, representing a scenario where high-risk events are eliminated but smaller and lower-risk activities are allowed to resume. This model construction focuses on reducing superspreading events in which a single individual infects many others over a short time period (e.g., due to simultaneous contact with many susceptible people); it does not address situations where individuals act as superspreaders because they have a longer-than-average infected duration. We investigate the absolute impact of this superspreading prevention strategy on epidemic control, and compare its impact on epidemic dynamics (and ) to test-and-isolate and strong social distancing (e.g., shelter-in-place) interventions. Using this comparison we highlight non-pharmaceutical intervention strategies that are less extreme than shelter-in-place but still expected to reduce both epidemic growth (i.e., keep below one) and the probability of explosive resurgence.

Methods

Model structure

We developed a stochastic compartmental model using an SEIR (Susceptible, Exposed, Infectious, Recovered) framework to modelCOVID-19 transmission, which was first described in Childs et al. (2020). Our model divides the population into the following classes: susceptible (S); exposed but not yet infectious (E); infectious and presymptomatic (IP), asymptomatic (IA), mildly symptomatic (IM), or severely symptomatic (IS); hospitalized cases that will recover (HR) or die (HD); recovered and immune (R); and dead (D). We assume an underlying, unobserved process model of SARS-CoV-2 transmission described by Eqs. (1)–(10) and shown in Figure S1, where each term denotes the transition from compartment to . Eqs. (11)–(18) describe in detail the stochastic rates used to approximate the transition terms in Eqs. (1)–(10). Transitions between compartments are simulated as binomial () or multinomial () processes to account for competing rates (implemented with the rbinom and reulermultinom functions, respectively, in the pomp package (King et al., 2016)). We use an Euler approximation of the continuous time process with a time step of 4 h. To produce more realistic latent and infectious periods we divide each infectious class and the exposed period into multiple sub-stages (not shown in the equations below), which results in Erlang distributed periods within stages (Anderson and Watson, 1980, Lloyd, 2001). Specifically, we use three sub-stages for the exposed class, seven sub-stages for the asymptomatic infectious class, two sub-stages for the presymptomatic infectious class, five sub-stages for the mildly symptomatic infectious class, and five sub-stages for the severely symptomatic class. We translate durations into rates for our model with sub-classes and a Euler approximation using the method described in He et al. (2010).

By including asymptomatic and presymptomatic individuals, we are able to track “silent spreaders” of the disease, which have been shown to contribute to COVID-19 transmission (Ferretti et al., 2020, Li et al., 2020b). Mildly symptomatic cases are defined as those people that show symptoms but do not require hospitalization. We assume that all severely symptomatic cases will eventually require hospitalization and that no onward transmission occurs from hospitalized individuals. The time-varying transmission parameter, , describes the average per capita rate of contact between susceptible and infectious people at time , multiplied by the per-contact transmission probability. We modeled as a function of human movement using the scaling function: which treats as an exponentially decreasing function of physical distancing (; on a scale of 0–1 where 0 is no physical distancing, and 1 is maximum physical distancing). Here, is the estimated minimum possible transmission rate given minimal human movement (i.e., maximal physical distancing) and thus contact rate, while is the estimated transmission rate given maximal human movement (i.e., no physical distancing). To measure physical distancing we use SafeGraph’s “Shelter in Place Index” (Shelter in Place Index, 2020), which calculates the proportion of cell phone devices that are staying home.

To model individual heterogeneity in SARS-CoV-2 transmission rate, we allow individuals to vary over time by modeling an individual’s transmission rate in each time step as a Gamma distributed random variable with a dispersion such that the sum of an individual’s transmission rates over the duration of their infection approximates a Gamma distributed random variable with dispersion equal to previous Negative Binomial parameterizations for SARS-CoV-1 () (Lloyd-Smith et al., 2006), which closely approximates estimates of overdispersion for SARS-CoV-2 (Endo et al., 2020a, Grantz et al., 2020). Because we model the transmission rates as the multiplication of contact rate and infection probability, this heterogeneity implicitly considers both variation among individuals in infectiousness and contact rate, and can be thought of as modeling superspreading periods or events—windows in time when an infected individual has a particularly high transmission rate. Some other forms of superspreading, such as when individuals have unusually long infectious durations, are not captured by this formulation. Additionally, this formulation assumes that infectious individuals’ transmission rates are independent draws, while events that involve large aggregations of people inherently imply dependence. For relatively small superspreading events likely to be driven by one infectious individual, this independence assumption will have no effect. However, our model may underestimate the occurrence of the most extreme superspreading events, which could be characterized by multiple superspreading events occurring in the same 4h time window and would require multiple independent samples in the tail of the individual-level gamma distributions. To incorporate individual variation into an average time step, we model the population level transmission rate as the average of the transmission of all infected individuals at time . To do so we apply the property of Gamma distributions that the mean and variance of samples from a Gamma distribution with defined rate and scale is itself a Gamma distribution with mean equal to that of the original Gamma distribution and variance equal to the variance of the original Gamma divided by . A full derivation of the equivalence between the individual time step transmission rate distributions (which we will hereafter refer to with ), the individual infectious period transmission rate distributions, and the population-level transmission rate distribution is presented in the Appendix A.

We assume that observed deaths are a Negative Binomial random variable with a mean equal to total new deaths accumulated over the observation period (i.e., one day for this analysis), and a dispersion parameter that we estimate. We also assume that daily observed cases are a Negative Binomial random variable, but have a mean equal to the daily accumulated number of new symptomatic infections multiplied by a daily detection probability that we estimate from the data. We model daily detection probability as a monotonically increasing logistic function: where is the maximum fraction of symptomatic cases detected, is the logistic growth rate, and gives the location of the inflection point (where the probability of detection equals one half of the maximum detection probability, ). Because can be estimated to be in the future, the probability of detection of an infected case in the present can be any value between 0 and . We estimate newly observed cases to be a fraction of all new symptomatic infections at time . Though this ignores testing asymptomatic infections, any detection of asymptomatic infections will be captured as a higher estimated detection fraction of symptomatic infections.

Fitting the model

We used COVID-19 death and case data from The New York Times, based on reports from state and local health agencies (available at https://github.com/nytimes/covid-19-data). Using these data, which are available for all counties in the US, and any form of human movement data that can be scaled to 0–1, our model can be used to fit infection dynamics in any county.

We sampled over uncertainty in a series of parameters (Tables 2, 3) for which there was large uncertainty in the literature. To do so, we drew 800 sobol sequences, an efficient method for sampling input parameters (Burhenne et al., 2011), across a range of plausible values for each parameter in order to form 800 plausible parameter sets. For computational efficiency we assumed point estimates for other parameters (Table 1). These parameters most notably include the average duration individuals spend within infectious categories, for which some consensus has been reached in the COVID-19 modeling literature. While the impact of changes in these assumed means would be somewhat dampened because of our use of Erlang distributed periods (which leads to variation in transition times), more work could be done to precisely identify the sensitivity of model results to these point estimates. We also assumed equivalence in the relative infectiousness of each infected category except completely asymptomatic infections, given a lack of any published estimates for these parameters. Our range of plausible epidemic start dates (Table 3) spanned a period from the first of the year to the last day before the first recorded case in each location.

Table 2

Parameters for which we used a range of values and those ranges. Parameter ranges are assumed to not vary by location.

Parameter	Lower bound	Upper bound	Description	Estimates and sources
κ_A	0.4	0.8	Relative infectiousness of asymptomatic infections	0.6 (Perkins et al., 2020), few direct estimates, but many examples of asymptomatic transmission potential less than, but potentially close to that of, symptomatic infected individuals (He et al., 2020b, Ye et al., 2020, Hu et al., 2020)
α	0.3	0.5	Proportion of infections that are asymptomatic	Mean of 43.3% (Lavezzo et al., 2020), 44% (Gudbjartsson et al., 2020)
δ	0.1	0.3	Fatality rate among hospitalizations	Demographic weighted average that will vary by location, see Verity et al. (2020), 5% (Li et al., 2020a)
ρ_D	13 days	20 days	Time from hospitalization to death	Mean of 16 days (Gaythorpe et al., 2020), mean of 17.8 days (Verity et al., 2020), highly variable by region (Rees et al., 2020)
1−μ	0.025	0.075	Proportion of symptomatic infections that require hospitalization	Demographic weighted average that will vary by location, see Verity et al. (2020)

Table 3

Location-specific parameter range estimates. Population sizes obtained from the US Census Bereau (County Population Totals 2010-2019, 2019).

Parameter	Santa Clara county, CA	Los Angeles County, CA	Miami-Dade County, FL	King County, WA	Fulton and DeKalb Counties, GA
Epidemic start date	01-Jan–05-Feb	01-Jan–31-Jan	01-Jan–16-Mar	01-Jan–04-Mar	01-Jan–Mar-07
Population size	1,927,852	10,039,107	2,716,940	2,252,782	1,755,830

Table 1

Parameters for which we used point estimates and their values. These parameters are assumed to not vary by location.

Parameter	Value	Description	Estimates and sources
κ_P, κ_M, κ_S	1	Relative infectiousness of presymptomatic, mild symptomatic, and severe symptomatic	Assumed
γ	3.5 days	Preinfectious period	One meta-analysis (McAloon et al., 2020) found a mean incubation period of 5.8 days, another found a median of 5.1 days (Lauer et al., 2020). We use a shortened duration because we assume 2 days of presymptomatic transmission
λ_P	2 days	Presymptomatic duration	Range of 1–3 days (Wei et al., 2020), mean of 3.8 days (Zhang et al., 2020a), viral shedding estimated to begin 2.3 days prior to symptom onset (He et al., 2020a), many articles find presymptomatic infection is likely but do not estimate duration (Zhang, 2020, Gao and Li, 2020, Arons et al., 2020, Lee et al., 2020)
λ_A	7 days	Infectious period for asymptomatic infections	Mean seroconversion after 7 days (Wölfel et al., 2020)
λ_S	5 days	Time from symptom onset to hospitalizations (severe cases)	Mean of 5.5 days (Sanche et al., 2020), median of 4 days (To et al., 2020), mean of 5.6 days (Verity et al., 2020)
λ_M	5.0 days	Time from symptom onset to recovery (mild cases)	Infectiousness based on viral shedding estimated to decline substantially within 7 days (He et al., 2020a, Wölfel et al., 2020) *Note, He et al. (2020a) takes samples from hospitalized patients; we assume similar viral shedding in mild infections
ρ_R	13.3 days	Time from hospitalization to recovery	Mean of 13.3 (Tindale et al., 2020), mean of 20.51 (Verity et al., 2020), highly variable by region (Rees et al., 2020)

For each parameter set we used the package pomp (King et al., 2016) in the statistical programming language R (R Core Team, 2020) to estimate the following parameters: and , transmission rate parameters (see Eq. (19)); , number of exposed individuals that initiate the epidemic; , m, and , maximum, slope, and inflection point day of the sigmoidal case detection function; , Negative Binomial dispersion parameter for deaths; and , Negative Binomial dispersion parameter for cases. We assume that all individuals are initially susceptible except for those who are exposed (). We fit all parameters to daily deaths, cases, and mobility in four steps. First, for each of the 800 parameter sets from the sobol sequences, we used the mif2 function in pomp to perform iterated filtering, beginning with random starting conditions for the remaining parameters that were fit and using 90 filtering iterations with 4000 particles. We then continued to fit the top 20% (i.e., 160) parameter sets with the highest log likelihoods for an additional 200 iterations using 4000 particles. Finally the top 10% of those parameter sets were filtered twice more using 200 iterations at 4000 particles. At the end of each of the four filtering steps, log-likelihoods were calculated ten times using 5000 particles to produce standard errors. From these four filtering sets, the top ten parameter sets by log likelihood were selected. In total, each county took approximately 24 h to fit using twenty cores.

We focus on and refrain from drawing conclusions directly about, for example, the estimated case detection probability (), the effectiveness of physical distancing (), or infection fatality ratio (IFR), because of some identifiability issues due to correlated parameters (Figures S2–S6). We calculated at each time as estimated times the median proportion of the population remaining susceptible on each day times the average infectiousness over an infection (as defined by our model structure). For model fits and parameter estimates (e.g., ) shown in results figures, we combine two sources of uncertainty: parameter uncertainty, which we include by using the 10 parameter sets with the largest log likelihoods; and stochastic epidemic trajectory variation, which we include by simulating 500 epidemics for each parameter set. Taken together, this uncertainty captures the widest plausible range of epidemic dynamics given our model and data. To avoid implausibly wide variation we discard simulated epidemics that did not reach at least a total of 500 infected individuals.

Simulating epidemics under interventions

Any intervention type, intensity, or duration can be modeled using this framework and open-source code (available at https://github.com/morgankain/COVID_interventions) given that it can be written as a function that modifies either human movement or (e.g., social distancing or a pharmaceutical intervention that reduces the probability of infection). Here we consider interventions that reduce the skew of the individual time step transmission rate distribution (), and thus the average time-varying transmission rate ; this is our mathematical representation of reducing highly infectious contact periods or events, which for COVID-19 tend to occur in crowded enclosed environments (e.g., church choirs and exercise classes). Specifically, we model truncation of the distribution by assuming that all samples within the top X% of the distribution are resampled. In particular, we consider truncation strategies that act in addition to moderate background social distancing measures. In doing so, we model the effects of social distancing as reducing the mean transmission rate (Eq. (19)) which still allows for individual heterogeneity in the resultant transmission rate. Truncation then operates to remove the highest transmission rates after considering background social distancing. Because of the skew of this distribution, truncating just the upper 0.1% of the distribution can have a profound effect on the mean transmission rate (Figure S7). To visualize the dynamics of interventions, for each location we simulate 500 epidemics from the maximum likelihood estimate across the 800 parameter sets. The uncertainty band we plot represents the central 95% range of outcomes seen across all stochastic realizations that resulted in epidemics for this parameter set, and thus should not be taken as representation of uncertainty in parameter values or model structure.

Results

Epidemic trajectories

The model produced realistic fits to five contrasting epidemiological settings—King County, Washington; Los Angeles County, California; Santa Clara County, California; Fulton and Dekalb Counties, Georgia; and Miami-Dade County, Florida (hereafter, Seattle, Los Angeles, Santa Clara County, Atlanta, and Miami; Fig. 1A). Among these locations, we estimated that prior to interventions, ranged between approximately 2 and 4 (Fig. 1B). We also estimated that dropped below one following initial shelter-in-place orders in all counties, though only briefly in some locations. In particular, in Miami, Los Angeles, and Atlanta climbed above 1 by mid-May, while in Seattle and Santa Clara County remained below 1 until early June. After approximately another month of rising , increased social distancing led to a second decline in , which once again dropped below one in all locations by mid August (Fig. 1C).

Fig. 1

Model estimated daily deaths and reported cases (A), and reproduction number (B, C) for five locations: Atlanta (red), Seattle (green), Miami (gold), Santa Clara County (blue), and Los Angeles (purple). Los Angeles is displayed on a different -axis due to differences in magnitude of reported deaths and cases. For each county, we show the 10 model fits with the best log likelihoods. Panel C show the same results pictured in B, but are zoomed in to April 15–Aug 14 to better show the dynamics around . Black points are observed daily deaths (A, top row) and reported cases (A, bottom row) in each county. Solid lines display median of model simulated trajectories (A) and mean of 7-day smoothed (B, C). Ribbons overlay the central 95% of all stochastic simulations for each of the 10 best-fit parameter sets for epidemic dynamics (A) and full range of estimated (B, C) over the 10 best-fit parameter sets. We highlight the large differences in stochastic epidemic trajectories for the single best parameter set in Figure S8. Vertical axes in panel A are pseudo-log transformed for visibility. Goodness of fit (R2 values) and log likelihoods for each of these fits are given in Figure S9 and Table S1 respectively . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

From the ten best fitting parameter sets in each location, our model estimated moderate variation in infection fatality rate (IFR), case detection probability, and other model parameters among locations (Figure S10, Figure S11); however we caution against interpreting these parameter estimates in isolation because of their correlation with other parameters (Figures S2–S6). While our model is able to recreate plausible epidemic dynamics and sensible estimates for (Fig. 1), model fitting indicated that there are a number of plausible parameter combinations that produce epidemic dynamics that match the observable data. As data continues to increase and more is learned about the disease, continued model fitting may be able to eliminate some of these alternative possibilities.

Interventions

As a basis for comparison, focusing on just three locations – Los Angeles, Miami, and Seattle – if all social distancing was simply to end, a rapid increase in cases and deaths would be inevitable (Fig. 2, blue). However, non-pharmaceutical interventions, including continuing strong social distancing, infected isolation with intermediate levels of social distancing, or averting superspreading with intermediate levels of social distancing are capable of limiting epidemic growth (Fig. 2) and keeping near or under 1. Here, we consider intermediate levels of social distancing that correspond to mobility levels that are an average of baseline mobility prior to the first wave of shelter-in-place orders and final mobility levels observed in the last week of data. Either an infected isolation strategy that reduces to intermediate levels of social distancing and catches 90% of all mild and severe cases of COVID-19 before they are able to transmit (Fig. 2, green), or a truncation strategy that similarly reduces to intermediate levels of social distancing and further removes the top 1% of the individual time step transmission rate distribution () with 75% efficiency (Fig. 2, purple) are able to suppress epidemic growth (and reduce to below one) in all locations.

Fig. 2

Maintaining strong social distancing (orange), test-and-isolate with intermediate social distancing (green), or superspreading aversion with intermediate social distancing (purple) are necessary over long periods to prevent a major epidemic resurgence (blue) in each location where we fit our model. Daily deaths are shown in (A - C) and daily reported cases are shown in (D - F). Continuing social distancing at the levels seen as of August 14 (orange) will lead to nearly steady daily reported cases and deaths in Los Angeles, CA (A, D), and slow declines in cases and deaths in Miami, FL (B, E) and Seattle, WA (C, F) due to estimated very near one. For both test-and-isolate and truncation interventions we assume an intermediate level of mobility (an average of baseline mobility prior to the first shelter-in-place orders and final mobility levels observed in the last week of data). Bands show the central 95% of stochastic simulations of daily cases and deaths for the single maximum likelihood parameter set and solid lines show the median among simulations. Dates range from February through October of 2020. Vertical axes are pseudo-log transformed for visibility. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Curtailing superspreading

Limiting opportunities for superspreading by “chopping off the tail” of the contact rate or infectiousness distributions can be highly effective at epidemic control (Fig. 2), driving epidemic growth to be negative and bringing the average number of secondary cases () below 1. An example truncation intervention is illustrated in Fig. 3: because the individual transmission rate distribution, , over a 4-hour time period is so skewed (Fig. 3A; see Appendix A for derivation), truncating the upper 0.1% yields a large reduction in the mean and a moderate reduction in the variance of the population-level average transmission rate (Fig. 3B, shifting from red to blue distribution). A variety of possible truncation strategies exist, including eliminating varying proportions of (e.g., upper 0.5%, upper 1%) with varying levels of efficiency (ranging from 50%–100%) (Figure S12).

Fig. 3

Example of how truncating the individual-level transmission rate distribution, , (A) affects the population-average transmission rate (B), and combinations of pSIP and truncation strategies that reduce to one in a fully susceptible population (C). The three panels in C show the combinations of truncation and pSIP that produce an of one for three levels of truncation efficiency. (A) Truncation at the upper 0.1% of (sampled over a 4-hour time step), in which truncation occurs with 100% efficiency for values above the dashed line. (B) Resulting effect on the population-level average infection rate when there are 1000 infected people currently in the population, where the original distribution is in red and the truncated distribution is in blue. The distribution is shown over 10,000 simulations for a population characterized by an individual reproduction number distribution with mean of 2.5 and overdispersion parameter, . Horizontal and vertical axes in A and B are square root transformed for visibility. In C, the triangles show baseline pSIP in each location and circles show max pSIP reached during social distancing. Solid lines indicate the mean over the ten best fits, and the ribbon is the full range of estimates from these fits. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

An alternative measure of the impact of averting superspreading (i.e., truncation interventions) is how much general social distancing can be avoided by instead truncating the transmission rate distribution. Prior to social distancing orders, the estimated proportion sheltering in place (the proportion of cell phone devices staying home; hereafter, pSIP) ranged from 20%–25% across our focal locations (Fig. 3C, triangles), indicating the baseline level of mobility. By increasing pSIP in combination with truncation interventions, a variety of combinations are predicted to provide epidemic control (for example, by reducing transmission rates such that in a fully susceptible population would be 1; Fig. 3C). If the truncation intervention is 100% effective, truncating only approximately the upper 0.15% of individual transmission rates, , (Fig. 3A) is effective enough to maintain transmission rates such that would be 1 in a fully susceptible population, while allowing mobility levels to return to baseline (Fig. 3C). Alternatively, if truncation interventions are only half as effective, the same 0.15% truncation intervention would require moderate-strong social distancing (pSIP from 30%–45%; Fig. 3C). The nonlinear effects of social distancing and truncation on transmission make the combination of interventions needed to maintain epidemic control sensitive to the efficiency and strength of each intervention mechanism.

Superspreading and epidemic resurgence

Even if the epidemic is brought almost entirely under control (e.g., to within 1–5 infected individuals remaining in the population), epidemic resurgence remains a possibility if interventions wane, allowing to increase above one. As we show in (Fig. 3), many different combinations of pSIP and truncation can be used to produce the same (in Fig. 3, an of 1); however, epidemic dynamics will vary by combination because of the variation in individual time-step transmission rates, . If rises above one because interventions are relaxed, the specific combination of pSIP and truncation that remains in place will determine the resulting dynamics. Here we examine how different truncation interventions will affect epidemic extinction probability and the size of epidemic resurgence when it does not go extinct. We compare the full effect of truncation interventions (which influence both the mean of the transmission rate distribution and its shape) to the effects of truncation when is held constant by scaling pSIP, reflecting only truncation effects on the shape of the transmission rate distribution (variance, skew, etc.).

With few infected individuals and , stochasticity and heterogeneity in can either lead to extinction, moderate resurgence, or explosive resurgence. Keeping interventions in place that remove even a tiny percent of the largest values can help to avoid the more explosive events (Fig. 4). Truncation markedly reduces the probability of explosive epidemic resurgence (Fig. 4A) both by increasing the extinction probability (Fig. 4B) and by reducing the magnitude of resurgent epidemics when they do occur (Fig. 4C). While epidemic size was less sensitive to the number initially infected when resurgences do occur (Fig. 4C), the stochastic extinction probability was extremely sensitive to the difference between even one, three, or five remaining infections (Fig. 4B). Much, but not all, of the benefit of truncation comes from changing the mean transmission rate (and therefore ). When is held constant by adjusting pSIP, effects of truncation are more moderate. An increase in efficiency at truncating the top 0.1% of the distribution noticeably decreases the number of infected 42 days after interventions are relaxed (Fig. 4D,F). However, because of the need to slightly reduce pSIP to hold constant under truncation, truncation of marginally decreases the extinction probability of the epidemic, which remains much more sensitive to the number initially infected (Fig. 4E). The highly stochastic nature of epidemic growth when cases are rare, combined with the fact that each truncation leaves behind highly skewed distributions regardless of the truncation parameters, results in even 10,000 epidemic simulations producing noisy patterns across intervention scenarios. Similar patterns are seen as more of the distribution is truncated (Figure S13).

Fig. 4

Effects of transmission rate truncation on epidemic die-out and explosive resurgence. With skewed individual variation in transmission rate, relaxing social distancing interventions when infections become rare (allowing to increase above one) may lead to explosive stochastic epidemic resurgence. Top panels (A–C) show the overall effect of truncation interventions, including effects on both the mean and shape of the transmission rate distribution, and resulting . Bottom panels (D–F) show the effect of truncation when is held constant by rescaling pSIP at the time of intervention relaxation. Specifically, for a 0% truncation efficiency we simulate epidemic resurgence assuming , which results in an at the time of resurgence, which will vary by simulation (where is the number of susceptible individuals and is the total population size). In panels (A–C) as truncation efficiency increases decreases; in panels (D–E) we scale pSIP to retain an average across truncations. Simulations are performed with varying efficiencies of truncation of the top 0.1% of the distribution. Envelopes in (A) and (D) show the central 98% of resurgent simulations (across 10,000 total simulations) for three efficiencies of truncation (0% in orange, 60% in green, 100% in blue). The proportion of epidemic simulations that go extinct within 42 days of intervention relaxation for thresholds of 1 (red), 3 (gold), and 5 (blue) infected individuals is shown in (B) and (E). The upper 99th percentile of concurrent infections 42 days after intervention relaxation in resurgent simulations for the same thresholds is shown in (C) and (F).

Discussion

Understanding local epidemiological dynamics of COVID-19 – and the impact of heterogeneity on those dynamics – remains a challenge due to both limited and imperfect data in most regions and ever evolving interventions and adherence. Reported cases are only a small fraction of all infections, and the proportion of symptomatic cases that are detected remains highly uncertain and variable over space and time. Our approach takes an important step toward capturing locally-specific epidemic dynamics and the impact of heterogeneity across settings by providing a platform (including a mathematical model and open access code) for estimating time-varying transmission rates () from death, mobility, and imperfectly observed case report data, all of which are publicly available. The model can estimate epidemic dynamics and transmission rates over time across epidemiological settings that vary in population size, demography, and control. By incorporating individual variation in contact rates (or, equivalently, infectiousness) into time step transmission rate distributions, we incorporate some of the known effects of heterogeneity without requiring detailed information on population mixing, structure, social networks, or movement patterns. We find that initial shelter-in-place measures in March of 2020 rapidly brought the average reproduction number——from 2–4 to below 1 in all locations we considered in early April. However, from early April through early July, relaxed social distancing led to a steady increase in in all locations, which we estimate to have peaked around 1.2, with some variation among the locations we considered. In response to the second wave brought on by , increased social distancing in all locations starting in early July once again successfully reduced to just below one by mid August.

Non-pharmaceutical interventions have been and will continue to be necessary to control COVID-19 in all settings until better pharmaceutical options (in particular, effective vaccines) are widely available. Social distancing in the general population is effective but costly: it is a blunt and imprecise tool. The social and economic necessity of relaxing social distancing demands safe long-term non-pharmaceutical intervention strategies based on more precise, targeted interventions to reduce transmission. Testing and isolating symptomatic people, combined with contact tracing, remains the gold standard intervention for limiting onward transmission as social distancing is lifted, but it is expensive and capacity remains limited in many settings. Our model shows that it is possible to target interventions even without precise information on specific population mobility, mixing, and infectiousness patterns, by limiting just the most high-risk activities. Specifically, we find that our method of truncating the upper percentiles of individual-level transmission rates, which models a restriction in the ability of individuals to contact a large number of susceptible people in a short period of time (such as by restricting their access to large gatherings and indoor events that have many close contacts), and thereby reduces population-level transmission rates, can be highly effective at maintaining epidemic control (Fig. 2), particularly when combined with mild to moderate social distancing (Fig. 3). Importantly, even after epidemic control is achieved and case numbers drop very low, “chopping off the tail” can provide powerful insurance against explosive resurgence after social distancing interventions are otherwise lifted (Fig. 4).

Five types of factors tend to promote superspreading: (1) high rates or intensity of contact between people or with surfaces; (2) large aggregations of people; (3) poorly ventilated physical environments, especially indoors (Lu et al., 2020); (4) highly infectious individuals; (5) highly susceptible recipient population (Kumar et al., 2020, Althouse et al., 2020, Frieden and Lee, 2020). Many settings where SARS-CoV-2 superspreading has occurred—including nursing homes (Arons et al., 2020), exercise classes, bars and restaurants (Lu et al., 2020), funerals, churches (Shim et al., 2020), meat-packing plants (Dyal, 2020)—combine multiple risk factors. For example, choir practices combine high densities of people, a high-risk activity (singing) (Stadnytskyi et al., 2020), and potentially poorly ventilated indoor spaces; long-term care facilities combine mobile, high-contact caregivers with highly vulnerable residents, often in high-density indoor spaces. In practical terms, “chopping off the tail” as used here aggregates these individual risk factors into the single metric of transmission rate, and models the reduction in transmission rate over four hour periods which we describe as superspreading events—short periods of time characterized by very high rates of transmission.

Some superspreading events may be easier to eliminate than others. Clearly, healthcare and long-term care facilities serve critical functions despite their high-risk nature, and taking all possible steps for decontamination and personal protection in these facilities is critical to mitigate this risk (Frieden and Lee, 2020). On the other end of the spectrum, voluntary, large, indoor events that are mainly for entertainment and could be postponed – gyms, clubs, sporting events, concerts, large lectures – may be the most viable option to reduce superspreading and “chop off the tail” of the contact rate distribution (Benzell et al., 2020). While these common sense interventions are not novel suggestions (Althouse et al., 2020), and are already part of reopening plans in almost all locations, our work allows a direct comparison of how much general social distancing is avoided by eliminating a fraction of these high-risk events (Fig. 3). Truncation strategies are even more desirable in light of their effectiveness at preventing explosive resurgence after controls are otherwise lifted (Fig. 4). Mapping actual event types onto the contact rate distribution to determine how particular superspreading reduction policies would affect control remains an important next step. Importantly, associating superspreading with events and locations, rather than specific people, can avoid the stigma sometimes associated with being identified as a superspreader (Kumar et al., 2020). While superspreading individuals (e.g., those with consistently high viremia, contact rates, or duration of shedding) likely contribute to the spread of COVID-19, our model cannot compare the relative importance of superspreading events and superspreading individuals nor discount the importance of superspreading individuals. We emphasize instead that removing just the highest transmission rates by all individuals can have a large positive impact on the population level transmission rate without the great difficulty and cost of identifying and sequestering specific individuals.

The impact of truncation interventions is two-fold. First, removing the upper tail of the individual transmission rate distribution reduces the population-level mean, often dramatically (Fig. 3A,B). If the mean transmission rate already placed near 1 (for example, due to other interventions), then additional truncation could be enough to cross this critical threshold. However, most intervention strategies that bring to 1 already include prohibiting large gatherings, especially indoors, so additional truncation may not be possible within the context of first-wave interventions. However, truncation also acts on the variance and skew of the transmission rate distribution, though these effects are smaller than the effect on the mean (Fig. 4D–F compared to A–C). Given that super-spreading events are particularly dangerous when cases are few (in the early or late phases of the epidemic) (Lloyd-Smith et al., 2005), sustained truncation interventions could be extremely important for preventing explosive stochastic re-emergence when low case numbers allow general social distancing to be lifted (Fig. 4). In this scenario, resurgence remains rare (Fig. 4B) but possible because individual variation in transmission rates is large; most of the time infectious people transmit to few others, but occasionally someone infects dozens (Fig. 4A), quickly overwhelming testing, contact tracing, and isolation efforts. Sustained truncation dramatically reduces the probability of explosive resurgence, and constrains incipient transmission chains to be smaller and more manageable.

One limitation on understanding the effect of heterogeneity in transmission in particular locations is the challenge of estimating epidemiological parameters from noisy and imperfect data: necessarily a balancing act between model simplicity and complexity. Here, we rely on metrics of heterogeneity previously estimated for SARS-CoV-1 and SARS-CoV-2 (Lloyd-Smith et al., 2005, Endo et al., 2020a, Grantz et al., 2020) instead of estimating them directly from data; we focus our parameter estimation on the mean of the transmission rate distribution. Heterogeneity in contact rates or infectiousness, and the resulting distributional variance and skew, may vary based on local patterns of movement, contact, behavior, and population demography. This heterogeneity can have important consequences: in some cases epidemics with low mean can actually infect a larger proportion of the population than epidemics with higher mean – as was the case for the 1918 influenza pandemic as compared to the 2014 Ebola outbreak – due to the heterogeneity in transmission rates, as described by higher moments of the secondary case distribution (Hébert-Dufresne et al., 2020).

Despite our focus solely on transmission rate heterogeneity, many sources of individual-level heterogeneity are likely to affect transmission, including variation in susceptibility, which lowers the herd immunity threshold and therefore produces a smaller epidemic than would be predicted in a homogeneously susceptible population (Britton et al., 2020, Gomes et al., 2020). Differences in the length of the viral shedding period, where individuals are infectious for a long period of time, could also create superspreading individuals and impact transmission, and possibly impact the efficacy of the truncation intervention we propose here. For example, because we assume that truncation is only partially effective, if that effectiveness is heterogeneous across individual types (i.e., worse at reducing spread by individuals with long-term infectiousness), then this extra layer of heterogeneity could cause us to over-estimate the effectiveness of our modeled interventions. Finally, it is also possible that we overstate the effectiveness of tail truncation because infectious individuals behave differently than how we model them. In particular, we truncate by re-sampling values from the individual-level transmission rate distribution for any values that fall above a given threshold. This re-sampling can be interpreted as individuals substituting the highest transmission rate activities (e.g., a concert) with activities that are represented by the remaining transmission rate distribution (e.g., staying at home, or meeting a couple of friends). If instead individuals substitute this unavailable activity with the highest transmission rate activity that is available, our estimates for the effectiveness of truncation will be overstated.

The true epidemiological parameters in any given location, and the extent of our uncertainty in these parameters, also remain unknown because of the computational challenges of parameter estimation given the limited information contained in noisy case, death, and mobility data. For example, depending on how a particular candidate parameter combination weights the noisiness of cases and deaths and estimates initial conditions, transmission rate estimates can vary substantially (Fig. 1). Fully characterizing uncertainty in model structure and parameter values proved difficult; some pairs of parameters were highly correlated (e.g.,  and : see Figures S2–S6), and the estimated values for many parameters from the best fitting parameter sets were indistinguishable from the full range of parameter values sampled over to fit the model (Figure S10). Future work that directly estimates case ascertainment rates (e.g., through metrics of percentage of tests that are positive, age distributions of positive tests, epidemiological contact information on cases, and analysis of viral genome sequences (Miller et al., 2020), as well as more detailed mobility and contact network information (Ferretti et al., 2020) could help to improve the model fit to the full shape of the transmission rate distribution.

First-wave interventions that eliminated large social gatherings and indoor activities and mandated mask-wearing and physical distancing have likely already affected the heterogeneity in transmission rates, by eliminating many of the high-risk events likely to fall into the upper tail of the distribution. It is important to recognize that as social distancing interventions relax, sustaining such truncation interventions may be critical for keeping transmission down to levels manageable through testing, contact tracing, and isolation. This truncation strategy can potentially reduce the social and economic costs of non-pharmaceutical interventions on the general populace, and facilitate sustained adherence by allowing lower-risk activities to resume while insuring against a resurgence. Ultimately, an unmitigated epidemic, whether as a first, second, or third wave, would kill thousands to tens of thousands of people in each of the locations we studied, reinforcing the point that aiming for population herd immunity through naturally acquired infections is not a viable public health strategy. Instead, exit strategies that can sustain epidemic control after shelter-in-place orders end, including truncating the transmission rate distribution, will be necessary until an effective vaccine can be widely distributed.

CRediT authorship contribution statement

Morgan P. Kain: Conceptualization, Methodology, Software, Investigation, Writing - original draft, Writing - review & editing. Marissa L. Childs: Conceptualization, Methodology, Software, Investigation, Writing - original draft, Writing - review & editing. Alexander D. Becker: Software, Investigation, Writing - review & editing. Erin A. Mordecai: Conceptualization, Methodology, Writing - original draft, Writing - review & editing, Supervision.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.