Using influenza surveillance networks to estimate state-specific case detection rates and forecast SARS-CoV-2 spread in the United States

2.1 Influenza like illness surge

We identified excess ILI cases by first subtracting cases due to influenza and then subtracting the seasonal signal of non-influenza ILI (Figure 1). Many states, including Washington, New York, Oregon, Pennsylvania, Maryland, Colorado, New Jersey, and Louisiana, have had a recent surge in number of non-influenza ILI cases far in excess of seasonal norms. For example, in the second week of March, 2020, Oregon saw 50% higher non-influenza ILI than it had ever seen since the inception of the ILINet surveillance system within the US. We find that with 95% probability approximately 4% of all outpatient visits in Oregon during this time were for ILI that could not be explained by either influenza or the normal seasonal variation of respiratory pathogens. We find that as the seasonal surge of endemic non-influenza respiratory pathogens declines, this excess ILI correlates more strongly with state-level patterns of newly confirmed COVID cases suggesting that this surge is a reflection of ILI due to SARS-CoV-2 (Pearson ρ > 0.35 and p < 0.01 for the last three weeks; Figure S1A).

• Download figure
• Open in new tab

Figure S1:

(A) Excess ILI correlates strongly with patterns of newly confirmed COVID cases. This correlation is strongest for the last three weeks of data, when other seasonal respiratory pathogens are at their lowest. (B) Of all ED visits for ILI in New York City, the proportion (black line) of those severe enough to warrant admission to the hospital has increased in the past month.

• Download figure
• Open in new tab

Figure 1:

The excess non-influenza ILI is extracted from all non-influenza ILI by identifying the amount of non-influenza ILI in excess of seasonal norms (blue point and error bars represent the posterior median and 95% credible set for ILI not explained by non-COVID endemic respiratory pathogens). A binomial logistic-normal non-linear regression model was fit to non-influenza ILI data from 2010-2018 (grey lines). The model predicted the expected amount of non-influenza ILI in the 2019-2020 season (grey ribbons represent the 95% and 50% credible sets; the black line represents the posterior median). Observed non-influenza ILI beyond seasonal norms are shown as a blue line (posterior median) and blue ribbon (posterior 95% credible set). A number of regions are not represented due to insufficient laboratory influenza data to complete our analysis (see Methods for full details).

To equate this surge to state-wide or national case counts, we assume that the average number of patients seen per week by sentinel providers is representative of their respective states that week. Using this assumption, the total excess non-influenza ILI across the US was approximately 28.5 million excess individuals in the 3 week period between March 8 and March 28, 2020 compared to the same period in 2019 (95% credible interval of 25.9-30.5 million). Notably, we find that the ILI surge appears to peak during the week starting on March 15 and subsequently decreases in numerous states the following week; the notable exceptions being New York and New Jersey, two of the states that have been the hardest hit by the epidemic.

2.2 Investigating ILI Admission Rates

Our prevalence estimates could be falsely elevated if patient behavior has recently changed leading to increasing detection of mild ILI. If the ILI surge reflected higher rates of detection of typically mild ILI, we would expect emergency department ILI rates would increase yet the proportion of those ILI cases admitted to the hospital would decrease. We were able to obtain data to evaluate this hypothesis from New York City’s influenza surveillance network [9]. In the month of March, the daily number of ILI visits to emergency departments across New York City increased while the proportion of those who went on to be admitted also increased by as much as 3-fold compared to the baseline rate prior to March (Figure S1). This suggests that patients are presenting less often for mild ILI, and such decrease in care-seeking behavior, if similar across the US, could be deflating the size of the ILI surge in later weeks of March.

2.3 Syndromic Case Detection Rate

The rate at which SARS-CoV-2+ patients with ILI symptoms are identified as having COVID varies by state and over time (Figure S2). Our estimated syndromic case detection rates have been increasing over the month of March, which can be expected given increases in testing capacity across the US since the February 28 detection of community transmission in Washington State. For the week ending March 14, COVID cases in the states with the highest estimated syndromic case detection rate (Washington, Nevada, and Michigan) are only capturing approximately 1% of ILI surges in those states (Figure 2). In the last week ending on March 28, we estimate the detection rate across the US increased to be 1.5% (95% credible interval 1.2%-1.9%).

