In December 2019, a novel coronavirus emerged in Wuhan, Hubei Province, China [1]. The pathogen causes a respiratory illness, now known as the coronavirus disease (COVID-19) [2,3]. From its original epicenter in Wuhan, the virus spread rapidly within 30 days to other parts of Mainland China and exported to other countries [4-8]. As of April 10, 2020, 210 countries and territories have reported 1,673,423 confirmed cases of COVID-19 and 101,526 deaths [9]. Due to the spread across multiple countries and the large number of people impacted, on March 11 the World Health Organization recognized the novel severe acute respiratory syndrome (SARS) coronavirus 2 as a pandemic that poses a major global public health threat [10,11].

Although the effects of the COVID-19 pandemic are experienced worldwide, many key health policy decisions designed to reduce transmissions are determined at national and regional levels. These critical policy decisions must be implemented quickly and evaluated continuously so they can be adapted to the local context, recognizing the clear effect that geography, community context, density, and social determinants of health have on COVID-19 outcomes. In North Carolina, the first COVID-19 case was reported on March 2, 2020, and cases increased to 3,963 total confirmed cases as of April 10 [12]. To slow the rapidly increasing transmission rate, within a few weeks after the first case was detected, North Carolina state officials promoted social distancing strategies (ie, deliberately increasing physical space), banned large social gatherings, and closed public schools and universities. Subsequently, a stay-at-home order, which only allows for essential travel outside the home, was issued in the southwestern part of the state by Mecklenburg County effective at 8 am, March 26, lasting through April 16 (since extended to April 29), while a statewide stay-at-home order was issued effective at 5 pm, March 30, lasting to April 29.

Because the COVID-19 landscape evolves rapidly due to the confluence of locally relevant factors, appropriate modeling using timely infection prevalence to drive decision making around containment, treatment, and resource planning is critical. Forecasting models are used to generate early warnings to identify how a pandemic might evolve. During the early stages of the COVID-19 pandemic, forecasting was frequently applied to predict national and international infection transmission trends [13,14]. Local communities and health systems turned to these national and international models for their own planning; however, the generalizability of such models to the local situation is limited and ignores important community-level population characteristics and transmission dynamics [3,15-17]. An objective of this study was to understand how spatial differences impact model results and their interpretation.

In response to the need for actionable data insights in our community and health system, investigators from the Atrium Health Center for Outcomes Research and Evaluation developed a series of COVID-19 forecasting models, which were used to guide Atrium Health’s initial proactive response to ensure sufficient capacity to treat the expected surge in patient care demands. In this study, we present an initial susceptible-infected-removed (SIR) epidemic model and its evolution to the susceptible-infected-removed-social distancing-detection rate (SIR-int) model. In this paper, we describe and compare these models, the spatial differences in a pandemic, the significance of observed cases versus actual prevalence in the setting of rapidly evolving testing strategies, the current epidemiological trends, and the potential effects of nonpharmaceutical interventions applied locally (eg, social distancing).

The observed cumulative case and death counts were obtained daily at noon starting March 2, 2020, when the first COVID-19 case was reported, from the North Carolina Department of Health and Human Services website for all 100 counties [12]. Data collection for this manuscript ended on April 7, just prior to submission. To accurately estimate the actual latent prevalence at time t, the cumulative case counts were adjusted for imperfect detection by dividing them by 0.14. Although estimates of detection probability for the coronavirus, also known as the ascertainment rate, vary in the literature, ours is in line with those reported [18-22]. Modeling only the observed prevalence will give an inaccurate timeline of pandemic behavior. Cumulative deaths were then subtracted from adjusted cumulative cases. We also adjusted cumulative cases for recoveries by removing cases after 20 days, the estimated median duration of viral shedding from illness onset [23]. Daily incremental incidence was obtained by subtracting the estimated latent prevalence at time t – 1 from that at time t. Crucially, in our research, we model estimated latent prevalence as constructed here, not observed prevalence. For the sake of brevity moving forward, we use the terms “latent prevalence” and “prevalence” interchangeably.

