Determining the best mathematical model for implementation of non-pharmaceutical interventions

1.   Introduction

SARS-CoV-2 is a virus that affects the respiratory tract and can lead to death or hospitalization ^[1]. Since its initial discovery in Wuhan, China, SARS-CoV-2 has come to spread across all corners of the world ^[2,3]. Before the development of effective vaccines, the primary defenses against the virus were behavioral interventions such as masking, lockdown, and quarantine measures ^[4].

When a novel virus emerges, preventative measures outside of antiviral treatment or vaccines are the first line of defense, so it is essential that we understand the aggregate effect of such practices on infectious diseases ^[5]. The spread of disease is dependent on a variety of factors and, likewise, our methods of decreasing the transmission of infectious disease are manifold. Battling the early stages of future pandemics and epidemics will stem from our ability to implement and enforce effective preventative measures and non-pharmaceutical interventions (NPIs) prior to vaccines and effective treatments. In order to plan for any future spread of disease we must understand how effective we are at curbing the transmission of disease so that we can develop sufficient public health infrastructure, responsible practices as citizens, and overall plans to fight infectious disease ^[6]. We need to understand the mechanism by which these non-pharmaceutical preventative measures affect infection rates for future outbreaks of infectious disease to better predict and plan for potential cases ^[7]. By doing this, it will allow us to ensure our public health infrastructure is prepared and can adequately deal with the influx of cases in future pandemics, which will be essential in managing the spread of a tarrataca21future pandemic and for ensuring general public safety ^[8].

Across the US, many of the SARS-CoV-2 responses were local or state-based rather than a singular federal response, and this has led to varying results and effectiveness ^[4,9]. Early in the pandemic, many states implemented lockdown measures, shutting down both private and public institutions and asking people to stay home if they could ^[10]. Less stringent measures were introduced as well, like mask mandates and social distancing policies in public places. These preventative measures typically revolved around changing the way in which people behaved. Generally, this caused transmission rates to drop ^[11,12]. Additionally, as the specifics of lockdowns, like timings and measures enforced varied by region, effectiveness also varied by region ^[13]. This state-to-state variation provides an opportunity to examine how NPIs might best be incorporated into epidemiological models.

Compartmental models for the spread of a disease through a population are a well-established tool for analyzing outbreaks of infectious disease ^[14]. The basis of most modern compartmental models was first proposed in the 20$^\mathrm{th}$ century by Kermack and McKendrick ^[15]. With respect to SARS-CoV-2 and its global pandemic, compartmental models have been an excellent tool for the analysis of the spread of SARS-CoV-2 at all scales and for forecasting how certain measures like NPIs might affect the spread of SARS-CoV-2 ^[14,16]. In particular, compartmental models have been used to investigate the effect of lockdowns ^{[17,18,19,20]}, masking ^[21,22], social distancing ^{[23,24,25,26,27]}, or non-specific NPIs ^[28,29,30]. The majority of these studies model NPIs as changing the effective transmission rate, with many assuming an abrupt change in this parameter when the intervention is applied ^{[31,32,33,34]}, examining the effect of NPIs in places as varied as Brazil ^[35], Germany ^[36], China ^[37], Korea ^[38], the UK ^[39] and Africa ^[40]. However, human behavior is rarely so abrupt, and implementation of policies throughout a region or country can take some period of time. Some modelers have attempted to account for this by incorporating a smoothing function to transition from one transmission rate to the next ^{[41,42,43,44,45,46]}, although it is not clear what type of function should be used to model the change in the transmission rate.

In our study, the goal was to test several transition model formulations to establish which one best represents reality. The primary goal of this study was not to identify the exact effect by which each NPI affected the infection rate and total progression of SARS-CoV-2 through the region, but rather the mode by which the transmission rate changed after the implementation of NPIs. To do this, we used the cumulative case data for US states to fit the different models of decay in the transmission rate. We found that the exponential model is the best model for the most states, followed by the logistic model. The instantaneous and linear transition models are rarely the best transition model.

2.   Materials and methods

2.1.   Epidemiological model

