Abstract
Retrospective analyses of interventions to epidemics, in which the effectiveness of strategies implemented are compared to hypothetical alternatives, are valuable for performing the cost-benefit calculations necessary to optimize infection countermeasures. SIR (susceptible-infected-removed) models are useful in this regard but are limited by the challenge of deciding how and when to update the numerous parameters as the epidemic progresses. We present a method that uses a “dynamic spread function” to systematically capture the continuous variation in the population behavior throughout an epidemic. There is no need to update parameters as the effects of interventions are gradually manifested in the infection dynamics. We use the tool to quantify the reduction in infection rate realizable from the population of New York City adopting different facemask strategies during COVID-19. Assuming a baseline facemask of 67% filtration efficiency, calculations show that increasing the efficiency to 75% could reduced the roughly 5000 new infections per day occurring at the peak of the epidemic by about 40%. The turn-around time for the epidemic decreases from around 37 days to 31 days. Mitigation strategies that may not be varied as part of the retrospective analysis, such as social distancing, are automatically captured as part of the calibration of the dynamic spread function.
1. Introduction
Retrospective analyses of strategies used to contain epidemics such as COVID-19 are valuable for countering seconds waves of the infection, selecting countermeasures for future epidemics, and educating the population regarding the efficacy of implementing certain behavior modifications. In particular, public-health agencies responsible for recommending types of personal protective equipment (PPE) to stockpile in anticipation of a future pandemic can benefit greatly from the cost-benefit information yielded by retrospective analyses. Mathematical models, including those of the SIR type, can be helpful in providing a quantitative framework for the analysis. SIR models have been applied during the COVID-19 pandemic (Stutt et al. 2020, Giordano et al. 2020, Bertozzi et al. 2020, Cooper et al. 2020), primarily in a predictive capacity. To our knowledge, no studies have attempted to re-create an actual scenario (e.g. CNYC) with different interventions, though in Cooper et al. (2020) the influence of different mask types on the reproduction number was computed using a generic infection scenario. A formidable challenge in applying SIR models is prescribing the values of the numerous parameters, and updating them to simulate changing infection dynamics, e.g. in response to interventions. The challenge is accentuated by the high sensitivity of the predictions to some of the parameter values (Giordano et al. 2020). Often, parameter choices are based upon best guesses, or closeness of fit (sometimes visual) of computed profiles with published curves (Cooper et al. 2020). In this paper we introduce a modification of traditional SIR models that incorporates a “dynamic spread” function that captures changes in population behavior in a continuous manner. There is no need to adjust parameters as interventions are implemented during the course of the infection. The spread function is calibrated using published infection curves for a specific epidemic, then modified to analyze the effect of alternate intervention strategies. We illustrate the process using the COVID-19 crisis in New York City (CNYC). The reduction in infection rate realizable in CNYC with an increased level of mask usage, and with deployment of masks offering higher levels of filtration, is estimated.
2. Methods
We illustrate the technique using a 3-equation SIR model. It assumes that the infection dynamics are dominated by one transmission mode (e.g. airborne particulates), and that the parameter values are appropriate for all particle sizes contributing to that mode (though interpretation of the resulting equations as an average for a broad particle distribution is possible (Myers et al. 2016)). More complicated SIR models can be useful, particularly if it is desired to model the details of the infection dynamic, e.g. symptomatic and asymptomatic individuals (Stutt et al. 2020). Our intention is to use the simplest model that can capture the known infection dynamics in a general sense. As noted by Siegenfeld et al. (2020), simpler models can prove more useful than complex ones, in part because accurate data is often not available to inform complicated formulations. We also note that much of the technique we present can be extended to more complex models.
2.1 Overview of Strategy
The dynamic-spread function is part of a systematic procedure for calibrating SIR models. The 5 steps in the procedure are listed below and implemented subsequently.
Use the rate of change (measured by the number of new infections per day), dS/dt, of the susceptible population, as the variable for calibrating the model, rather than S. The derivative profile, which we call T(t), was felt to be the quantity known most accurately. Its determination does not require the number of recovered patients to be tracked.
Normalize variables and identify critical dimensionless parameters. Formulating the model in terms of dimensionless quantities reduces the number of property values that must be determined, which can be numerous in SIR models (Stilianakis and Yannis Drossinos 2010).
Allow the dimensionless parameter δ, which is essentially the product of the infection transmission rate and the virus production rate, to vary with time, and account for its time dependence in the governing differential for T(t). We denote δ(t) the “dynamic-spread” function, as it contains the elements that both vary with time and govern the rate of spread of the infection.
Derive the governing equation for δ(t). Provide the required coefficient functions using published T(t) profiles for a given infection scenario.
Designate as the time origin for the dynamic analysis the point of the first intervention into the epidemic. For CNYC, we identify this as day 17, when shelter-in-place was instituted. Prior to that point, it is assumed that δ is constant in time, and a standard SIR model applies. The parameters for the standard SIR model can be estimated from the published growth rate and reproduction number. The resulting values serve as initial conditions for the dynamic analysis.
2.2 Development of Governing Equations
The SIR equations under these conditions listed in sections 2 and 2.1 are (Stilianakis and Drossinos 2010, Myers et al. 2016) as follows. The primes denote dimensional quantities and will be dropped following nondimensionalization.
Here S’ is the number of susceptible individuals in the total population N, I’ the number of infected individuals, D’ the total number of droplets, 
the transmission rate, μI the infection recovery rate, κ the droplet production rate, and v−1 the droplet removal rate. We next introduce the maximum number of newly reported infections (roughly 5000 per day for CNYC) α, and a time scale Δ, which we take to be the turn-around time (time required to reach dT’/dt’ = 0, about 37 days for CNYC.) The function T’ is scaled by α, I’ by αΔ, and D’ by καΔ2 Additionally, we differentiate Eq. 1a (after nondimensionalization), allowing the spread parameter
to vary with time. Also, for huge populations as in CNYC, we ignore terms proportional αΔ/N, including S’/N – 1. Under these assumptions, the SIR equations become:
where γ = ΔμI is the dimensionless infection recovery rate and 
is the dimensionless droplet decay rate. Inserting (2d) in (2b) and using (2a) yields the 2-equation system:
As noted above, we take the time origin to be the time of first intervention.
To determine the dynamic spread function, we reformulate Eq. (3a) as an equation for δ(t), assuming T(t) and I(t) to be known. We then obtain the T(t) profile from the published number of new infections per day (Johns Hopkins Resource Center 2020) in the locale of interest (e.g. New York City). We label this published profile Tp(t) and the resulting (from Eq. (3b)) infection profile Ip(t), and we insert them into the equation for δ(t). The resulting equation for the dynamic spread function is:
Because the governing equation for δ(t) is informed by the published Tp(t) profile, solving Eq’s 3 using this dynamic spread function will exactly reproduce the published Tp(t) curve. The utility of δ(t) derives from modifying it to model alternative intervention strategies and solving Eqs (1) to determine the modified infection rate. Modifications to account for protective strategies were performed in the following manner.
2.3 Accounting for Protective Equipment
We build upon a previously developed SIR model (Myers et al. 2016, Yan et al. 2018) that systematically accounts for the presence of protective equipment. Differentiating Eq. (1d) with respect to time yields
Apportioning a fraction ϵκ (e.g. ½) of the change in δ to changes in droplet production, we set
Since from (1d)
which can be integrated to
In Myers et al (2016), it was shown that the production rate in the presence of protective equipment can be written as
Here FE is the filtration efficiency (e.g. the FE for an N95 respirator is 95%) of the mask for the dominant droplet size, and fi is the fraction of the infected population wearing the covering at any given time. To perform a retrospective analysis in which a barrier material of different capturing efficiency is investigated, the new FE value would be used in Eq. (10) which, with Eq. (9), would be used to create a new dynamic-spread function. The modified spread function would then used (in Eqs (3)) to estimate the change in infection rate.
To analyze scenarios where different fractions of the infected population deploy a given mask type, (9) and (10) can be combined to give
and
Eq. (12) can be used to estimate the fraction of infected individuals deploying a mask of a given filtration efficiency for a baseline case (known δ(t)). The effect of different fractions of the infected population deploying the mask can be quantified by prescribing fi(t), solving for δ(t), and using this modified spread function in Equations (3).
2.4 Solution Technique
The initial conditions for Equations (3) are obtained by simulating the dynamics of the infection prior to any intervention, days 1 – 17 for CNYC. In that case, the derivative of the spread function is zero and Equations (3a,b) revert to a traditional SIR model. Seeking solutions that have an exponential time dependence of the form exp(Mt) result in the algebraic equation
The subscript “0” on δ implies that the value applies to the initial period of the infection, before intervention occurs. The other parameters do not vary during the course of the epidemic and are not subscripted. An exponentially growing solution will occur when
The ratio δ0/(γ λ) is the reproduction number R0 (Myers et al. 2016) for the standard SIR model. The growth rate M can be obtained from infection rates published during the beginning of the infection. Estimates of the reproduction numbers R0 for the early stages of epidemics are also published. In the simulations, a range of recovery times μI (dimensionless recovery times γ) ranging from 2 days to 10 days were considered. For any given value of γ, λ and δ0 were obtained from Eqs (13) and (14) and used in the solution of the dynamic equations (3). The initial value for T(t) was obtained from the published profile Tp(t) (published number of new infections per day, Johns Hopkins Coronavirus Research Center, 2020), evaluated at day 17. Ip(t) was derived from Tp(t) using Eq. (3b), rather than using a published infection profile, so that it was not necessary to ascertain how well recoveries were tabulated in the published infection curves. The Ip(t) function evaluated at day 17 was used to provide the initial condition for I(t). Equations (3) were solved using a Runge-Kutta method (Matlab ode45, Mathworks Inc.).
3. Results
We performed a retrospective analysis of CNYC during days 17 - 37. This interval was chosen because day 17 is the day of the first intervention (shelter in place), and day 37 is the time of maximum new infections per day, based upon a 7-day average (Johns Hopkins Coronavirus Research Center, 2020). Initial reproduction numbers between 2 and 10 were considered, along with recovery times between 2 days and 10 days. The fraction ϵκ was ½. We analyzed scenarios where the infected population deployed different types of masks. For baseline, it was assumed that the FE was 67%. This value is representative of homemade masks (Howard et al. 2020), though the filtration capability of homemade masks spans a wide range. Higher-efficiency masks with FE’s of 75%, 80%, and 90% were considered for the hypothetical scenarios. Uncertainty in the calculated results was obtained by performing simulations for different combinations of reproduction numbers and recovery rates and computing the mean and standard deviation for the ensemble of parameter combinations.
Figure 1a shows the dynamic spread function as a function of time. A sharp decrease is seen initially, owing to the shelter-in-place restriction. For larger FE, a sharper decrease in the spread function is observed. A marked decrease in new infections results from an increase in FE (Fig. 1b). Increasing FE from 67% to 75%, for example, reduces the maximum number of new infections per day by about 40%, and decreases the turn-around time from 37 days to about 31 days. The number of infected individuals (Fig. 1c) is reduced by about 30%. The uncertainty is considerably larger for the spread function (Fig 1a) and the infected population (Fig. 1c) than the number of new infections per day (Fig. 1b), as the spread function and infected population are much more strongly influenced by the recovery time than the number of new infections. The recovery time was allowed to span a factor of 5.
Dynamic spread function (1a), change in new infections (1b), and number of infections (1c), for CNYC with infected population deploying masks with different filtration efficiencies. Shadowed regions denote values within a standard deviation of the mean.
An additional set of simulations was performed in which the fraction of the CNYC population wearing a mask (FE=67%) was increased, by 2%, 5%, and 10%. The baseline fraction was determined from Eq. (12), with δ(t) (baseline curve in Fig. 1a) derived from Eq. (4). The fraction deploying masks increased from 0% at day 17 to roughly 75% on day 37. Increasing the fraction by 10% (for all times) reduced the number of new infections per day from about 5100 to 3600 (Fig. 2). The turn-around time is reduced from approximately 37 days to 32 days.
Number of new infections per day for scenarios where various fractions of the infected population in CNYC deploy masks with a 67% filtration efficiency.
4. Discussion
While days 17 to 37 were featured in our simulations, the dynamic-spread-function technique can be applied to any period where reliable numbers of new infections are available. Calibration of the model is based upon number of new infections, rather than the size of the infected population, to eliminate the uncertainty associated with how well recoveries are tracked in the calibration dataset.
Since the dynamic spread function is the product of the transmission rate times the droplet production rate, it quantifies the ability of the infection to spread. The ability of the infection to spread decreased rapidly from day 17 (Fig. 1a). Because adoption of masks in CNYC likely occurred on a continuous basis over the weeks succeeding day 17, the curves in Fig. 1a are smooth and monotonic. For more abrupt changes in behavior, the spread function need not be monotonic.
The salutary effects of protective equipment in CNYC were investigated without having to update the SIR parameters during the epidemic. The continuous adoption of masks would be difficult to simulate by updating coefficients in standard SIR models. With the dynamic-spread approach, the gradual adoption of masks is captured very naturally.
The dynamic-spread approach also allowed other behaviors affecting the infection dynamics, such as social distancing, to be captured without being specifically modeled. The pattern of social distancing in CNYC was retained for all the simulated scenarios involving different facemasks. This commonality is largely responsible for the similar shapes of the curves in Fig. 1a. The only decision made relative to social distancing was that the factors contained in the transmission rate 
, which includes social distancing (Myers et al. 2016), were responsible for roughly half (ϵκ = ½) of the reduction in the spread function (shown in Fig. 1a). Other fractions would result in different numbers in Figs 1 and 2, with higher values of ϵκ resulting in larger reductions in infection rate, and vice versa.
The simulated scenarios addressed only changes in masks worn by the infected population. No change in protection for the susceptible population was assumed. The susceptible population deployed masks, but the type was not varied between scenarios. As shown in Myers et al. (2016), the effect of mask deployment by the susceptible population can be simulated by modifying the transmission term 
in Eq. (3), using equations analogous to (9) and (10).
For the conditions of the simulations, a slight increase in facemask efficiency resulted in a larger benefit than a commensurate increase in compliance. At day 37, for example, a fractional increase in compliance of 0.1 resulted in a reduction in new infections of about 1500 per day (Fig. 2), while a fractional increase in FE of 0.1 reduced the number of new infections by over 2000 (interpolating Fig. 1b). For a higher baseline FE than 67%, increasing the compliance rate would produce a larger decrease in new infections. This comparison between filter efficiency and population compliance illustrates the utility of the model for determining how resources devoted to countermeasures can be optimally spent. In this case, the model can help inform the choice between 1) producing and distributing barriers of higher FE compared, and 2) educating and incentivizing the population to deploy barriers more readily available.
A noteworthy conclusion emerging from the simulations is that considerable benefit can be obtained from higher FE masks without requiring N95 levels of efficiency (Fig 1). It is important to emphasize that for the benefits to be realized, the filtration efficiencies for the barrier material must be attainable for the particle-size range of the dominant transmission mode for the given scenario. One way of assuring this is for the barrier to provide the given FE across the spectrum of particle sizes. Otherwise, knowledge of the material filtration efficiency for the intended application (e.g. reducing airborne particulates generated by coughing or sneezing by infected persons indoors) is required in order to generate useful estimates. The complex issues of dominant transmission mode for COVID-19, and the filtration efficiency of different masks designs for the different modes, will be addressed in future applications of the model.
Data Availability
No human, animal, or in-vitro data is included in the manuscript, which describes a risk-assessment model. Algorithms used to derive the mathematical results have been referenced.