In addition to North Carolina, interest also lay in the subpopulation served by Atrium Health’s greater Charlotte market. For convenience, we make use of a designation of a group of counties that constitute the greater Charlotte area used by the North Carolina Department of Health and Human Services. Specifically, the US Centers for Disease Control and Prevention’s Cities Readiness Initiative (CRI) is a federally funded program designed to enhance preparedness in the nation’s largest population centers to rapidly and effectively respond to large public health emergencies such as an act of bioterrorism. This also allowed us to model on a region that harmonized with the state’s approach to disaster planning in case statewide coordination of resources would be required. Within North Carolina, 11 counties are grouped into a CRI region that includes Anson, Cabarrus, Catawba, Cleveland, Gaston, Iredell, Lincoln, Mecklenburg, Rowan, Stanly, and Union. Collectively, we henceforth refer to these counties as “the CRI” (Figure 1). Because the CRI closely mirrors the large area served by Atrium Health’s greater Charlotte market, we used this population base for our local modeling efforts. The CRI includes over 2.5 million residents (24% of the North Carolina population) and ranges from rural settings like Anson County to Mecklenburg County, which contains North Carolina’s largest city, Charlotte [24]. To understand how spatial differences impact model results and their interpretation, we compared the CRI to North Carolina throughout the early phases of this pandemic.

We introduce the SIR deterministic compartmental model originally described by Kermack and McKendrick [25] and depict it in Figure 2.

S is the number of individuals that are susceptible to infection in the population; I is the number of individuals that are infected; R is the number of individuals that are removed from the population via recovery and subsequent immunity or death from infection. This mutually exclusive and exhaustive partition is such that S + I + R = N, where N is the closed population size. We further assume all uninfected individuals are susceptible to infection. The transition flow is described by the arrows in the figure labeled with two rates. The parameter β is the infection rate and can be further decomposed as the product of the probability of transmission per contact and the rate of contact per person per unit time. γ is the removal rate.

More formally, the SIR model is a system of three ordinary differential equations (ODEs) involving two unknown parameters.

Note that all of S, I, and R, and their derivatives are functions of time t, such as S = S(t), although we do not denote this notationally here. By how the model is constructed, the first equation in the system returns a number less than or equal to zero, the second equation returns any real number, and the third equation returns a number greater than or equal to zero.

All data analysis was done using R statistical software, version 3.6.2 (R Foundation for Statistical Computing). As described in Churches [26], we used the ode() default solver from the desolve package to solve the system of ODEs defining the SIR model. Next, we used a quasi-Newton method with constraints to find the optimal values for β and γ on (0, 1) by minimizing the square root of the sum of the squared differences between I, which is our prevalence, and its prediction Î over all time t [27]. To establish initial conditions for model fitting, we estimate the population size of North Carolina and the CRI to be 10,488,084 and 2,544,041, respectively, using information taken from census estimates [24]. After obtaining the estimates and , to help assess model goodness-of-fit, we define the following statistic:

Time is indexed from i 1, ... , n, and n is the number of prevalences in the sample. Note that Ī is the average of the I_i’s.

To compare different scenarios for both North Carolina and the CRI, we define an SIR model (SIR-pre) fit to the data from the time of the outbreak until the time of the March 26, 2020, Mecklenburg County stay-at-home order. Since Mecklenburg County is the state’s second largest county, this could have a strong effect on the pandemic trajectory, both in the CRI and the state; therefore, we have used this date to delineate the date of the significant public health intervention. We further define an SIR model (SIR-post) fit to the data from the time of the outbreak until the end of data collection.

Given the major public health intervention implemented on March 26, 2020, we modified the SIR model for both the CRI and North Carolina to accommodate this (denoted SIR-int). SIR models with interventions can be simulated using the EpiModel package. This package provides tools for building, simulating, and analyzing several classes of models for the population dynamics of infectious disease transmission in epidemics. These include not only deterministic compartmental models, but stochastic individual contact models and network models. We first fitted the SIR model as before to the data up until March 26 and extracted the estimates of β and γ. After March 26, we retained the removal rate but modified the infection rate. We set the preintervention probability of transmission equal to 0.015, which is consistent with other viral infectious diseases like SARS and AIDS [28,29]. We then set the rate of contact so that the probability of transmission multiplied by the rate of contact equaled . To simulate the observed intervention, using the default fourth-order Runge-Kutta Method (RK4) ODE solver, we affected the probability of transmission by iteratively decreasing the hazard ratio of infection, given exposure to the intervention (step size of 0.0001) compared to no exposure, until the fitted infection curve yielded a maximum .