To model SARS-CoV-2, we use a standard SEIR (Susceptible Exposed Infected Recovered) model as the basis for our work. A SEIR model was used in particular because of the known latency period between exposure and infection in cases of SARS-CoV-2 ^[47]. The model is composed of four differential equations representing individuals who are susceptible ($S$), exposed ($E$), infected ($I$), and recovered ($R$). The model has three parameters that describe spread of the infection: $k$ is the rate at which individuals move from exposed to infected compartments and is given by the inverse of the incubation period; $\delta$ is the rate at which individuals move from infected to recovered and is given by the inverse of the recovery time; and $\beta$ is the transmission rate that moves individuals from susceptible to exposed individuals. $N$ is the total population for each region and is assumed to be constant. The model is given by the equations

All parameters in the SEIR model, as well as the parameters in our models of the decay in $\beta$, are described in Table 1. A flow chart of the SEIR model is given in Figure 1.

Our implementation of the SEIR model assumes a constant population of $N$ inhabiting a homogeneous area, and that individuals flow from one compartment to another, only flowing in one direction. Individuals start in the susceptible compartment and then flow into the exposed compartment. The rate at which the number of susceptible individuals is decreasing at any given point in time is determined by the product of susceptible individuals, infectious individuals, and $\beta$. These people then enter the exposed compartment, where they remain for an average time $1/k$ before moving into the infected compartment. They remain in the infected compartment for an average time $1/\delta$ before finally flowing into the recovered compartment. The SEIR model ends in the recovered compartment, and once an individual has entered this compartment, they will remain there. This means that the SEIR model assumes no possibility of reinfection. The cumulative flow of individuals between compartments can do well to represent the spread of a disease like SARS-CoV-2 through a population.

2.2.   Modeling non-pharmaceutical interventions

We assume that NPIs change the value of the transmission rate ^[11]. That is, we assume a transmission rate $\beta_1$ before the NPIs and a transmission rate $\beta_2$ after the NPIs take effect. It is important to note that we measure the cumulative effect of multiple NPIs. We can not differentiate the effects of each individual NPI on the transmission rate. We test four different mathematical expressions to model the transition from $\beta_1$ to $\beta_2$.

● Instantaneous change in $\boldsymbol\beta$: The most common assumption is that there is an immediate change in the infection rate at the time when NPIs are imposed.

● Linear decrease in $\boldsymbol\beta$: This model allows for a more gradual change in the transmission rate but includes an additional parameter, $t_{2}$, the time at which $\beta$ stops decreasing. We assume that $\beta$ starts decreasing at $t_{ld}$. The linear transition is modeled by the equation

● Exponential decay in $\boldsymbol\beta$: This model assumes an exponential decay to $\beta_2$ beginning after $t_{ld}$. There is not a $t_{2}$ parameter in this model; instead, the parameter $\tau$ is a constant that describes how quickly the change between different $\beta$s occurs. The exponential transition is modeled by the equation

● Logistic change in $\boldsymbol\beta$: This model allows for a possible delay in the onset of the effect of the NPIs as well as a gradual decay to the lower transmission rate. This model has both the constant $\tau$, which describes how fast $\beta$ decays, and $t_2$, a constant that represents the midpoint time over which the change occurs. The logistic transition is modeled by the equation

Mathematical analysis of the SEIR model, even with a time-dependent (or non-autonomous) transmission rate, is covered in other manuscripts ^{[50,51,52,53,54,55]}.

2.3.   Fitting the models to the data

To fit our models, we used cumulative case data from all 50 states in the United States plus the Washington D.C. area. This gave us data from a total of 51 regions that we could fit with the different models. These 51 regions represent a broad range of responses to the pandemic, which might allow us to make general conclusions about how the SARS-CoV-2 transmission rate was affected by the widespread use of NPIs. We used the SARS-CoV-2 case data provided by state governments through their websites and online databases. The calendar start date of each data set is given in a table in the supplementary material.

We minimized the sum of squared residuals (SSR) to fit the SEIR model incorporating the various possible transition models to the cumulative case count of each state, which is given by

