Comparison of effectiveness of enhanced infection countermeasures in different scenarios, using a dynamic-spread-function model

1.   Introduction

In designing intervention strategies to combat future waves of an infection, such as those occurring in the COVID-19 pandemic, it is advantageous to incorporate known behavioral tendencies of the affected population into the analysis. Relevant behaviors include willingness to deploy a certain type of personal protective equipment (PPE), and level of compliance with social-distancing (quarantines, isolations, lockdowns…) mandates. Of particular value are tools that can translate a given level of change in behavior to changes in infection rate. Such information can enable policy makers to estimate the public-health benefit obtainable from resources spent on PPE acquisition and stockpiling, public education, and enforcement of legal orders.

SIR (Susceptible-Infected-Removed) type models represent a family of possible tools for assessing the influence of behavioral characteristics on infection dynamics. SIR models have been used throughout the COVID-19 crisis for predicting various effects of intervention strategies ^{[1,2,3,4,5,6,7,8]}. To incorporate behavioral characteristics into the model, one could adjust the model parameters to best reproduce infection curves that are already available for the scenario of interest, e.g. from a first wave of the infection. One challenge to such an approach is that data available to inform the model is often limited, and parameter selections are sometimes "educated guesses" or values from other infection scenarios (e.g. different type of pathogen) ^[3]. Another challenge to such an approach is that different sets of parameter values can yield the same quality of fit with the recorded data (infection rates, death rates, hospitalization rates…). When the chosen parameter values are varied to model behavioral changes, the outcome can be dependent upon the particular set ^[2]. That is, the sensitivity of the model to changes in the parameters varies for different parameter selections.