For exploratory data analysis, we generated time plots for prevalence, incidence, and both daily and cumulative deaths. The basic reproduction number R₀ is the average number of secondary cases of disease caused by a single infected individual over his or her infectious period in a population where all individuals are susceptible to infection. To estimate R₀, we compute:

and are estimates taken from the model fit. Since the SIR model is fully parameterized by β and γ, we also obtain predictions and over all time t. The percentage of infected at peak prevalence was computed by dividing the maximum Î by the population size N, while the final percentage of infected was computed as the limit 1 − (∞)/N. To estimate doubling time and compute a 95% confidence interval, we modeled incidence growth by fitting a loglinear model as a function of time t using the incidence package.

Figure 3 shows time plots of prevalence, cumulative deaths, incidence, and daily deaths for North Carolina from the start of the outbreak on March 2 up to and including April 7, 2020. The first death was recorded in North Carolina on March 24.

Figure 4 shows time plots of prevalence, cumulative deaths, incidence, and daily deaths for the CRI from the start of the outbreak on March 11 up to and including April 7, 2020. The first death was recorded in the CRI on March 25.

Notably, the prevalence and cumulative death curves for both figures look exponential. Although both incidence curves are increasing, the incidence curves become volatile after the stay-at-home order went into effect. Prior to March 26, 2020, doubling time was estimated to be 2.56 days in the CRI (95% CI 2.11-3.25) and 2.94 days in North Carolina (95% CI 2.33-4.00). Once data after March 26 were included, the doubling times increased and were estimated to be 4.70 days in the CRI (95% CI 3.77-6.22) and 4.01 days in North Carolina (95% CI 3.43-4.83).

Tables 1 and 2 gives a synopsis of the model fits for each location and model type. The estimated R₀ of 2.36 for the CRI prior to March 26, 2020, is more typical of the range of R₀ values given in the literature for COVID-19, while the value of 1.79 for North Carolina is substantially lower [5,30]. After the intervention, the estimated R₀ values for both locations drop to a similar value, although this result was affected by a reduced model fit. A comparison of the efficacy of intervention, defined as 1 – the hazard ratio of infection, gives 0.25 for North Carolina and 0.43 for the CRI. Using these hazard ratios to compute estimates of R₀ from March 26 onward (_{0, post}), we derive 1.34 and 1.33 for North Carolina and the CRI, respectively. This suggests that the COVID-19 outbreak is rapidly decelerating in North Carolina and the CRI after the aggressive public health intervention.

Figures 5 and 6 show plots of the three fitted models’ infection curves for North Carolina and the CRI, respectively, out to April 7, 2020. The behavior in the two plots is the same. The SIR-post model clearly demonstrates a lack-of-fit to the data. For the SIR-int model, we noted the hinge point induces a change of behavior from March 26 onward. The dotted orange line represents the SIR-pre forecast projections from March 26 onward. They are much larger than the actual data.

Figures 7 and 8 show plots of the three fitted models’ infection curves for North Carolina and the CRI, respectively, projected out to the beginning of August. In both plots, we see the dramatic effect of the public health intervention; that is, the so-called “flattening of the curve.” There are two important differences to note between North Carolina and the CRI region. First, the CRI visibly shows relatively more flattening. This effect can be best observed in Table 3 in the peak infected and final infected columns. Moving from the pre to post to int models within a location, the drop in percentage infected is more pronounced in the CRI. In fact, for the SIR-int model, the percentages are virtually the same for both locations; that is, the CRI has “slowed down” to the state as a whole. Second, the date of peak prevalence was initially 8 days earlier for the CRI compared to North Carolina. However, using the current SIR-int model, although both locations showed their infection curves shifting forward in time, the date of peak prevalence is now 3 days later in the CRI (Table 3). To put this into context, for North Carolina, the time duration from the start of the outbreak to the peak prevalence has gone from 49 days to 70 days (43% increase). However, for the CRI, the time duration from the start of the outbreak to the peak prevalence has gone from 32 days to 64 days (100% increase).

Figures 9 and 10 show plots of the three fitted models’ removal curves for North Carolina and the CRI, respectively, projected out to the beginning of August. These plots support what we have observed so far. With the continued intervention, the removal curves are beginning to collapse, which is a behavior we would expect. For the SIR-int model, both locations show a removal plateau being reached roughly around the beginning of July.