Note that the assumption here is that all infected people are identified and accounted for in the cumulative case count, so this really represents the total number of people who have transitioned to the infected state at that point in time. Although lockdowns often had set dates by the government, they might not have necessarily followed those dates in actual effect, so we allowed for this in our model by allowing the time of lockdown, $t_{ld}$, to be a parameter we were fitting. This fitted $t_{ld}$ value represents the effective date of lockdown or NPI implementation in each state rather than the official date announced by state governments. Ultimately, this choice was made to try and best represent the complex nature of how individuals and governments responded to SARS-CoV-2 and account for additional factors that might have facilitated or inhibited the implementation of preventative measures. By fitting this additional $t_{ld}$ parameter, we can also calculate a difference between the effective lockdown date and the declared lockdown date. The declared lockdown date in the scope of this paper refers to when each region's respective governing agency declared a state of emergency for the SARS-CoV-2 pandemic. These dates were sourced from government websites and documents. This difference between official and effective $t_{ld}$ dates gives a measure of how long it may have taken NPIs to go into effect after they were implemented.

To accommodate the relatively well-known values for the incubation period and the recovery time, we fixed the parameters of $k$ and $\delta$ as constant values in our simulations. We assumed that the death rate ($\delta$) and the incubation period ($1/k$) are the least likely to vary between states, since these are largely (though not entirely) determined by virus host interactions. The transmission rate ($\beta$), however, is much more dependent on environmental and social factors, and is more likely to vary from one state to another. Thus when deciding which parameters to fit, we fitted the transmission rate and fixed incubation period and death rate, assuming that state-to-state variation in the latter two was much smaller than state-to-state variation in the transmission rate. The values are given in Table 1. In addition to $k$ and $\delta$, we used one other parameter that we did not fit for each state, namely $N$, the population of the state, which is given in the supplementary material. Using the data gathered from each state, we set each individual $t_0$ as the day case(s) were first recorded in the given state. In all, this left us with the free parameters $\beta_{1}$, $\beta_{2}$, $t_{ld}$, $t_{2}$, and $\tau$ for all the models.

When we fitted each model to find parameters that yielded the smallest sum of squared residuals (SSR), no bounds were set for the parameter values. After finding the best-fit values for each model for each state, we determined 95% confidence intervals for those parameter values using bootstrapping ^[56]. For each model, we ran 1000 bootstrap iterations to determine our 95% confidence intervals. For the bootstrapping, we made use of general parameter bounds to ensure that we achieved plausible values. Those bounds are listed in Table 2. Bootstrapping results were used to make corner plots of the parameter distributions, which are available in the supplementary material.

We used the Bayesian information criterion (BIC) to determine the best-fit model for each state. The BIC attempts to quantify the amount of information gained by adding extra parameters and applies a penalty if the addition of extra parameters has not improved the fit commensurate with the additional available information. The BIC provides a metric for the comparison of models with different numbers of parameters when fitted to the same data ^[57]. Models with lower BIC values are considered to be better representations of the system. The BIC is given by

In the equation, SSR is the sum of square residuals for our best-fit parameter values, $n$ is the size of the data set, and $k$ is the number of parameters in the model. Computing these values allows us to make a value judgment on the relative effectiveness of each of these models in representing the change in $\beta$.

3.   Results

3.1.   Instantaneous change in the transmission rate

We first looked at the assumption of an instantaneous change in the transmission rate. The instantaneous model fit to the data, along with a corner plot showing parameter correlations and distributions for the state of New York, is shown in Figure 2. We use New York as an example since we cannot include graphs for all 51 regions in the main manuscript, but it is not necessarily representative of the fits achieved for all 51 regions. Model fits, corner plots, and 95% confidence intervals for all 51 regions are included in the supplementary material. Histograms of the estimated parameter values for $\beta_{1}$, $\beta_{2}$, and $t_{ld}$ for all 51 regions are also available in the supplementary material. The corner plots summarize all the possible parameter combinations found through bootstrapping. The graphs along the diagonal show the distributions of the parameter indicated on the left of the row or bottom of the column. The remaining plots are correlation plots, where we have plotted combinations of parameters pairwise, with the combination indicated on the left of the row and bottom of the column. Ideally, if parameters are uncorrelated, these plots should be circular. Correlation plots that are linear or show some other curve indicate a correlation between those two parameters, i.e., when one parameter changes, the other parameter adjusts in a known way to maintain the original curve.