• Download figure
• Open in new tab

Figure S2:

Since March 1, 2020, the case-detection of symptomatic COVID patients has increased by a factor of ≈ 100. This likely represents increased awareness of community transmission within the US combined with increased availability of testing. Still, the syndromic case detection rate remains below 1% for most states with many states with detection rates closer to 0.1%.

• Download figure
• Open in new tab

Figure 2:

Assuming the non-influenza excess ILI for the week starting March 22 consists entirely of patients with COVID, the probability that a symptomatic COVID+ patient will be detected varies by state but even the highest syndromic case detection rates are likely below 10%. Across the US the syndromic case detection rate is 1.5% (95% credible set 1.2% to 1.9%). In Figure S2 we show how the syndromic case detection rate varies over time across states.

2.4 Epidemic Growth Rates and Clinical Rates

The true prevalence of SARS-CoV-2 is unknown. However, if we assume the excess non-influenza ILI is almost entirely due to SARS-CoV-2, an assumption that becomes more valid as the virus becomes more prevalent, we can use the excess non-influenza ILI to define bounds and understand the mutual dependence of exponential growth rates, the rate of subclinical infections, and the time between the onset of infectiousness and a patient reporting as ILI Figure 3. With a January 15 start date of the US epidemic [10], allowing early stochasticity from start-time to the onset of regular exponential growth, we find that it’s impossible to explain the ILI surge with an epidemic whose doubling time is longer than 3.5-days, as such slow growth scenarios fail to produce enough infected individuals to match the observed excess ILI.

• Download figure
• Open in new tab

Figure 3:

(a) The excess ILI estimated falls within the range of what one could expect from a US epidemic growing at the exponential growth rate defined by the growth of new deaths. Adjusting the ILI surge based on decreased care seeking in New York does not reconcile the difference bewteen the ILI surge and the epidemic curve from US growth rates, suggesting additional forces are at play. The cause of the apparent deceleration in the ILI surge is hypothesized to be due to some combination of successful interventions, faster decreases in care-seeking behavior changing than measured in New York, and/or other possibilities including faster growth and higher subclinical rates. (b) If the growth rate in the US is faster than US deaths suggest, such as a growth rate observed in Italy prior to the Italian lockdown (2.645 day doubling time), it could provide alternative explanations of the curvature of the excess ILI through a larger subclinical rate and epidemic curves near their peak at the time of the peak of the ILI surge. Serology or other measures of prevalence are needed to reconcile these alternative hypotheses. (c) More generally, the ILI surge forces a dependence between growth rate (doubling time), the clinical rate, and the lag between the onset of infectiousness and presentation to the doctor with ILI, where faster growth implies a slower clinical rate.

Across the entire US, the doubling rate for deaths due to COVID is 3.01 days (± 0.001, p-value of test that doubling rate is less than 3.5 days approximately 0). Under a 1-day lag from the onset of infectiousness to reporting as ILI, the doubling time of deaths in the US imply an expected 98.6% clinical rate (the proportion of patients who have symptoms for which they would present to a health care provider) if the entirety of the first week of ILI surge is comprised of COVID patients (Figure 3A). Adjusting the ILI surge to account for decreased care-seeking does not produce congruence between the epidemic curve and the ILI surge, suggesting additional factors can be affecting the ILI surge, such as successful interventions, even faster decreases in care-seeking than observed in New York, or epidemic growth rates faster than 3 days.

Faster growth rates, however, require lower clinical rates to explain the ILI surge. If the US epidemic prior to March 14 grew at the rate of deaths in Italy, doubling every 2.65 days, it could better match the curvature of the ILI surge and would imply a clinical rate of 16.5% (Figure 3B). For a four-day lag between the onset of infectiousness and presentation with ILI, the doubling time of US deaths produces, on average, too few COVID cases to explain the excess ILI on March 14. However, 29.8% of the stochastic simulations with a growth rate similar to that of US deaths produced enough COVID cases to explain the ILI surge and thus suggest either secondary introductions, super-spreading, or rapid transmission events in early transmission chains to exceed the ILI surge [11]. On the other hand, the doubling time of deaths in Italy could capture the US excess ILI with a 38.8% clinical rate. If researchers produce estimates of growth rates for the US epidemic, the ILI surge can be used to estimate bounds and ranges of possible clinical rates (Figure 3C). If the entirety of the ILI surge is attributable to COVID, it suggests a slowest-possible doubling time of 3.5 days for the US epidemic starting on January 15.