An alternative approach that reduces these difficulties was recently introduced by Osborn et al. ^[8]. In that approach, the population behavior as a function of time is governed by a differential equation, along with the standard S, I, and R dependent variables. The "dynamic spread function", which contains the transmission rate (influenced by social distancing) and production rate (influenced by PPE's adopted) is determined from published infection curves for previous waves of the infection. In the governing equation for the dynamic spread function, the published curves (for the infected population, and the number of new infections) enter as variable coefficients. To study the effect of behavioral changes away from the baseline behaviors exhibited in the previous wave, the spread function can be modified in a systematic way, and the system of equations can be solved to provide estimates of the change in infection rate (over the time course of the previous waves) realizable from the proposed interventions. An advantage of the dynamic-spread approach is that parameters that are common to both the baseline and modified scenarios do not need to be specified. The dynamic-spread approach also captures in a natural manner the gradual way in which behavioral changes (e.g. gradual adoption of facemasks) by a population occur.

The following sections summarize the dynamic-spread formulation derived in ^[8]. Calculations are then made to illustrate the significance of the dynamic spread function. Subsequently, four COVID-19 scenarios are considered, and common modifications to the baseline behaviors are imposed. The resulting change in the COVID-19 infection rate is then estimated. The scenarios considered are the COVID-19 outbreaks in New York City, New York State, the Diamond Princess cruise ship, and Harris County, Texas. The effectiveness of different strategies for the different scenarios is discussed in the final section.

2.   Methods

2.1.   Basic methodology

In ^[8], the dynamic-spread methodology is derived from a 4-equation SIR model (^[9,10,11]), which tracks the susceptible, infected, and removed populations, along with the droplet transmission. In the formulation, the transmission rate $\widetilde {\beta }$ and the production rate $\kappa$ are allowed to vary with time. The resulting system of equations can be written as ^[8]:

In the equations, all the variables are made dimensionless based on a scheme discussed in ^[8]. Time is normalized by the temporal window of interest Δ, which we often take to be the time from the first lockdown order to the first maximum in the recorded number of new infections per day. U is the number of new infections per day, normalized by α. We note that Osborn et al. ^[8] derived the governing equations in terms of T = dS/dt = -U, where S is the susceptible population. For convenience, we work here in terms of U rather than T, since T is always negative, i.e. the negative of the number of new infections per day. I is the infected population, scaled by (α Δ). λ = Δ/ν is the dimensionless droplet removal rate, and ${\gamma = {\mathit{\Delta}} \mu }_{I}$ is the dimensionless infection recovery rate, ${\mu }_{I}$ being the dimensional infection recovery rate. The function δ(t) is the dynamic spread function, defined by

The functions I_b(t) and U_b(t) are the published infected population and number of new infections for the baseline scenario (denoted by subscript "b"). If the baseline spread function δ_b(t) is used as the spread function in the governing equations (3a, b), the published baseline profiles U_b(t) and I_b(t) are recovered (within numerical tolerances). The utility of δ(t) derives from modifying it to model alternative intervention strategies and solving Eqs. (3a, b) to determine the new infection curves (i.e., U(t) and I(t) profiles). Modifications to account for new strategies, derived in ^[8], are presented below.

2.2.   Modeling alternative intervention strategies

Three types of enhanced intervention strategies were investigated.

2.2.1.   Enhanced outward filtration efficiency

We assume that the infected population initially deploys a barrier of outward filtration efficiency, FE_{out, b}, where the subscript "b" denotes baseline. If a different mask design characterized by an enhanced efficiency FE_{out, mod} is then deployed by the infected population, then the modified spready function becomes ^[8]

Here ϵ_κ is the fraction of the change in δ(t) due to changes in droplet production ^[8].

2.2.2.   Increased compliance

Eq. (2a) assumes the population compliance is unchanged when the new mask is adopted. To study a scenario where the mask is the same but the compliance rate changes, the following modified spread function applies ^[8]:

Here F is the factor by which the baseline compliance changes. F can be greater than 1 or smaller than 1. The compliance rate itself is prescribed by ^[8]:

2.2.3.   Higher levels of social distancing

To analyze changes in social distancing, say from a baseline level L_b to a modified level L_mod, the modified spread function is ^[8]

In practice, only the ratio L_mod/L_b is important.

Using one of the modified spread functions (Eqs. (2a, b, d)) in Eqs. (1a, b) enables estimation of the change in infection rate for the modified scenario of interest. Detailed examples, using different scenarios, are provided in Section 3.

2.3.   Solution procedure

To solve Eqs. (1a, b, c), values of the droplet removal rate λ, the infection recovery rate γ, and initial value of the spread function δ₀ are required. Information to assist in prescribing these parameters can be obtained from known features of the infection dynamics prior to the first intervention. Prior to the first intervention, the standard SIR formulation given by Eqs. 1a, b with δ(t) set to zero, is applicable. Assuming exponential growth (with exponent M) in the number of new infections per day, it was shown in ^[8] that

The growth rate M can be obtained from infection rates published during the beginning of the epidemic (prior to any interventions). An exponentially growing solution will occur when the reproduction number R₀, given by ^[10]

is larger than 1. Estimates of R₀ for the early stages of epidemics are also published. Details of the initial R₀ estimates (including published literature sources) for the four scenarios are provided in the description of the uncertainty analysis in the Appendix section. Relations (3a, b), along with prescription of the recovery rate ${\mu }_{I}$ (with corresponding dimensionless recovery rates $\gamma$) enable the simulations to proceed. For the four scenarios considered in this study, a common range of recovery times 1/${\mu }_{I}$ between 2 and 10 days (or recovery rates ${\mu }_{I}$ between 0.1 and 0.5 per day) was assumed based on reported observations from the COVID-19 pandemic ^[1,5,12]. This range resulted in a mean recovery rate ${\mu }_{I}$ of 0.3 per day with a standard deviation of 0.1 per day (assuming a 95% confidence interval (CI) for the range of recovery rates). The range of recovery rates along with the mean and standard deviation are considered in the uncertainty analysis (outlined in the Appendix).

The governing equations (1a, b) containing the modified spread function were solved using a Runge-Kutta method (Matlab ode45, Mathworks Inc.).

2.4.   Scenarios considered

We considered four COVID-19 infection scenarios – New York City, New York State, Diamond Princess cruise ship, and Harris County, Texas – to investigate how differences in the behavior of the population affected the manner in which different intervention strategies influenced the course of the infection. As noted in the previous section, computations commenced at the day of the first lockdown order. The day of the first lockdown order in the four scenarios was taken to be: New York City – March 17, 2020 ^[13], New York State – March 17, 2020 ^[14], Diamond Princess cruise ship – February 5, 2020 (^[15,16]) Harris County, Texas – March 24, 2020 ^[17]. Different time intervals were considered for the different scenarios. For New York City, New York State, and the Diamond Princess cruise ship, simulations were performed until the first maximum in the profile for new infections per day, averaged over 7 days. The duration was 15 days for New York State, 20 days for New York City, and 12 days aboard the Diamond Princess. For Harris County, Texas, the scenario of interest begins on May 31, 2020, roughly the beginning of the surge in cases after Memorial Day in 2020. The Harris County scenario ends 26 days later, at the first large peak in new infections (averaging over 7 days) after Memorial Day. Infection data used to derive the baseline spread function was obtained from the following published online databases: New York City ^[18], New York State (^[19,20]), Diamond Princess cruise ship ^[21], Harris County, Texas ^[22]. For all of the scenarios, we investigated the effect of: ⅰ) increasing the outward filtration efficiency of the mask used by the infected population; ⅱ) increasing the level of compliance for mask usage, while maintaining the baseline filtration efficiency; and ⅲ) increasing the level of social distancing. For each of these three studies, the baseline spread function was determined using Eq. (1c) informed by published infection data for the four scenarios. The baseline spread function was then modified using Eq. (2a) to consider different filtration efficiencies, Eq. (2b) to study different compliance rates, and Eq. (2d) to see the effect of changes in social distancing.