The fits for the instantaneous model were decent, but tended to fit poorly during the period of transition. As in the corner plots of Figure 2, many instantaneous model fits suggested a correlation between $\beta_1$ and $\beta_2$. Interestingly, the two values of $\beta$ are negatively correlated, so that an increase in one is linked to a decrease in the other. This means that the change in the transmission rate is not what determines how well the model fits the transition.

3.2.   Linear change in the transmission rate

The second model we used to model NPIs was a linear change in $\beta$ from a higher $\beta_{1}$ to a lower $\beta_{2}$ starting on the day $t_{ld}$ and ending on the day $t_{2}$. The linear model fit to the data from New York state and a corner plot for the parameters in the linear model fit for New York state are shown in Figure 3. Best fits, corner plots, and 95% confidence intervals for all 51 regions are included in the supplementary material. Histograms of the estimated parameters across all 51 regions are also available in the supplementary material.

For the linear model, we again see reasonable fits over most of the time course, except near the time of lockdown, when the transmission rate is changing. We are easily able to calculate the duration over which $\beta$ was decreasing. Across all 51 regions, the median duration of the change in $\beta$ was 24 days, and the mean duration was 27.4 days. The histogram of the duration across all states can be found in the supplementary material. Generally, this shows that in the linear model, it took roughly 3–4 weeks for $\beta$ to change. We again see a negative correlation between $\beta_1$ and $\beta_2$.

3.3.   Exponential change in the transmission rate

The third model we used to describe the change in the transmission rate caused by NPIs was an exponential model, where the transmission rate changed from $\beta_{1}$ to $\beta_{2}$, starting on the day of lockdown and decaying with a time constant of $\tau$. The exponential model fit to the data from New York state and a corner plot for parameters in the exponential model fit for New York state are shown in Figure 4. Best fits, corner plots, and 95% confidence intervals for all 51 regions are included in the supplementary material. Histograms of the estimated parameters across all 51 regions are also available in the supplementary material.

This model generally captures the change in the transmission rate better than the previous two models, with a better fit of the data around the bend in the curve. In the fits for this model, there is no correlation between $\beta_1$ and $\beta_2$, although both $\beta_1$ and $\beta_2$ are negatively correlated with $\tau$. Thus, as the time constant for the transition increases, both $\beta_1$ and $\beta_2$ decrease.

3.4.   Logistic change in the transmission rate

The final model used to fit the data was a logistic change in $\beta$ between $t_{ld}$ and the end of the fitting. This logistic change in $\beta$ was centered on the day $t_{2}$, and the coefficient $\tau$ describes the rate at which the decay takes place in the logistic curve. The logistic model fit to the data from New York state and a corner plot for the parameters in the logistic model fit for New York state are shown in Figure 5. Best fits, corner plots, and 95% confidence intervals for all 51 regions are included in the supplementary material. Histograms of the estimated parameters across all 51 regions are also available in the supplementary material.

One thing to note about the logistic model is that in certain circumstances, it can look very similar to the exponential model. Since it has the most free parameters of any of our models, it also generally fits the data best, since it has the ability to fit the data through the bend in the curve. This model also shows no correlation between $\beta_1$ and $\beta_2$. In fact, $\beta_2$ appears to be uncorrelated with other parameters. The correlations for this model are between $\beta_1$ and $\tau$, and between $\tau$ and $t_2$.

3.5.   Comparison of model fits

To quantify the effectiveness of each model in comparison with each other, we used the Bayesian information criterion or BIC, given by Eq (2.7). The model that best represents the data is the one with the lowest BIC. The BIC for each state and each model is given in the supplementary material. Using each state's fits, we calculated the BIC for each model and then deduced the best model for each state. We then tabulated the best model for each state, which is seen in Figure 6. The exponential model has the highest number of best fits according to the BIC calculations. The logistic model was next, with three-quarters the number of best BIC values. Significantly less than this is the linear model, with only five of the best BIC values. The instantaneous model was the worst of the four models, having four states where it achieved the best BIC value. The exponential and logistic models were both significantly better than the other two models.

We can compare the graphs of how $\beta$ changes in the different models. For each region, we can plot all the different models on one graph to compare them (Figure 7). For all states, these combined $\beta$ plots can be found in the supplementary material. We use the states of Iowa and Nebraska in Figure 7 as examples. For most states, the different transition models start at different $\beta_1$ values but all tend to converge to a specific $\beta_2$ value. These graphs also clearly show the stark difference between using an instantaneous transition and the other more gradual transitions.