4.1 Data Sources

Since 2010 the CDC has maintained ILINet for weekly influenza surveillance. Each week approximately 2,600 enrolled providers distributed throughout all 50 states as well as Puerto Rico, the District of Columbia and the US Virgin Islands, report the total number of patient encounters n_it and the total number of which met criteria for influenza-like illness (ILI – defined as a temperature 100F [37.8C] or greater, and a cough or sore-throat without a known cause of than influenza; y_it) [13]. Let d_it denote the number of reporting providers in state i in week t. For scale, in the 2018-2019 season ILINet reported approximately 60 million outpatient visits. Coupled to these data are weekly state-level reports from clinical and public health labs detailing the number of patient samples tested for influenza as well as the number of these samples which are positive for influenza . Therefore ILINet data can be thought of as a weekly state-level time-series representing the superimposed prevalence of various viruses which can cause ILI. ILINet data was obtained through the CDC FluView Interactive portal [14].

In addition to ILINet data, US State population data for the 2020 year was downloaded from https://worldpopulationreview.com/states/. The number of primary care providers in each state per 100,000 residents b was obtained from Becker’s Hospital Review [15]. COVID confirmed case counts were obtained from The New York Times’ database maintained at https://github.com/nytimes/covid-19-data. This dataset contains the daily cumulative confirmed case count for COVID for each state z_il for day l. The dataset of deaths in Italy was downloaeded from https://github.com/pcm-dpc/COVID-19 on April 6, 2020.

4.2 Data Processing

Within the ILINet dataset, New York City and New York were summed into a combined New York variable representing both New York city and the surrounding state. Due to incomplete data in one or more of the data-sources described above the Virgin Islands, Puerto Rico, The Commonwealth of the Northern Mariana Islands, and Florida were excluded from subsequent analysis. In addition, daily cumulative confirmed COVID cases were converted to weekly counts of new cases by

4.3 Extracting non-influenza ILI signal

To subtract influenza signal from y_it we assume that the population of patients with ILI within a state are the same population that are potentially tested for influenza. This assumption allows us to calculate the number of non-influenza ILI cases as

The resulting time-series are shown in Figure S3.

• Download figure
• Open in new tab

Figure S3:

Once the signal attributable to influenza is extracted, the proportion of Patient encounters in which patient had non-influenza ILI displays strong seasonal trends. The most notable deviations from these trends occur around Febuary to March of the 2019-2020 flu season and align with the onset of the COVID epididemic within the US.

4.4 Identifying ILI Surges

We identified ILI surges in by training a model on for all data prior to July 21, 2019. We then used this model to predict the prevalence of non-influenza ILI for dates after and including July 21, 2019. We calculated the ILI surge as the difference between the observed proportion of non-influenza ili and .

More specifically, to account for variation in the number of reporting providers, we trained the following binomial logistic-normal model

We made the following prior specifications: We set the bandwidth parameter for the squared exponential kernel as ρ = 3 representing a strong local correlation in time that died off sharply beyond 3 weeks, α = 1 representing a signal to noise ratio of of approximately 1, v = 1 and ξ = 1 representing weak prior knowledge regarding the overall scale of variation in the the latent space. Finally, we set θ = − 2.197 representing an off-season prevalance of 0.1% non-influenza ILI. Samples from the posterior predictive density were collected using the function basset from the R package stray [16]; a total of 4000 such samples were collected in this analysis. We define the prevalence of non-influenza ILI in excess of normal seasonal variation .

To exclude variation attributable to unseasonably high rates of other ILI causing viruses (such as the outbreak of RSV in Washington state in November-December 2019) we only investigate for weeks after March 7th 2020 as only these later weeks had high correlation to the COVID confirmed case rate (Figure S1).