The baseline outward filtration efficiency was taken to be 67% (^[8,23]). In the modified scenarios, higher-efficiency masks with outward FE's of 75%, 80%, and 90% were considered. Only outward efficiency was modified, i.e., the ability of the different masks to perform source control was enhanced.

The baseline compliance rate, which varied with time as more people adopted mask use, is determined from Eq. (2c) once the fraction ϵ_κ is prescribed. We took ϵ_κ to be 1/5, essentially saying that changes in droplet transmission $\widetilde {\beta }$ (weight of 1- ϵ_κ) is four times as important as the droplet production $\kappa$ (weight of ϵ_κ). The transmission $\widetilde {\beta }$ is influenced primarily by social distancing and the production $\kappa$ primarily from the use of PPE. The effect of different choices of ϵ_κ is discussed in ^[8]. Here we comment that the results are not highly sensitive to the choice of ϵ_κ, provided ϵ_κ is small compared to 1. Increases to the baseline compliance rate were affected by adjusting the factor F in Eq. (2c) so that the final compliance (at the end of the time interval) was 50%, 60%, and 70%.

Changes in social distancing were implemented in a gradual manner, as follows. A hypothesized decrease in the relative level of social contact (Eq. 2d), i.e. an increase in social distancing, was gradually applied through a temporal factor of the form (δ_b(t)/δ₀)^x, where δ_b(t) is the baseline spread function and x is an exponent that produces the target value L_mod/L_b at the end of the time window. That is, (δ_b(t_end)/δ₀)^x = L_mod/L_b. L_mod/L_b values of 90%, 80%, and 70% were considered which are representative of a 10%, 20%, and 30% increase in social distancing, respectively.

2.5.   Uncertainty analysis

The method for ascribing uncertainty to our estimated infection quantities in described in detail in the Appendix. Essentially, the uncertainty is derived by evaluating the variation in the input parameters γ and R₀ and observing the level to which the uncertainty propagates to model outputs. The space of possible parameter values is sampled using a Monte-Carlo-based technique.

3.   Results

As with most infection-spread models, the uncertainty in the calculations can be significant. To minimize clutter due to overlap of error bars with neighboring plots, we provide uncertainties for a small number of time points, and sometimes for only some of the curves. Unless otherwise stated, the uncertainties provided can be taken as representative of the error at other time points and for other curves in the given figure.

In the plots of numbers of daily new infections, values are normalized by the number at the beginning of the scenario. For New York City, New York State, and Diamond Princess, this also coincides with the initial time of the computations. As noted above, for Harris County, the beginning of the scenario of interest was approximately 70 days beyond the initial time for the computations.

3.1.   Baseline dynamics

The function δ(t)/Δ², which is the spread function in dimensional form, is plotted in Figure1a for all four of the scenarios considered. The dimensional form is provided for ease in interpreting (in Section 4) the role of the spread function δ(t).

The rates of new infections associated with these spread functions are plotted in Figure 1b. The rates are normalized by the number of new infections occurring at the beginning of the scenario.

3.2.   Population compliance in mask usage