3.6.   Parameter comparisons between regions

In addition to comparing the values of $\beta$ between the models, we can compare the values of $\beta$ between states. We can make a histogram of the fitted values of $\beta_{1}$, $\beta_1-\beta_2$, and $\beta_2$ for each model. These histograms can be seen in Figure 8 with $\beta_1$ in the left column, $\beta_1-\beta_2$ in the center column, and $\beta_2$ in the right column. Each row is a different transition model.

As noted earlier, $\beta_1$ tended to vary more than $\beta_2$. This can be seen in the histograms in the range covered by $\beta_1$, which ranges across about 0.5–3.5 /d, while $\beta_2$ ranges across about 0–0.2 /d. The instantaneous model tends to have the lowest values of $\beta_2$, while the linear model tends to have the highest $\beta_1$ values. We can also consider the difference in $\beta$ across the different regions and states. The difference between $\beta_{1}$ and $\beta_{2}$ gives us a measure of the maximum effectiveness of NPIs in reducing the transmission rate of SARS-CoV-2 during the pandemic. The general conclusion made from the histograms is that across all models, there was roughly a 90% decrease from one $\beta$ to the other. This suggests that the cumulative effect of NPIs was significant in decreasing the SARS-CoV-2 infection rate.

We can also examine the differences between the declared date of lockdown and the effective date of lockdown ($t_{ld}$) calculated by our fits. The official date of lockdown is defined to be when states declared an emergency in response to SARS-CoV-2. In the histograms of Figure 9, we can see the distribution of the differences between the models. For most states, the declared date of lockdown came before the estimated date of lockdown. The instantaneous model was the only model in which all fits were after the declared date of lockdown. All other models had a mean difference of 1 to 4 days. Thus, the declared dates of lockdown were close to the effective dates of lockdown, but were generally slightly before the effective lockdown. This suggests that in many states, it took some time for the NPIs to actually take effect. Interestingly, there are some cases where the estimated time of lockdown is before the declared state of emergency, which could suggest that individuals were already beginning to self-isolate before governments officially acted. One important thing to consider about the effective lockdown concerning every model except the instant one is that this $t_{ld}$ marks the beginning of when $\beta_{1}$ goes to $\beta_{2}$, so it may not begin to decrease rapidly at $t_{ld}$ and may stay relatively near $\beta_{1}$ for some time.

4.   Discussion

Using data gathered in the first 180 days of the COVID-19 pandemic in the 50 states as well as the Washington D.C. area, we fitted cumulative case data to several different models of the transmission rate's transition to demonstrate how the transmission rate of SARS-CoV-2 was cumulatively affected by preventative measures like government mask mandates, quarantines, and lockdowns. We found that the exponential model performed the best using our data. The exponential model had 24 states for which it had the best BIC value, which represented a little under half of all the datasets. The logistic model performed the second best among the models used, having 18 states for which it produced the lowest BIC. The linear and instantaneous models performed poorly compared to these two, having five and four best BIC values, respectively. We can conclude that the instantaneous and linear models are generally not good models for reproducing changes in transmission due to implementation of NPIs, with the exponential and logistic models being better able to capture the epidemic's dynamics.

One notable observation was that all models generally approached similar values of $\beta_{2}$, while the estimates of $\beta_1$ varied over a wider range. During the early stages of the pandemic, appropriate testing was not available in the United States, leading to serious under-counting of the number of cases ^[58,59]. Since testing was initially limited to people with symptoms, asymptomatic cases were almost certainly missed ^[60]. This could have contributed to the variability of the estimates of the initial transmission rate. The consistancy of $\beta_2$ estimates points to the broad effectiveness of NPI measures, since they reduced the transmission rate to a similar value across a variety of geographic, population, and socioeconomic conditions.

We also found that the estimated date of lockdown was generally after the declared date of lockdown. This difference was greatest in the instantaneous model. This could be particularly problematic from a modeling standpoint if an instantaneous transition model is used, particularly when it is assumed that the transition happens on the day of a government declaration of NPI measures. Our results suggest that this assumption does not replicate actual case counts during the transition very well. It is important to consider, however, that the difference between estimated and official lockdown could also depend on the exact measures states took when they declared an emergency, as these measures varied from state to state ^[59,61]. There was also some spatial heterogeneity in these measures, as some local governments introduced NPI measures before a state-wide emergency was declared. Despite this, we were still able to see that the effective date of lockdown generally lagged behind the declared date.