4.5 Calculating scaling factors to relate ILINet data to COVID case numbers

As COVID new case counts represent the number of confirmed cases in an entire state and ILINet data represents the number of cases seen by a select number of enrolled providers, we must estimate scaling factors w_i to enable comparison of ILINet data to confirmed case counts at the who state level. Let denote the probability that a patient with ILI in state i has COVID as estimated from ILINet data. Let p_i denote the population of state i and let b_i denote the number of primary care providers per 100,000 people in state i. We simulated the number of COVID cases (excess ILI meeting criteria above) as

That is we translate the inferred proportion of individuals with ILI due to COVID to the state level by considering the average number of patients seen by each provider in the study and the number of primary care providers in state . Notably to account for potential errors in these scaling factors, we add propogate uncertainty into our calculation by using Monte-Carlo simulation of the average number of patients seen by each provider in the study.

4.6 Estimating syndromic case detection rates

Assuming that the majority of SARS-CoV-2 testing within the US has been directed by patient symptoms[17], the pool of newly diagnosed SARS-CoV-2+ patients is a subset of the pool of SARS-CoV-2+ patients who are identified as having ILI. Therefore, we calculate the probability that a symptomatic SARS-CoV-2+ patient will be identified as having SARS-CoV-2 as (Figure S2).

4.7 Growth Rate estimation

Growth rates were estimated for the US and Italy by poisson generalized linear models predicting new deaths with date. Data on COVID deaths in the US were obtained from https://github.com/nytimes/covid-19-data on April 6, 2020 and all deaths from March 5, 2020 to April 1, 2020, were summed by date. Initially, April 2-5 were included but were found to have extremely high leverage and were hence excluded from our analysis. Data on COVID deaths in Italy were obtained from https://github.com/pcm-dpc/COVID-19. The same procedure was applied, focusing on deaths from February 24 until March 12. The slope from poisson regression was used as the estimated exponential growth rate, yielding a US growth rate r_US = 0.23 or a 3.01 day doubling time and r_IT = 0.26 or a 2.65 day doubling time.

4.8 Epidemic simulations and clinical rates

SEIR models, were parameterized for the US to a timescale of units days by setting ζ = 3.23× 10⁻⁵ corresponding to a crude birth rate of 11.8 per 1000 per year, a baseline mortality rate ω_b = 2.38× 10⁻⁷ corresponding to 8.685 per 1000 per year, an infectious mortality rate ω_i = 2.62 ×10⁻⁷, incubation period γ⁻¹ of 3 days, infectious period ν⁻¹ of 10 days, and β parameterized to ensure I(t) grew with a specified exponential growth rate early in the epidemic. A total of 2,000 simulations were run for each of the two growth rates (US and Italy) analyzed. Growth rates were drawn at random with r_US N (r_US, 0.1) and r_IT N (r_IT, 0.1). To illustrate the mutual dependence between estimates of growth rate, clinical rate, and the lag between the onset of infectiousness to presentation to a doctor with ILI, 2,000 simulations with uniform growth rates in the interval [0.173,0.365] corresponding to a range of doubling times between 1.9 days and 4 days.

Each simulation was initialized with (S, E, I, R, t) = (3.27 ×10⁸, 0, 1, 0, 0) where time 0 was January 15. To simulate the stochastic time it took from the first case to the onset of regular exponential growth, a Gillespie algorithm was used from the initial conditions until either t = 50 (March 5, 2020) or E(t) + I(t) = 100. The initial Gillespie algorithm was implemented on the assumption that a large amount of variation in the epidemic trajectory stems from uncertainty in trajectory of early transmission chains. The output from Gillespie simulations was input as an initial value into the system of differential equations and integrated until the August 5, 2020. The number of infected individuals on a given day was the last observed I(t) for that day, and a weekly pool of infected patients was computed by a moving sum over the number of infected individuals every day for the past week,.

Defining as the national excess ILI, the clinical rate implied by a given simulation was estimated as for a given time delay t_d it takes from the onset of infectiousness to a patient reporting to the doctor with ILI.

4.9 Code Availability

All code and data required to reproduce our results is publicly available at https://github.com/jsilve24/ili_surge.