Using the spread function for the four scenarios in Eq. (2c), and assuming an average outward filtration efficiency of 67%, the rate of mask adoption over time was estimated. The results, plotted in Figure 2, indicate a compliance of about 40% in New York State at the end of the first peak in new infections. At a given instant of time, the calculated compliance in New York State was about 1.3 times that of New York City, though at the end of each scenario (the time of peak number of new infections), the ratio was about 1.2 For the Diamond Princess, the calculated compliance at the end of the first peak (~ 20%) was about half of the calculated compliance for New York State at the end of its peak. In Harris County, the gradual acceptance of facemasks apparent in the other scenarios had already occurred prior to the scenario, and the compliance consequently varied more slowly and non-monotonically (around a mean of about 20%) for the duration of the scenario of interest.

3.3.   Effect of increasing outward filtration efficiency

Increasing the FE from the baseline 67% to 90% reduced the calculated number of new infections (at the end of the time window) by about 24% for New York State, 40% for New York City, 15% for the Diamond Princess, and 19% for Harris County (Figures 3a, b, c, d).

3.4.   Effect of increasing compliance

In Figures 4a, b, c, d, the computed number of new infections per day (at the end of the time window) decreases by about 50% (relative to baseline) in New York State, 75% in New York City, 70% aboard the Diamond Princess, and 84% in Harris County, when the compliance is increased from baseline to 70%. The large effect in moving from baseline to 50% compliance in the Diamond Princess and Harris County scenarios is due to the relatively low baseline compliance. Increasing compliance to 70% decreases the computed time of peak number of new infections, to about half its baseline value for New York State (i.e. from day 15 to day 7) and New York City (i.e. from day 20 to day 9). For the lower-compliance Diamond Princess and Harris County scenarios, the peak is reduced less significantly, between 2 and 5 days. The 70% compliance also reduces the number of new infections at the end of the scenario close to the initial value of 1.0.

3.5.   Effect of enhanced social distancing

The change in calculated daily new infections corresponding to increases in social distancing (reduced social contact) of 10%, 20%, and 30% are shown in Figures 5a, b, c, d.

For the 30% increase in social distancing, the calculated number of new infections decreased by about 56% in New York State, 76% in New York City, 63% aboard the Diamond Princess and 82% in Harris County. The predicted day of peak new infections dropped to about a third of its baseline value in New York State (i.e. day 15 to day 6) and New York City (i.e. day 20 to day 8), and decreased less than 10% (about 2 days) aboard the Diamond Princess.

4.   Discussion

Compared with other SIR-type models, the dynamic-spread model possesses an advantage related to the uniqueness of the result when a parameter changes value. To illustrate mathematically, we consider the number of new infections per day U to be a function of a parameter p (e.g. the filtration efficiency of a mask):

Incrementing the baseline parameter value p₀ to p₀ + dp to simulate an enhanced intervention strategy, and assuming a relatively small change (dp/p₀ ≪ 1), gives

or

Different SIR models will give comparable values for the baseline number of new infections per day, $f\left({p}_{0}\right)$, since they are all calibrated with published infection data. However, the parameter set selected (e.g. by a least-squares algorithm, or a "best guess") in order to achieve proximity to the published infection rates is non-unique. Depending upon the values of the parameters, and to some extent the functional form of the parameter in the model, different SIR models that agree closely on $f\left({p}_{0}\right)$ can yield significantly different values of $\frac{\partial f}{\partial p}\left({p}_{0}\right)$. As an example, in ^[5], the filtration efficiency is multiplied by an effective contact rate (as well as other parameters and independent variables), which must be prescribed. Different prescriptions yield different values of dU (Eq. 4c). We remark that even low-dimensional systems, containing 3 or 4 equations (the present model is derived from the 4-equation model in ^[9]), usually contain clusters of parameters composed of numerous other parameters that must be prescribed. The model in ^[9], for example, requires that roughly (depending upon the droplet size of interest) 20 parameters be determined. With the dynamic-spread model, the number of new infections per day U can considered to be a function of the spread-function δ (i.e. U(δ)). Then dU is (using the chain rule) given by