Other studies corroborate the idea that there is a delay between the announcement of an NPI policy and the appearance of its effect in the case count data; one study used machine learning to calculate transmission rates ^[62], and one study used deterministic models to calculate transmission rates ^[63], both noting a delay between the application and effect of NPIs. The delayed time of lockdown could be due to a delay in the actual implementation of the government's policies ^[64], but there is also a natural lag in reporting cases ^[65], which could contribute to the later values of the estimated time of lockdown. Another possible contribution to the delayed time of lockdown estimate is that the effect of NPIs is not static. That is, there is not really a single final value of $\beta$, since people's adherence to policies is likely to fluctuate with local disease conditions and will also change for events such as holidays ^[66,67]. We are also limited by the assumption of a single inflection point when, in fact, there are likely to be several points in time when the different measures really took hold ^{[32,42,68,69]}. This could be mitigated by finding the inflection points ^[70] from the data and using appropriate transition functions at each time point.

One reason we likely could not decisively find the best model is due to the spatial heterogeneity in many of the factors contributing to the transmission rate ^[71] and NPI adherence ^[72]. Our model also does not distinguish between different preventive measures, so it does not allow us to determine whether specific NPI interventions are associated with a particular transition function. Further studies could potentially decrease the scale at which the modeling was done or increase the number of areas modeled to better determine trends in the transmission rate. There is evidence suggesting that local lockdowns were preferable and more effective than a global lockdown in controlling the later stages of the pandemic ^[73], so there is an incentive to study local transmission dynamics in smaller regions.

There were several other limitations to our methods. The primary limitation is in the use of cumulative case counts. As mentioned earlier, problems with testing led to under-reporting of cases early in the pandemic ^[58], with asymptomatic cases largely being missed ^[60]. A possible improvement would be to use hospitalizations or deaths for parameter estimation ^[74], both of which tend to be more accurately captured by the available data. This would require the addition of extra compartments to our model. Our model only employed four compartments, but several other compartments could be incorporated to better represent the reality of the spread of SARS-CoV-2 throughout the US. Unfortunately, this would also require the addition of new parameters that would need to be estimated from the data. Some compartments that have been examined in other models include asymptomatic cases ^{[34,60,75,76]}, quarantined patients ^[34,76,77], and hospitalized patients ^[34,78]. Our model also did not include the chance of reinfection of recovered individuals, which could have contributed to the infection rate ^[79,80]. We also did not consider movement of people between states, which would also affect our estimates of the transmission rate ^[81].

Another major constraint of this study is the areas fitted. Large state-sized entities are relatively diverse and are typically spatially heterogeneous, which may have led to some inaccuracies in our fits. Other factors that affect particularly local transmission, such as socioeconomic factors ^[82], population density ^[83], or access to healthcare ^[84] were not explicitly considered in the model. We also assumed that each state's population was constant and insular from other states, which does not necessarily reflect the reality of the SARS-CoV-2 pandemic, where individuals had the mobility to move between areas ^[85]. The consideration of more states or data from smaller, more homogeneous regions, such as cities or counties, might have revealed clearer trends in the results.

5.   Conclusions

Despite these limitations, our study indicates that an exponential change in transmission rate is the best assumption of the four models examined here for modeling the implementation of NPIs. Other transition models should continue to be explored over a larger diversity of areas to further extend knowledge on NPIs' effects on transmission rates during the spread of infectious disease. Future work could also include extending the basic SEIR model used here to include more compartments, such as asymptomatic and hospitalized patients. While we focused on the use of NPIs before vaccines were introduced during the COVID-19 pandemic, the interplay of NPIs and vaccination is also of interest, particularly if the vaccine is only moderately effective. Nonetheless, our study provides a foundation for future work that incorporates realistic changes in the transmission rate as NPIs are introduced.

Use of AI tools declaration

The authors declare that they have not used artificial intelligence (AI) tools in the creation of this article.

Conflict of interest

Hana M. Dobrovolny is a guest editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.