where pagain is a parameter such as filtration efficiency or level of social distancing. The quantity dU/dδ is determined from the governing equations (1a, b, c) informed by published infection time traces. No non-uniqueness is introduced by the process of determining dU/dδ. The ratio dδ/dp is the dependence of the spread function upon the parameter of interest. This is typically a linear dependence that arises naturally, e.g. the spread function depends linearly upon the level of social distancing and filtration efficiency of the barrier ^[8]. While there is still uncertainty associated with predictions (associated with uncertainties in the published values of reproduction number R₀, the exponential grown rate M, and the recovery rate μ_Ι) using the dynamic-spread model, the systematic approach of the dynamic-spread model avoids the non-uniqueness problems associated with many SIR models. As discussed in ^[8], the dynamic-spread approach also naturally simulates the continuous process by which populations change behavior. No decision need be made regarding when to update parameters (i.e. p₀ to p₀ + dp in Eq. (4)) to simulate changes in population behavior, as is the case with most SIR-based models.

Actual measurements of mask compliance in the locations considered, which would be valuable for comparing with the results in Figure 2, are difficult to perform. One observational survey conducted in New York City ^[24] at the end of July 2020 found that 75% of persons surveyed at 14 locations within the city were wearing masks. The compliance varied widely between locations, from 20% to 99%. The study organizers indicated that the sampling locations were not chosen at random and were not necessarily representative of the whole city. They noted that the city-wide mask-wearing prevalence was probably lower. Another survey ^[25], this one self-reporting, found that 67% of respondents within New York City always wore masks when going out in public. For the United States as a whole, the value was 48% to 51%. A second nationwide poll, performed by Gallup ^[26], found that 44% of people across the United States claimed to wear masks in public all the time. To compare with the survey results specific to New York City, our simulation in the New York City scenario was continued two additional months, to the end of July (Figure 6)

As shown in Figure 6, calculated mask compliance plateaued at approximately 38% about one month into the epidemic. Variations around the 38% value are likely due to reporting uncertainties and the method of averaging employed (7-day in this case) and are likely not significant. One reason why the calculated compliance may be low compared with survey values is the relatively small value (1/5) of ϵ_κ. As explained above, this value represents the fraction of the temporal change in the spread function due to changes in mask-wearing behavior. Intuitively, we expect increased social distancing to decrease the spread function more intensely than improved barrier usage, especially when social distancing includes complete lockdown. Thus, we make ϵ_κ less than 1/2. However, the value of the fraction is somewhat uncertain. Increasing ϵ_κ to 1/3 gives a compliance of about 62% in New York City in July of 2020. This is comparable to the survey values, given the uncertainties associated with the computations and the surveys. The uncertainty in the calculated compliance (Figure 2) is on the order of 15%. Another important factor relevant to the comparison of compliance values is the filtration efficiency. The calculated compliance is based upon an outward filtration efficiency of 67%. While even many homemade mask designs can achieve this level in laboratory tests ^[27], the calculations represent an average over all times, including when masks might not be worn. These times include situations when the user is singing, conversing while eating or drinking, or taking relief from mask fatigue. During such situations, the filtration efficiency is zero, effectively lowering the laboratory value. Thus, an FE lower than 67% may be more appropriate for comparison with measurements. When detailed information regarding the average filtration efficiency of barriers deployed by the community is available, along with the compliance for the barrier use, Eq. (2c) can be used to prescribe ϵ_κ.

The rate of spread of infection and the effectiveness of countermeasures are a function of both the spread function and its derivative. A large value of the spread function, e.g. New York City (Figure1a), is associated with a high infection rate (Figure1b). For New York state as a whole, the spread function was lower (Figure1a) and the per-capita number of new infections was lower (Figure1b). The exponential growth rate M (Section 2.3) was determined to be 0.46/day in New York City, compared with 0.36/day in New York State. In epidemiological terms, the results indicate that the residents of New York City engaged in more behaviors conducive to the spread of infection than did the residents of the state in its entirety. Regarding the derivative of the spread function, even a small positive value is associated with a rapid increase in infection rate. In Harris County, the spread function is relatively flat until around day 4 (June 3, 2020), due to measures taken in the months prior to the Memorial Day weekend. After opening of businesses and social gatherings around Memorial Day, the spread function began to increase. The result was the large increase in new infections per day visible in Figure1b. For the Diamond Princess scenario, a slight positive derivative of the spread function also occurred, commencing around day 7. A rapid increase in new daily infections (Figure 1b) matched this increase in the value of the spread function. Enhanced interventions that convert the spread function from increasing to decreasing are predicted to have a profound effect on the infection rate. For example, increasing the level of social distancing in the Harris County scenario by 10% nearly halved the number of new infections, and increasing social distancing by 30% completely eliminated the infection surge after Memorial Day (Figure 5d). For the Diamond Princess, the 30% increase in social distancing is predicted to cut the number of new infections in half (Figure5c). However, the intervention must be strong enough to convert the derivative of the spread function from positive to negative. Increasing the filtration efficiency to 90% had a relatively minor effect in the Harris County (Figure 3d) and Diamond Princess (Figure 3c) scenarios. This modest decrease is due to the fact that mask usage is a weaker countermeasure than social distancing (ϵ_κ ≪1), and the mask compliance was relatively low (on the order of 20%) in both the Harris County and Diamond Princess scenarios (Figure 2). An increase of compliance to 50% or more in the two scenarios would produce a much more profound effect (Figures 4c, 4d); however, this would require a substantial change in population behavior. In cases where the derivative of the spread function is negative and large, i.e. in cases where the population responds quickly to mask recommendations/orders, enhanced interventions can have a pronounced effect. An example is the 40% reduction in number of new infections in New York City when the filtration efficiency is increased from baseline (67%) to 90% (Figure 3a.)

The insights just presented regarding the infection dynamics for the four specific scenarios can be generalized to the following statements about the dynamics-spread-function approach:

1) Scenarios in which the spread function is large will be associated with large infection rates.

2) When the derivative of the spread function is large and negative (population highly adaptive to proposed countermeasures), enhanced intervention strategies will have a pronounced effect.

3) When the derivative of the spread function is positive, the infection spread is rapid.

4) Enhanced countermeasures that change the derivative of the spread function from positive to negative will produce a significant decrease in the predicted infection rate.

5) The spread function can be used to compute mask compliance throughout the entire period for which calibration data (number of new infections, total infected population) is available.

6) When the computed compliance is low (e.g. below 20%), enhanced mask filtration efficiency will be unlikely to change the sign of the spread-function derivative, and the reductions in infection rate will be modest.

7) Statements 1–6 do not involve new simulations of the infection dynamics (e.g. the affected populations), only computation of the spread function for the time interval of interest. The spread function is thus a useful dimensionless parameter for understanding, and predicting the general features of, the infection profile when new intervention strategies are contemplated, without having to perform calculations.

We reiterate the caveat made in ^[8] that the dynamic-spread model is not a forecasting model, in terms of predicting the course of an ongoing epidemic assuming a given type of intervention. The dynamic-spread model is informed by the infection profile for a scenario that has already occurred. The infection rate depends upon the specific: pathogen, population, time interval, set of environmental conditions, and intervention strategy. All of these influencing factors are accounted for, indirectly, in the spread function. The model can then estimate the modified course of the infection for that same scenario, assuming a different intervention strategy. The tradeoff for lack of forecasting capability is the ability to capture the baseline population behavior retrospectively in a precise and continuous way, which enables studies of hypothetical intervention strategies to proceed in a natural manner. The resulting information can be subsequently used for designing interventions to future scenarios.

5.   Conclusion

The dynamic-spread model was used to compare enhancements to baseline intervention strategies in four COVID-19 scenarios. The effectiveness of the same enhancement was strongly dependent upon the scenario involved. When the average outward filtration efficiency for facemasks worn in New York City during the first wave of COVID-19 was increased from an assumed baseline of 67% to a hypothesized 90%, the calculated peak number of new infections per day decreased by 40%. For the same baseline and hypothesized filtration efficiencies aboard the Diamond Princess cruise ship, the calculated peak number of new infections per day decreased by about 15%. An important factor contributing to the difference between the two scenarios is the lower mask compliance (derivable from the spread function) aboard the Diamond Princess. The dynamic spread function can predict the gross features of the infection profile when alternative intervention strategies are implemented, without requiring computation of the relevant populations. The dynamic spread function is thus a useful tool for enabling public-health officials and policy makers to estimate the benefits of a certain type of intervention, for their specific population.

Conflicts of interest

The authors declare there is no conflict of interest.