A mathematical model for assessing the effectiveness of protective devices in reducing risk of infection by inhalable droplets

A mathematical model for assessing the effectiveness of protective devices in reducing risk of... Abstract Respiratory protective devices (RPDs) are critical for reducing the spread of infection via inhalable droplets. In determining the type of RPD to deploy, it is important to know the reduction in the infection rate that the RPD enables for the given pathogen and population. This paper extends a previously developed susceptible-infected-recovered (SIR) epidemic model to analyse the effect of a protection strategy. An approximate solution to the modified SIR equations, which compares well with a full numerical solution to the equations, was used to derive a simple threshold equation for predicting when growth of the infected population will occur for a given protection strategy. The threshold equation is cast in terms of a generalized reproduction number, which contains the characteristics of the RPDs deployed by the susceptible and infected populations, as well as the degree of compliance in wearing the equipment by both populations. An example calculation showed that with 50% of the susceptible population deploying RPDs that transmit 15% of pathogens, and an unprotected infected population, an otherwise growing infection rate can be converted to one that decays. When the infected population deploys RPDs, the transmission rate for the RPDs worn by the susceptible population can be higher. 1. Introduction Respiratory protective devices (RPDs) such as respirators constitute an important element of the control strategy against the transmission of infectious diseases. In determining whether to deploy a given RPD to protect a specific population against a certain pathogen, including whether to allocate substantial resources to stockpile RPDs against a possible future threat, realistic estimates of the reduction of risk associated with use of the RPD in the scenario of interest is important. In evaluating the reduction of risk associated with deployment of a given RPD, important parts of the assessment include determination of the intrinsic penetration of the device (Committee on the Certification of Personal Protective Technologies 2011; Mniszewski et al., 2014), which is governed by the microscopic properties of the barrier and the amount of leakage due to imperfect fit (Coffey et al., 2004; Lee et al., 2005) between a particular facial profile and the RPD. Another important part of the risk assessment for RPDs is an estimate of the change in the infection rate associated with the change in pathogen transport introduced by the presence of the barrier. The change in pathogen transport could apply to either an uninfected person attempting to reduce the intake of pathogens or an infected person attempting to reduce the output. While numerous models of various types have been developed to simulate the spread of epidemics (Furuya, 2007; Chen & Liao, 2008; Coburn et al., 2009; Sze To & Chao, 2010), to the authors' knowledge no systematic treatment of the effect of RPDs has been incorporated into the models. Of particular interest for RPD evaluation is the transmission of disease by inhalable respiratory droplets. A model that is well suited to estimating the spread of disease by inhalable respiratory droplets is the susceptible-infected-recovered (SIR) epidemic model. In the SIR models, the influence of RPDs can be accounted for in a systematic manner in the parameters within the system of differential equations. Stilianakis & Drossinos (2010) developed an SIR model for disease transmission with inhalable droplets as the transmission vector, assuming a closed, homogeneously mixed population. Epidemics mediated by both airborne and settled droplets were considered, and the time course for the susceptible, infected and recovered populations derived using parameters derived from experimental evidence. Relevant parameters of the model include the droplet production rate, gravitational settling rate, deposition probability and number of pathogens per droplet, all of which are a function of droplet size. Stilianakis & Drossinos (2010) do not address the effect of protective equipment. Of considerable importance in evaluating risk in the presence of protective equipment is the polydisperse nature of the droplets. Even without considering protective equipment, distinguishing between different sizes can be important, as studies have found that cough droplets possess a size distribution ranging often from less than 1 μm to over 100 μm (Chao et al., 2009). Nicas et al. (2005) estimated that the droplet production rate at 4 μm is over 20 times that at 8 μm. When protection is considered, the strong dependence of RPD performance on size must be taken into account. Huang et al. (2007) showed that the penetration of 1 micron particles for some respirators could be an order-of-magnitude higher than the penetration at 5 microns. In this paper we modify the SIR theory as presented by Stilianakis & Drossinos (2010) to account for the influence of RPDs. We account for the fact that only a fraction of the population will likely deploy protective measures by dividing the susceptible population into two groups, one of which deploys RPDs and one which does not. Similarly, a fraction of the infected population utilizes RPDs, thereby reducing the source of pathogens. We first implement our modifications under the idealized assumption of a monodisperse aerosol distribution, and then extend the formulation to the more realistic case of a polydisperse distribution. 2. SIR Model in the presence of protective equipment 2.1 Governing equations for a monodisperse aerosol The original SIR model equations for a monodisperse aerosol, in the absence of protective equipment, are (Stilianakis & Drossinos, 2010):   dSdt=−β˜DSN, (2.1a)  dIdt=−dSdt−μII, (2.1b)  dDdt=κI−1νD, (2.1c) where   β˜=cBVclτpqNp (2.1d) and   v−1=(1+cτ)BVclq+μd+θ. (2.1e) 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, μI the infection recovery rate, κ the droplet production rate, c the contact rate between a susceptible person and an infected person, B the breathing rate, Vcl the volume of the personal cloud of an infected person, τ the characteristic breathing time, p the probability of infection by an inhaled pathogen, q the inhaled-droplet deposition probability, Np the number of pathogens per droplet, μd the inactivation rate of airborne pathogens, and θ the gravitational settling rate. The boundary conditions are   S(0)=S0, (2.1f)  I(0)=I0, (2.1g)  D(0)=0. (2.1h) Equations (2.1) describe the evolution of the susceptible and infected populations, as a function of the relevant parameter groups. The processes described by the evolution equations are illustrated in Fig. 1a and b. The protected population will be introduced in the next section; for now this population can be ignored, as can minor notational differences between the figure and (2.1). (Notation in Fig. 1 corresponds to that used in the RPD model below.) The susceptible population converts to the infected population with a time rate of change proportional to β˜DS. The infected persons convert to the recovered population at a rate proportional to μII.The total rate of change of infected individuals (2.1b) is equal to the rate at which susceptible persons convert to infected minus the rate at which infected persons recover. The total number of droplets, described by (2.1c) and shown in Fig. 1b, increases by one mechanism and decreases by four mechanisms. The generation mechanism is production by infected persons, with a rate proportional to κwI. The attenuation mechanisms are: (1) inhalation of an infected person's own droplets (rate proportional to BqD/Vcl), with the droplet ultimately residing with the infected person generating the droplet; (2) inhalation of someone else's droplet (rate proportional to cτ BqD/Vcl), with the destination of the droplet being potentially any of the three populations; (3) loss of infectivity (rate proportional to μDD), with the droplet remaining in the global environment but not being infectious and (4) droplets settling onto a surface under the influence of gravity (rate proportional to θD). Fig. 1. View largeDownload slide (a) Flux of people between the susceptible, infected and recovered populations, along with parameters characterizing the rate at which people convert from one population to the other. Dashed line symbolizes new features of the dynamics introduced by protective equipment. (b) Schematic of droplet dynamics. D(t) is total number of droplets. Droplet generation mechanism is shown to the left of D(t); annihilation mechanisms are shown to the right. Box on the far left denotes origin of the drops; boxes on the far right denote destinations. Next to the mechanism arrows are parameters characterizing the rate of droplet creation or removal Dashed arrow indicates feature that is affected by protective equipment. Fig. 1. View largeDownload slide (a) Flux of people between the susceptible, infected and recovered populations, along with parameters characterizing the rate at which people convert from one population to the other. Dashed line symbolizes new features of the dynamics introduced by protective equipment. (b) Schematic of droplet dynamics. D(t) is total number of droplets. Droplet generation mechanism is shown to the left of D(t); annihilation mechanisms are shown to the right. Box on the far left denotes origin of the drops; boxes on the far right denote destinations. Next to the mechanism arrows are parameters characterizing the rate of droplet creation or removal Dashed arrow indicates feature that is affected by protective equipment. Disease transmission is initiated by a small number of infected persons ( I0≪S0) that have not yet generated a significant number of droplets at the time the disease progression begins to be tracked ( D(0)=0). At any time during the analysis of the disease progression, the recovered population can be derived from the equations S+I+R=N. That is, the total population is assumed constant. The SIR equations can be modified to account for the presence of RPDs in the following manner. We consider a susceptible population that contains a group Sr deploying RPDs, and a group Snr not deploying RPDs. The two groups comprise the total susceptible population, i.e.   S=Sr+Snr. (2.2) The transmission rate β˜nr for the no-RPD group Snris identical to β˜ in (2.1d). For Sr, the reduction in pathogen transmission introduced by the RPD is accounted for by the introduction of the factor Tin representing the fraction of incoming (to the respiratory system of the susceptible population) pathogens transmitted by the RPD when the susceptible person breathes in. For example, for an N95 respirator capable of filtering out 95% of the pathogens, Tin=5%. The transmission rate in the presence of respiratory protection can be written   β˜r=TincBVclτpqNp=Tinβ˜nr. (2.3) This equation can be thought of as a pathogen transmission rate proportional to a reduced breathing rate TinB or a reduced probability of infection Tinp. While the protected group of susceptible individuals ( Sr) converts to the next phase (infected) at a different rate than the unprotected ( Snr) group, a similar dynamic does not hold true for the infected population. The infected population not deploying RPDs ( Inr) converts to the next phase (recovered) at the same rate as those deploying RPDs ( Ir). The two groups produce pathogens at different rates, however, and the droplet product term can be written   κrIr+κnrInr. (2.4a) Analogously with the reduction in incoming pathogens associated with a barrier having an incoming transmission rate Tin, we can model the reduction in outgoing airborne pathogens introduced (through coughing, sneezing, breathing …) by an infected person deploying an RPD using an outward transmission rate Tout. The effective respirable droplet production rate for an infected person wearing an RPD is   κr=Toutκnr. (2.4b) The infected population changes because susceptible people convert to the infected group, and infected persons convert to the recovered group. To prescribe what proportion of the susceptible population deploys RPDs after becoming part of the infected population, at least two possibilities arise. The first is that the susceptible people continue to wear RPDs in the same proportion after becoming infected as they did before they were infected. A second possibility is that they deploy RPDs in the same proportion as the rest of the infected population. We adopt the second choice, for the following mathematical simplification. Since the Ir and Inr groups convert to the recovered population at the same rate, and the newly converted (from the susceptible population) infected persons mimic the behaviour of previously infected persons, the initial fraction of infected people deploying RPDs maintains throughout time. We denote this initial fraction of the infected population by fi and write the infected population as   I=Ir+Inr=fiI+(1−fi)I. (2.5) The production rate in (2.4) can similarly be written   [κrfi+κnr(1−fi)]I=κwI, (2.6a) where κw represents a ‘weighted’ droplet production rate. Using (2.4b), the weighted droplet production rate can also be expressed as:   κw=κnr[Toutfi+(1−fi)]. (2.6b) By employing the modified variables and parameters in (2.2)–(2.6) in the original SIR equations (2.1a)–2.1h, we can compose the SIR equations for the case where a fraction of the population deploys RPDs. After additionally normalizing S, I and D by N to work in terms of fractional populations, we obtain   ddtS^r=−β˜rD^S^r, (2.7a)  ddtS^nr=−β˜nrD^S^nr, (2.7b)  ddtI^=−ddt(S^r+S^nr)−μII^, (2.7c)  ddtD^=κwI^−1νD^, (2.7d) where   S^r=Sr/N, (2.7e)  S^nr=Snr/N, (2.7f)  I^=I/N, (2.7g)  D^=D/N. (2.7h) The boundary conditions are:   S^r(0)=S^r,0=Sr,0/N, (2.7i)  S^nr(0)=S^nr,0=Snr,0/N, (2.7j)  I^(0)=I^0=I0/N, (2.7k)  D^(0)=0. (2.7l) Here S^r,0 and S^nr,0 are the normalized susceptible populations initially deploying RPDs and not deploying RPDs, respectively. It is important to note that it is not precisely true that S^r,0+S^nr,0=1, since S^r,0 and S^nr,0 are normalized by the total population N, not the total initial susceptible population S^0. (The precise relation is S^r,0+S^nr,0+I^0+R^0=1, where R^0 can assumed to be 0 and I^0 to be small.) For convenience in later expressions, we introduce the fraction of the total initial susceptible population, fs, that deploys RPDs. Then (2.7i) and (2.7j) can also be written   S^r(0)=fsS^0, (2.7m)  S^nr(0)=(1−fs)S^0. (2.7n) The additional features of disease-transmission dynamics introduced by the presence of protective equipment are denoted with the dashed lines in Fig. 1a and b. As shown in Fig. 1a, protective measures introduce a new population, the protected susceptible population. This group converts to the infected population at a slower rate ( β˜r<β˜nr) than the unprotected susceptible population. The dashed lines in Fig. 1b signify that, with the availability of protective equipment, the infected population will generate droplets at a slower rate than in the absence of protective equipment (κw<κnr). 2.2 Approximate solution for a monodisperse aerosol We now analyse (2.7a), with the goal of identifying conditions under which the infected population decays with time due to the presence of the RPDs. We begin by noting that (2.7a) may be rewritten   ddt¯(lnS^r)=−β˜nrD^ (2.8) Integration of (2.8) and application of boundary condition (2.7e) yield   S^r(t)=S^r,0exp−(β˜r∫0tD^(t′)dt′) (2.9) Similarly,   S^nr(t)=S^nr,0exp−(β˜nr∫0tD^(t′)dt′) (2.10) Differentiating (2.7d) with respect to t and using (2.7c) for dI^/dt, with ddt(S^r+S^nr) provided by differentiating (2.9) and (2.10), results in the second-order differential equation   d2D^dt¯2+1νdD^dt¯=κwD^[β˜rS^r,0 exp(−β˜nr∫0tD^(t′)dt′)+β˜nrS^nr,0 exp(−β˜nr∫0tD^(t′)dt′)]−κwμII^ (2.11) Using (2.7d) to recast the final term of (2.11) in terms of D^ results in   d2D^dt¯2(μI+1ν)dD^dt¯+μIνD^=κwD^[β˜rS^r,0 exp(−β˜nr∫0tD^(t′)dt′)+β˜nrS^nr,0 exp(−β˜nr∫0tD^(t′)dt′)] (2.12) Because D^ is initiated by the fractional infected population (2.7d) and (2.7h), which is typically small, the exponents in (2.12) will initially be small, and the exponentials can be approximated by 1 for small times. Additionally, if conditions can be found where the solution for D^ exhibits long-term exponential decay, approximating the exponentials by 1 may be valid for longer time periods. This approximation will be checked a posteriori, as well as by full numerical solution to the full equation (2.12). For now, we set the exponentials to 1.0 and obtain the equation   d2D^dt¯2+(μI+1ν)dD^dt¯+μIνD^−κwβ˜wD^=0 (2.13a) where   β˜w=[β˜rS^r,0+β˜nrS^nr,0]= β˜nrS^0[Tinfs+(1−fs)] (2.13b) Assuming solutions to (2.13a) of the form   D^(t)=exp[Mt] (2.14a) yields   M1,2=−(μI+1/ν)±(μI+1/ν)2−4(μI/ν−κwβ˜w)2 (2.14b) or   M1,2=−(μI+1/ν)±(μI−1/ν)2+4κwβ˜w2. (2.14c) A solution to (2.13a) satisfying (2.7l) is   D^(t)=C{exp[M1t]−exp[M2t]}, (2.15) where C is a constant to be determined. The solution for I^ can be derived directly from (2.7d). After performing the operations in (2.7d) and using (2.7k) to prescribe C, we have   I^(t)=I^0M1−M2{M1exp[M1t]−M2exp[M2t]+(1/ν)exp[M1t]−(1/ν)exp[M2t]} (2.16a) or   I^(t)=I^0(μI−1/ν)2+4κwβ˜w{M1exp[M1t]−M2exp[M2t]+1νexp[M1t]−1νexp[M2t]} . (2.16b) The expression (2.16b) for the infected population will have exponentially decaying solutions if both M1 and M2 are negative. From (2.14c), this will occur if   (μI+1ν)>(μI−1/ν)2+4κwβ˜w (2.17a) or   μIν>κwβ˜w. (2.17b) In terms of original variables, (2.17b) can be written (using (2.6b) and (2.13b))   [Tout fi+(1−fi)][Tin fs+(1−fs)]<μI(1/v)S^0β˜nrκnr. (2.18) This threshold equation, which provides the criterion for determining whether the infected population will grow with time, is analysed in Section 3. The assumption that the exponents in the full governing equation (2.12) are small (allowing (2.13) to apply) can be checked using the solution (2.15) for the normalized droplet number. Using (2.15) in the integral   −β˜r∫0tD^(t′)dt′ (2.19a) occurring in the exponents in (2.12), with   C=κwI^0M1−M2 (2.19b) derived from (2.16a), we obtain   −β˜r∫0tD^(t′)dt′=−β˜rκwI^0M1−M2[exp[M1t]−1M1−exp[M2t]−1M2] (2.19c) Since M1 and M2 are both negative, the exponential terms are small for times greater than about 1/(μI+1ν) (2.14c). For the conditions considered by Stilianakis & Drossinos (2010), this time scale is on the order of 10 days. Beyond this time, the approximation   −β˜r∫0tD^(t′)dt′≈−β˜rκwI^0M1−M2[M2−M1M2M1]=β˜rκwI^0M1M2 (2.20) can be made. The product M1M2 is given by (2.14c)   4μI⋅1ν−4κwβ˜w (2.21a) and thus   −β˜r∫0tD^(t′)dt′≈β˜rκwI^04(μI⋅1ν−κwβ˜w) (2.21b) The denominator is assumed positive from (2.17b). The entire expression will be small if   μI⋅1ν−κwβ˜wβ˜rκw≫ I^0/4 (2.21c) This result is analysed further in Section 4. 2.3 Governing equations for a polydisperse aerosol While the monodisperse aerosol model provides valuable insight into the role of RPDs in reducing the spread of infection, quantitative predictions should be based upon a polydisperse model, as most pathogen size distributions span the range from a few nanometers to a few micrometers. In the polydisperse model, a more realistic aerosol distribution consisting of multiple sizes is considered, and the size dependence of the parameters is incorporated. For example, the deposition rate within the lungs, the settling rate of the droplet, and the level of protection provided by the RPD are all a function of droplet size. In the case of a polydisperse aerosol, we follow Stilianakis & Drossinos (2010) and identify M discrete droplet size bins. The total number of droplets varies for each size bin; the number in the jth bin is designated Dj(t). The transmission rate β˜ is a function of droplet diameter, because the number of pathogens per droplet Npand the deposition probability q are functions of the aerosol size, as are the fractions of pathogens Tin and Tout transmitted by the RPDs. The production rate κ and droplet removal rate 1/ ν also depend upon droplet size. Writing the model equations, given by (2.7a)–(2.7d) for a single size, for all of the size bins yields the following system:   ddtS^r=−(β˜r,1D^1+β˜r,2D^2+⋯)S^r, (2.22a)  ddtS^nr=−(β˜nr,1D^1+β˜nr,2D^2+⋯)S^nr, (2.22b)  ddtI^=−ddt(S^r+S^nr)−μII^, (2.22c)  ddtD^1=κw,1I^−1ν1D^1,ddtD^2=κw,2I^−1ν2D^2,… (2.22d) where   β˜r,j=Tin,jcBVclτpqiNp,j (2.22e)  β˜nr,j=cBVclτpqjNp,j (2.22f)  κw,j=κr,jfi+κnr,j(1−fi) (2.22g)  vj−1=(1+cτ)BVclqj+μd+θj (2.22h) We consider boundary conditions analogous to those for the monodisperse case:   S^r(0)=S^r0, (2.22i)  S^nr(0)=S^nr,0, (2.22j)  I^(0)=I^0 (2.22k)  D^j(0)=0,j=1,2,⋯,M. (2.22l) 2.4 Approximate solution for a polydisperse aerosol We next determine whether (2.22a), like (2.7) for monodisperse droplet distribution, admit solutions where the number of infected individuals is exponentially decaying due to the presence of RPDs. To that end, we multiply the governing equation for D^j (2.22d) by β˜r,j, and add all the equations:   ddt(β˜r,1D^1+ β˜r,2D^2+⋯)=(κw,1β˜r,1+κw,2β˜r,2+⋯)I^−(β˜r,1ν1D^1+β˜r,2ν2D^2+⋯) (2.23a) An identical operation using the parameter β˜nr,j yields   ddt( β˜nr,1D^1+ β˜nr,2D^2+⋯)=(κw,1β˜nr,1+κw,2β˜nr,2+⋯)I^−(β˜nr,1ν1D^1+β˜nr,2ν2D^2+⋯) (2.23b) We now assume that the droplet removal rate 1/νj is the same for all size bins. This can be an unrealistic assumption, particularly because the gravitational settling rate θis different for different droplet diameters. Several strategies can be employed to overcome this shortcoming. The first is to use in the computations the value of 1/νj yielding the most stringent requirements on the RPDs necessary to avoid an epidemic. This conservative value would be the smallest of the 1/νj values, since 1/νj characterizes droplet removal. Another strategy for making use of a model with identical 1/νj values in all size bins is to use the smallest and largest values of 1/νj in the calculations, in order to bracket the actual infected rate. Assuming 1/νj=1/ν for all j, and setting   D˙r=β˜r,1D^1+β˜r,2D^2+⋯ (2.24a)  D˙nr=β˜nr,1D^1+β˜nr,2D^2+⋯ (2.24b) we obtain the following modification of (2.22a)–(2.22d) and (2.23a) and (2.23b):   ddtS^r=−D˙rS^r (2.25a)  ddtS^nr=−D˙nrS^nr (2.25b)  ddtI^=−ddt(S^r+S^nr)−μII^ (2.25c)  ddtD˙r=δrI^−1νD˙r (2.25d)  ddtD˙nr=δnrI^−1νD˙nr (2.25e) with   δr=κw,1β˜r,1+κw,2β˜r,2+⋯ (2.25f)  δnr=κw,1β˜nr,1+κw,2β˜nr,2+⋯ (2.25g) The relevant boundary conditions are:   S^r(0)=S^r0, (2.25h)  S^nr(0)=S^nr,0 (2.25i)  I^(0)=I^0 (2.25j)  D˙r(0)=0, (2.25k)  D˙nr(0)=0. (2.25l) Equations (2.25a) and (2.25b) can be integrated into forms similar to (2.9) and (2.10). Using those results in (2.25d) and incorporating this modified form of (2.25d) into the time derivative of (2.25e) yields (see steps leading to (2.12)):   d2D˙rdt2+(μI+1ν)dD˙rdt+μIνD˙r=δr[D˙rS^r,0exp(−∫0tD˙r(t′)dt′)+D˙nrS^nr,0exp(−∫0tD˙nr(t′)dt′)] (2.26a)  d2D˙nrdt2+(μI+1ν)dD˙nrdt+μIνD˙nr=δnr[D˙rS^r,0exp(−∫0tD˙r(t′)dt′)+D˙nrS^nr,0exp(−∫0tD˙nr(t′)dt′)]. (2.26b) As with the monodisperse case, we assume that the exponents in (2.26a) and (2.26b) are small and the exponential can be approximated by 1.0. The condition for the validity of this assumption is discussed below. Here we set the exponentials to 1.0 and obtain   d2D˙rdt2+(μI+1ν)dD˙rdt+μIνD˙r−δr[D˙rS^r,0+D˙nrS^nr,0]=0, (2.27a)  d2D˙nrdt2+(μI+1ν)dD˙nrdt+μIνD˙nr−δnr[D˙rS^r,0+D˙nrS^nr,0]=0. (2.27b) Equations (2.27a) and (2.27b) are coupled, precluding a simple solution similar to (2.15). However, a relatively simple solution can be obtained for three important special cases: (i) All susceptible persons initially deploy RPDs, i.e. fs=1. In this case, (2.27a) reduces to   d2D˙rdt2+(μI+1ν)dD˙rdt+μIνD˙r−δrD˙r=0. (2.28a) This equation is similar to (2.13a) and will have exponentially decaying solutions ( D˙r~exp(Mt) ,M<0) when   μIν>δr, (2.28b) i.e.   κw,1β˜r,1+κw,2β˜r,2+⋯<μI⋅1ν (2.28c) Since the simplification was made prior to (2.24a) that the 1/νj values were the same for all bins, in order for (2.28c) to be accurate the smallest (over all bins) 1/νj value, (1/ν)min, should be used. By making this change, and explicitly entering the features of the RPD (analogous to (2.18) in the monodisperse case), we can rewrite (2.28c) as   ∑jκnr,jβ˜nr,jTin,j[Tout,jfi+(1−fi)]<μI(1/v)min (2.28d) When this criterion for exponentially decaying solutions holds, (2.27b) will have solutions of the form D˙nr~texp(Mt) ,M<0, i.e. solutions manifesting some algebraic growth but also exponentially decaying. The number of infected people is (from (2.25d))   I^(t)=I^0(μI−1/ν)2+4δr{M1exp[M1t]−M2exp[M2t]+1νexp[M1t]−1νexp[M2t]} , (2.28e) where   M1,2=−(μI+1/ν)±(μI−1/ν)2+4δr2. (2.28f) (ii) No susceptible individuals initially deploy RPDs, i.e. S^r,0=0,S^nr,0=1 .In this case, (2.27b) reduces to   d2D˙nrdt2+(μI+1ν)dD˙nrdt+μIνD˙nr−δnrD˙nr=0. (2.29a) The criterion for exponentially decaying solutions, analogous to (2.28b), is   μIν>δnr (2.29b) i.e.   κw,1β˜nr,1+κw,2β˜nr,2+⋯<μI⋅1ν (2.29c) Again setting 1/ν equal to the minimum value of 1/νj from all of the bins, and entering the RPD characteristics explicitly, we obtain   ∑jκnr,jβ˜nr,j[Tout,jfi+(1−fi)]<μI(1/v)min. (2.29d) The solution for the number of infected persons is:   I^(t)=I^0(μI−1/ν)2+4δnr{M1exp[M1t]−M2exp[M2t]+1νexp[M1t]−1νexp[M2t]} , (2.29e) where   M1,2=−(μI+1/ν)±}(μI−1/ν)2+4δnr2. (2.29f) (iii) Half of the susceptible population initially deploys RPDs, i.e. S^r,0=1/2,S^nr,0=1/2. Adding (2.27a) and (2.27b) gives,   d2dt2(D˙r+D˙nr)+(μI+1ν)ddt(D˙r+D˙nr)+μIν(D˙r+D˙nr)−12(δr+δnr) (D˙r+D˙nr)=0. (2.30a) The combination field D˙r+D˙nr will have exponentially decaying solutions if   μIν>12(δr+δnr), (2.30b) i.e.   {κw,1(β˜r,1+β˜nr,12)+κw,2(β˜r,2+β˜nr,22)+⋯}<μI⋅1ν (2.30c) or   ∑jκnr,jβ˜nr,j[12Tin,j+12][Tout,jfi+(1−fi)]<μI(1/v)min. (2.30d) The evolution of the infected population is described by:   I^(t)=I^0(μI−1/ν)2+2(δr+δnr){M1exp[M1t]−M2exp[M2t]+1νexp[M1t]−1νexp[M2t]} , (2.30e) where   M1,2=−(μI+1/ν)±}(μI−1/ν)2+2(δr+δnr)2 (2.30f) The criterion for validity of the small-exponent approximation leading to (2.28), derived in the same manner as (2.21c), is   μI⋅1ν−δrδr≫I^0/4 (2.31a) or, writing δr explicitly,   μI•1ν−(κw,1β˜r,1+κw,2β˜r,2+⋯)κw,1β˜r,1+κw,2β˜r,2+⋯≫I^0/4. (2.31b) Similar criteria involving δnr and (δr+δnr)/2 apply for (2.29) and (2.30), respectively. 3. Results The threshold inequality, (2.18), is the essential result for the case of a monodisperse aerosol. For convenience in describing the results we recast it as:   Ωm=[Toutfi+(1−fi)][Tinfs+(1−fs)]{S^0β˜nrκnrμI(1/v)}<1, (3.1) where Ω is the ‘reproduction number’ (Stilianakis & Drossinos, 2010; Diekmann et al., 1990) and the subscript ‘m’ denotes monodisperse. The term is curly brackets is the reproduction number in the absence of protective equipment. Roughly speaking, it is the ratio of the mechanisms promoting infection to those preventing it. For high values of this number, conditions are favourable for the growth of infection. To prevent the growth of infection, the characteristics of the protective strategy, given by the two factors in square brackets, must be such that the overall reproduction number is less than 1. Examining the factors in square brackets, we see that the prevention strategy involves a high degree of compliance in deploying RPDs by the susceptible population ( fs≈1) and a low-inward-transmission-rate barrier for the susceptible population ( Tin≈0), or a high degree of compliance in deploying RPDs by the infected population ( fi≈1) combined with a low-outward-transmission-rate barrier for the infected population ( Tout≈0), or both. The threshold inequality (3.1) was analysed numerically in the computations described below. For the polydisperse case, generalized reproduction numbers and threshold inequalities, corresponding to (2.28d)–(2.30d), can be defined by   Ωp=∑jκnr,jβ˜nr,jTin,j[Tout,jfi+(1−fi)]μI(1/v)min<1(fs=1) , (3.2a)  Ωp=∑jκnr,jβ˜nr,j[Tout,jfi+(1−fi)]μI(1/v)min<1(fs=0) , (3.2b)  Ωp=∑jκnr,jβ˜nr,j[12Tin,j+12][Tout,jfi+(1−fi)]μI(1/v)min<1(fs=12) . (3.2c) These threshold inequalities were also analysed in the computations. Numerical solutions to the governing equations (2.7a)–(2.7h) for the monodisperse case, and (2.25a)–(2.25l) for the polydisperse case, were obtained to validate the approximate analytic solution. The differential equations were solved using a variable time-step Runge–Kutta method, as implemented in the Matlab (MathWorks, Inc., Natick, MA) mathematical software. The analytical approximations for the infected population are given by (2.16b) in the monodisperse case and (2.28d)–(2.30d) for a polydisperse aerosol. The threshold criteria involving the reproduction number (3.1) for the monodisperse case and (3.2a)–(3.2c) for polydisperse were used to prescribe the quality of protection required to avoid a growing infected population, under various conditions. Unless otherwise stated, the parameters used in the calculations are those provided by Stilianakis & Drossinos (2010) for an influenza epidemic. The values are provided in Table 1. For the polydisperse case, we considered two bins, one having a characteristic droplet diameter of 4 microns, and the other a diameter of 8 microns. In all computations, the initial infected population was 0.1% of the total population. Parameters characterizing the protective equipment were varied in the calculations. Table 1. Parameter values used in calculations Parameters  Values  Breathing rate, B ( m3/day)  24  Volume of the personal cloud of an infected person, Vcl ( m3)  8  Characteristic breathing time, τ(days)  0.0139  Probability of infection by an inhaled pathogen, p  0.052  Droplet diameters, d1 and d2(μm)  4 and 8  Inhaled-droplet deposition probability, q1 and q2  0.88 and 0.99  Number of pathogens per droplet, Np1 and Np2  9.15×10−4 and 8.2×10−3  c (per day)  13  Inactivation rate of airborne pathogens, μd (per day)  8.64  Gravitational settling rate, θ1 and θ2 (per day)  28.8 and 113.2  Infection recovery rate, μI (per day)  0.2  Droplet production rate for two particle sizes without respirator, κnr,1κnr,2  4.1×105, 1.92×104  Fraction of initial number of infecteds, I0  0.001  Fraction of incoming pathogens transmitted by the RPD, Tin  0.25  Outward transmission rate, Tout  0.1  Parameters  Values  Breathing rate, B ( m3/day)  24  Volume of the personal cloud of an infected person, Vcl ( m3)  8  Characteristic breathing time, τ(days)  0.0139  Probability of infection by an inhaled pathogen, p  0.052  Droplet diameters, d1 and d2(μm)  4 and 8  Inhaled-droplet deposition probability, q1 and q2  0.88 and 0.99  Number of pathogens per droplet, Np1 and Np2  9.15×10−4 and 8.2×10−3  c (per day)  13  Inactivation rate of airborne pathogens, μd (per day)  8.64  Gravitational settling rate, θ1 and θ2 (per day)  28.8 and 113.2  Infection recovery rate, μI (per day)  0.2  Droplet production rate for two particle sizes without respirator, κnr,1κnr,2  4.1×105, 1.92×104  Fraction of initial number of infecteds, I0  0.001  Fraction of incoming pathogens transmitted by the RPD, Tin  0.25  Outward transmission rate, Tout  0.1  View Large 3.1 Numerical results for monodisperse case In Fig. 2, the size infected population is plotted as a function of time, for a monodisperse aerosol. Different values of the reproduction number (3.1), generated by considering different values of the droplet production rate κnr, are considered. In the calculations it was assumed that 50% of the susceptible population deploys RPDs ( fs=50%), and that the RPDs allow 50% of the aerosols incoming into the RPD during breathing to be transmitted ( Tin=50%). A Tin value of 50% is representative of a respirator that has not been fit-tested (Coffey et al. 2004). The fraction of infected persons deploying RPDs was 10% ( fi=10%) and the outward transmission rate was 10% ( Tout=10%, Patel et al., 2015). Decaying infected populations can be seen for the reproduction numbers less than 1, while growing infected populations (at least initially) can be seen for reproduction numbers greater than 1. For reproduction numbers less than 1, the analytical approximation closely matches the numerical one for the entire 10-day time period. For Ωm=3, the analytical approximation also closely approximates the numerical one for the 10-day duration, even though the infected field is growing. For the higher reproduction number ( Ωm=6), the analytical approximation tracks the numerical one closely for about half the 10-day period, but then fails to predict the turnover in the infected-population count. Fig. 2. View largeDownload slide Comparison of the analytical and numerical solutions for the infected population, in situations where the monodisperse reproduction number is less than 1 (bottom curves) and greater than 1 (upper curves). Fig. 2. View largeDownload slide Comparison of the analytical and numerical solutions for the infected population, in situations where the monodisperse reproduction number is less than 1 (bottom curves) and greater than 1 (upper curves). A subsequent set of calculations was performed in which it was assumed that 50% of the susceptible population deploys RPDs ( fs=50%), and that the fraction of infected individuals deploying RPDs was zero ( fi=0). It was determined that for an inward transmission rate of slightly greater than 16%, the threshold inequality was just satisfied ( Ωm=1) and the infected population should not grow. In other words, as long as the RPDs deployed by the susceptible population were able to contain 84% of the pathogens, even if the infected individuals were not protected, a pandemic would not spread. This prediction was checked by full numerical solution, using Tin=16% as well as values slightly lower and higher, namely 15, 17 and 18%. The results are shown in Fig. 3. The dynamics of the infected population follow the prediction of the threshold inequality, with Tin values of 15 and 16% leading to decay, and 17 and 18% leading to initial growth. The initial growth rate, as well as the turn-over time, increased significantly with increasing barrier transmission rate. Fig. 3. View largeDownload slide Infected population as a function of time, when the susceptible population deploys RPDs having different inward transmission rates. Curves were computed from full numerical model for a monodisperse droplet. Analytical threshold inequality predicts that the transition between a growing infected population and a decaying one occurs for Tin slightly above 16. Fig. 3. View largeDownload slide Infected population as a function of time, when the susceptible population deploys RPDs having different inward transmission rates. Curves were computed from full numerical model for a monodisperse droplet. Analytical threshold inequality predicts that the transition between a growing infected population and a decaying one occurs for Tin slightly above 16. The threshold inequality (3.1) was used to identify the level of integrity required of the barriers worn by the susceptible population, in order to avoid a growing infection rate, as a function of the level of compliance in wearing RPDs by the infected population. The outward transmission rate was taken to be Tout=10%. Different levels of compliance by the susceptible population ( fs=0.1, 0.3, 0.5, 0.7 and 0.9) were considered. The results are shown in Fig. 4. As expected, when the compliance rate of the infected population increases (higher fi), lower-integrity barriers (higher Tinvalues) worn by the susceptible population can still prevent growth in the infection population. When the compliance of the infected population exceeds approximately 0.14 (for the parameter values considered), the number of pathogens produced is low enough that infection growth does not occur regardless of the behaviour of the susceptible population. That is, for any level of compliance by the susceptible population (all fs values), and the use of barriers that transmit all pathogens ( Tin=100%), infection growth does not occur. This is manifested in Fig. 4 in the convergence of all the curves to the point ( fi=0.14, Tin=100%). For a compliance rate by the susceptible group of 10% ( fs=0.1), the curve originating at ( fi=0.03, Tin=0) implies that for compliance rates by the infected group of less than 3%, a growth in the infected population cannot be prevented except by perfect protection ( Tin=0) of the susceptible population. Fig. 4. View largeDownload slide Level of barrier integrity (designated by the transmission rate Tin) required in order to avoid a growing infection rate, as a function of the fraction of infected persons deploying RPDs. A monodisperse droplet distribution was assumed. Curves are plotted for different fractions of the susceptible population initially wearing RPDs. Outward transmission rate for the infected population is 10% ( Tout=0.1). Fig. 4. View largeDownload slide Level of barrier integrity (designated by the transmission rate Tin) required in order to avoid a growing infection rate, as a function of the fraction of infected persons deploying RPDs. A monodisperse droplet distribution was assumed. Curves are plotted for different fractions of the susceptible population initially wearing RPDs. Outward transmission rate for the infected population is 10% ( Tout=0.1). 3.2 Numerical results for polydisperse case Figure 5 for the polydisperse case is analogous to Fig. 2 for the monodisperse. The barrier parameters were the same ( fi=10%, Tout=10%, fs=50%, Tin=50%) as in Fig. 2, and again the reproduction number Ωp was varied by adjusting the droplet production rate κnr The infected population decreases as a function of time for reproduction numbers less than 1 (bottom two sets of curves in Fig. 5), and increases for reproduction numbers greater than 1 (top two sets of curves). As with the monodisperse droplets, the analytical approximation matches the numerical one closely for the decaying infected populations. For the increasing populations, the analytical approximation matches the numerical one well for about half of the 10-day duration, but does not predict the maximum and eventual decrease in the population. Fig. 5. View largeDownload slide Comparison of the analytical and numerical solutions for the infected population, in situations where the polydisperse reproduction number is less than 1 (bottom curves) and greater than 1 (upper curves). Fig. 5. View largeDownload slide Comparison of the analytical and numerical solutions for the infected population, in situations where the polydisperse reproduction number is less than 1 (bottom curves) and greater than 1 (upper curves). Figure 6 is the polydisperse equivalent to Fig. 4 for the monodisperse case. As indicated above, for the polydisperse distribution, a second size bin (8-micron) was added to the first bin (4-micron). There are thus additional inhalable droplets for the polydisperse case, making the challenge against RPDs deployed by the susceptible population greater. The allowable transmission rate ( Tin) to prevent the growth of infection is therefore lower for the polydisperse case. When no susceptibles deploy RPDs ( fs=0), the transmission rate is immaterial, and hence the result is a vertical line. The fi value (roughly 0.45) at which the vertical line is plotted is the threshold value for compliance by the infected population (assuming a Tout of 0.1), above which there will be decaying infected populations. Fig. 6. View largeDownload slide Level of barrier integrity (designated by the transmission rate Tin) that is required in order to avoid a growing infection rate, as a function of the fraction of infected people deploying RPDs, for a polydisperse distribution. Curves are plotted for different fractions of the susceptible population initially wearing RPDs. Outward transmission rate for the infected population is 10% ( Tout=0.1). Fig. 6. View largeDownload slide Level of barrier integrity (designated by the transmission rate Tin) that is required in order to avoid a growing infection rate, as a function of the fraction of infected people deploying RPDs, for a polydisperse distribution. Curves are plotted for different fractions of the susceptible population initially wearing RPDs. Outward transmission rate for the infected population is 10% ( Tout=0.1). 4. Discussion This work extends the SIR model of Stilianakis & Drossinos (2010) for predicting the spread of infectious disease by inhalable droplets to include effect of protective measures. The most important use made of the extended equations was to derive an analytical threshold inequality for predicting the level of protection required from RPDs in order to avoid a growing infection rate. This inequality (3.1) for the monodisperse case and (3.2a,b,c) for polydisperse aerosols, is written in terms of a generalized reproduction number containing the properties of the RPDs and the level of compliance of the susceptible and infected populations. The predictions of threshold for growth from the analytical model agreed well with full numerical solutions to the SIR equations. Another useful result emerging from the analytical approximation of the SIR equations is the ability of the approximate model to track the initial dynamics of the infected population (Figs 2 and 5), even in the case of reproduction numbers substantially bigger than 1. This ability could prove useful, e.g. in identifying the time interval available for medical assistance to arrive before the infected population reaches a certain size, as a function of the properties of the pathogen, the level of protection, and the characteristics of the susceptible population. An initial normalized infected population of I^0=0.001 was used in our calculations. To the extent that the linearized equation (2.13a) is valid, the results (e.g. threshold equation) are independent of the initial infected population (solutions for different initial conditions are the same to within a scaling factor). For decaying solutions, the criterion for validity of the small-exponent approximation, and hence of the linear equation (2.13a), was shown to be (2.21c). For the parameters considered in the calculations (including Tin=Tout=50%, fs=fi=50%), this condition is roughly I^0≪1. For growing solutions, the analysis leading to (2.21c) is not applicable, as the exponential terms grow large with time. However, for sufficiently small times the small-exponent approximation is still valid with growing solutions. If we consider the exponent, provided in (2.19c), and take the largest of the terms exp[ M1] and exp[ M2], we can say that the exponent will be small if   β˜rκwI^0M2−M1exp[M1t]M1≪1 (4.1a) or   t≪1M1ln[M1(M2−M1)β˜rκwI^0] . (4.1b) For the parameters used in the calculations, this inequality becomes, roughly, t≪30 days. This criterion is consistent with the results of Figs 2 and 5, where the analytical approximation correctly predicted the infection dynamics for the growing population over a period of at least a few days. As noted above, this period is useful for identifying the initial growth rate of an epidemic, and devising an initial intervention strategy. The symmetry in (3.1) between the parameters ( Tin,fs) characterizing the level of protection for the susceptible population and the parameters ( Tout,fi) characterizing the level of protection for the infected population indicate that, formally, deployment of protective equipment by one population is as effective in reducing the likelihood of an epidemic as deployment by the other. However, as reflected in the parameter selection in the computations ( Tout=10%, Tin=50%) protective equipment may be more effective in reducing the outward flux of pathogens than the inward flux (Coffey et al., 2004; Patel et al., 2015). In such a case, preferential attention to the infected population (e.g. to bring about a higher level of compliance) may be a worthwhile strategy. In the derivation of the analytical approximation for the polydisperse aerosol, including the threshold equations (3.2a)–(3.2c), it was assumed that the droplet removal rate 1/ν was the same for all droplet sizes. This assumption led to a set of equations, (2.25), that was reduced in size and more mathematically tractable than the full set of polydisperse equations, (2.22). Choosing the smallest removal rate over all size bins as the one removal rate, as we did in our calculations, results in a conservative estimate of the infection growth rate. That is, the number of infected persons determined through a numerical or analytical approximation to (2.25) will be greater than or equal to the number of infected persons determined by solution to (2.22). For the bimodal set of data used in our polydisperse calculations, the difference between assuming a uniform removal rate and using different removal rate for the two size bins can be gleaned from Fig. 7. The conditions match those of Fig. 5, for a reproduction number of 4. The predicted number of infected individuals is substantially lower (about a factor of 4) for the variable removal rate, owing to the very large settling rate for the 8-micron size. In practice, if highly conservative values are to be avoided in favour of more accurate predictions, average droplet removal rate can be used instead of the smallest, or calculations using the lowest and highest removal rates can be performed to establish a range of infected-population values. In any event, we highlight the desirable property that the analytical approximation leads to a conservative estimate of the infected population. Fig. 7. View largeDownload slide Infected population as a function of time determined three ways: (1) numerical solutions the full polydisperse equations (25); (2) numerical solution of the polydisperse equations assuming uniform droplet removal rate (22) and (3) approximate analytical solution to the polydisperse equations for uniform droplet removal rate (22). Reproduction number is equal to 4. Fig. 7. View largeDownload slide Infected population as a function of time determined three ways: (1) numerical solutions the full polydisperse equations (25); (2) numerical solution of the polydisperse equations assuming uniform droplet removal rate (22) and (3) approximate analytical solution to the polydisperse equations for uniform droplet removal rate (22). Reproduction number is equal to 4. The majority of the effort in this study was devoted to the development of the equations and approximations for the polydisperse case. The value of the polydisperse formulation appeared in the example just given where a significant (factor of 4) difference in the infected population relative to the monodisperse model was found. Additionally, the polydisperse formulation was developed in order to handle more complicated distributions than the present bimodal one. Nicas et al. (2005) list more than 10 size bins in the 1−200 μm range where the number of cough droplets produced in each bin varies by more than 10% from the previous bin. There are 4 bins containing more than 50% of the maximum number of droplets. Even if the size dependence of some parameters is unknown and an average (effectively monodisperse) value must be used, the ability to leverage the known size dependence of important quantities such as the droplet production rate in a model provides more semblance of reality to the calculations. Many of the calculations performed were limited to a short time (e.g. 10 days in Figs 2, 5, 7) relative to the actual duration of a pandemic in which the protective equipment might be deployed. The duration was kept short in order to focus on the applicability of the analytical model in the challenging case of a reproduction number significantly larger than 1. There is no time limit at which the accuracy of the full model ((2.7) for the monodisperse case and (2.22) for polydisperse), or the approximate model for reproduction number less than 1, degrades. The calculations featured in Fig. 3, e.g. cover the full range of infection dynamics: growth, peak and decay. These aspects of disease transmission occur over a period of 100 days. Ideally, the results generated by the present model should be validated in experiments involving protective devices. It should be recognized that validating even standard SIR models, which involves deliberate infection of a susceptible population, is extremely difficult. Experiments involving controlled use of protective equipment are even more complicated. The feasibility of attaching protective devices to animals is likewise low. We do note that the SIR equations were modified in a very systematic manner to account for the effect of RPDs, by adjusting the flux of droplets input to and output from of the respiratory system. Thus, other methods devised to test SIR models that can also modify the inflow or outflow of pathogens, e.g. with pharmacological interventions or by controlling the air quality in an environment, can be used to indirectly validate the present model. In the absence of rigorous validation, the accuracy of the model's absolute predictions for the number of persons infected is not known. However, the accuracy of relative predictions is likely to be more reliable, and such relative predictions represent an important use for the model. In the purchasing and stockpiling of protective equipment for a possible spread of infection, example, a valuable piece of information is an estimate of the relative risk reduction achieved when one possible barrier is replaced by another, more expensive one with a lower transmission rate. Another useful estimate that the model can provide is the degree of compliance by the target population that would be required for an RPD to produce a decay (rather than growth) in the infected population. In such cases, sociological data may suggest that the level of compliance is not achievable, and resources would not be well spent pursuing a strategy involving the given RPD for the given population. The computations in this paper were performed primarily to validate the analytical approximation and explore the dynamics of the infected population near the threshold for infection growth. Much remains to be done using the generalized equations derived in this paper in order to assess the impact of protective measures during an epidemic (i.e. well above threshold). The full differential equations should be numerically solved with the relevant parameters informed by data from actual barrier tests (e.g. Guha et al., 2016). The uncertainty in the predictions must also be quantified. In this regard, the analytical approximation (e.g. (2.16)) is useful. The variability in a prediction, say for the infected field I, with respect to a given parameter p may be written   δI=∂I∂pδp (4.2) where δp is the uncertainty in the parameter p. The sensitivity ∂I∂p may be obtained through direct differentiation of the analytical approximation. Calculations of the infected and susceptible populations under epidemic conditions, as well as the uncertainties in the predictions, are currently underway. Funding This work was supported by the FDA Medical Countermeasures Initiative, Grant Numbers MCM2D20511 and MCM2J270HT. References Chao C. Y. H. Wan M. P. Morawska L. Johnson G. R. Ristovski Z. D. Hargreaves M. Mengersen K. Corbett S. Li Y. Xie X. Katoshevski D. ( 2009) Characterization of expiration air jets and droplet size distributions immediately at the mouth opening. J. Aerosol. Sci. , 40, 122– 133. Google Scholar CrossRef Search ADS   Chen S.-C. Liao C.-M. ( 2008) Modelling control measures to reduce the impact of pandemic influenza among schoolchildren. Epidemiol. Infect. , 136, 1035– 1045. Google Scholar CrossRef Search ADS PubMed  Coburn B. J. Wagner B. G. Blower S. ( 2009) Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1). BMC Med.,  7, 30. Google Scholar CrossRef Search ADS   Coffey C. C. Lawrence R. B. Campbell D. L. Zhuang Z. Calvert C. A. Jensen P. A. ( 2004) Fitting characteristics of eighteen N95 filtering-facepiece respirators. J. Occup. Environ. Hyg. , 1, 262– 271. Google Scholar CrossRef Search ADS PubMed  Committee on the Certification of Personal Protective Technologies. ( 2011) Certifying Personal Protective Technologies: Improving Worker Safety  Washington DC: Institute of Medicine, National Academy of Sciences. Diekmann O. Heesterbeek J. A. P. Metz J. A. J. ( 1990) On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J. Math. Biol. , 28, 365– 382. Google Scholar CrossRef Search ADS PubMed  Furuya H. ( 2007) Risk of transmission of airborne infection during train commute based on mathematical model. Environ. Health. Prev. Med. , 12, 78– 83. Google Scholar CrossRef Search ADS PubMed  Guha S. McCaffrey B. Hariharan P. Myers M. ( 2016) Quantifying the leakage through gaps in surgical masks and pediatric face masks. J. Occup. Hyg. (to appear). Huang S.-H. Chen C.-W. Chang C.-P. Lai C.-Y. Chen C.-C. ( 2007) Penetration of 4.5 nm to 10 μm aerosol particles through fibrous filters J. Aerosol Sci. , 38 719– 727. Google Scholar CrossRef Search ADS   Lee S. A. Grinshpun S. A. Adhikari A. Li W. McKay R. Maynard A. Reponen T. ( 2005) Laboratory and field evaluation of a new personal sampling system for assessing the protection provided by the N95 filtering facepiece respirators against particles. Ann. Occup. Hyg. , 49 245– 257. Google Scholar PubMed  Mniszewski S. M. Del Valle S. Y. Priedhorsky R. Hyman J. M. Hickman K. S. ( 2014) Understanding the impact of face mask usage through epidemic similation of large social networks. Theories and Simulations of Complex Social System (D. Vahid and M. Vijay Kumar eds), 97–115. Nicas M. Nazaroff W. W. Hubbard A. ( 2005) Towards understanding the risk of secondary airborne infection: emission of respirable pathogens. J. Occup. Environ. Hyg. , 2, 143– 154. Google Scholar CrossRef Search ADS PubMed  Patel R. B. Skaria S. D. Mansour M. M. Smaldone G. C. ( 2015) Respiratory source control using a surgical mask: an in vitro study. J Occup. Environ. Hyg. , 15, 1– 29. Stilianakis N. I. Drossinos Y. ( 2010) Dynamics of infectious disease transmission by inhalable respiratory droplets. J. R. Soc. Interface , 7, 1355– 1366. Google Scholar CrossRef Search ADS PubMed  Sze To G. N. Chao C. Y. H. ( 2010) Review and comparison between the Wells–Riley and dose-response approaches to risk assessment of infectious respiratory diseases. Indoor Air , 20, 2– 16. Google Scholar CrossRef Search ADS PubMed  Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications 2016. This work is written by (a) US Government employee(s) and is in the public domain in the US. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical Medicine And Biology: A Journal Of The Ima Oxford University Press

A mathematical model for assessing the effectiveness of protective devices in reducing risk of infection by inhalable droplets

Loading next page...
 
/lp/ou_press/a-mathematical-model-for-assessing-the-effectiveness-of-protective-51Rn0x8rmj
Publisher
Institute of Mathematics and its Applications
Copyright
Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications 2016.
ISSN
1477-8599
eISSN
1477-8602
D.O.I.
10.1093/imammb/dqw018
Publisher site
See Article on Publisher Site

Abstract

Abstract Respiratory protective devices (RPDs) are critical for reducing the spread of infection via inhalable droplets. In determining the type of RPD to deploy, it is important to know the reduction in the infection rate that the RPD enables for the given pathogen and population. This paper extends a previously developed susceptible-infected-recovered (SIR) epidemic model to analyse the effect of a protection strategy. An approximate solution to the modified SIR equations, which compares well with a full numerical solution to the equations, was used to derive a simple threshold equation for predicting when growth of the infected population will occur for a given protection strategy. The threshold equation is cast in terms of a generalized reproduction number, which contains the characteristics of the RPDs deployed by the susceptible and infected populations, as well as the degree of compliance in wearing the equipment by both populations. An example calculation showed that with 50% of the susceptible population deploying RPDs that transmit 15% of pathogens, and an unprotected infected population, an otherwise growing infection rate can be converted to one that decays. When the infected population deploys RPDs, the transmission rate for the RPDs worn by the susceptible population can be higher. 1. Introduction Respiratory protective devices (RPDs) such as respirators constitute an important element of the control strategy against the transmission of infectious diseases. In determining whether to deploy a given RPD to protect a specific population against a certain pathogen, including whether to allocate substantial resources to stockpile RPDs against a possible future threat, realistic estimates of the reduction of risk associated with use of the RPD in the scenario of interest is important. In evaluating the reduction of risk associated with deployment of a given RPD, important parts of the assessment include determination of the intrinsic penetration of the device (Committee on the Certification of Personal Protective Technologies 2011; Mniszewski et al., 2014), which is governed by the microscopic properties of the barrier and the amount of leakage due to imperfect fit (Coffey et al., 2004; Lee et al., 2005) between a particular facial profile and the RPD. Another important part of the risk assessment for RPDs is an estimate of the change in the infection rate associated with the change in pathogen transport introduced by the presence of the barrier. The change in pathogen transport could apply to either an uninfected person attempting to reduce the intake of pathogens or an infected person attempting to reduce the output. While numerous models of various types have been developed to simulate the spread of epidemics (Furuya, 2007; Chen & Liao, 2008; Coburn et al., 2009; Sze To & Chao, 2010), to the authors' knowledge no systematic treatment of the effect of RPDs has been incorporated into the models. Of particular interest for RPD evaluation is the transmission of disease by inhalable respiratory droplets. A model that is well suited to estimating the spread of disease by inhalable respiratory droplets is the susceptible-infected-recovered (SIR) epidemic model. In the SIR models, the influence of RPDs can be accounted for in a systematic manner in the parameters within the system of differential equations. Stilianakis & Drossinos (2010) developed an SIR model for disease transmission with inhalable droplets as the transmission vector, assuming a closed, homogeneously mixed population. Epidemics mediated by both airborne and settled droplets were considered, and the time course for the susceptible, infected and recovered populations derived using parameters derived from experimental evidence. Relevant parameters of the model include the droplet production rate, gravitational settling rate, deposition probability and number of pathogens per droplet, all of which are a function of droplet size. Stilianakis & Drossinos (2010) do not address the effect of protective equipment. Of considerable importance in evaluating risk in the presence of protective equipment is the polydisperse nature of the droplets. Even without considering protective equipment, distinguishing between different sizes can be important, as studies have found that cough droplets possess a size distribution ranging often from less than 1 μm to over 100 μm (Chao et al., 2009). Nicas et al. (2005) estimated that the droplet production rate at 4 μm is over 20 times that at 8 μm. When protection is considered, the strong dependence of RPD performance on size must be taken into account. Huang et al. (2007) showed that the penetration of 1 micron particles for some respirators could be an order-of-magnitude higher than the penetration at 5 microns. In this paper we modify the SIR theory as presented by Stilianakis & Drossinos (2010) to account for the influence of RPDs. We account for the fact that only a fraction of the population will likely deploy protective measures by dividing the susceptible population into two groups, one of which deploys RPDs and one which does not. Similarly, a fraction of the infected population utilizes RPDs, thereby reducing the source of pathogens. We first implement our modifications under the idealized assumption of a monodisperse aerosol distribution, and then extend the formulation to the more realistic case of a polydisperse distribution. 2. SIR Model in the presence of protective equipment 2.1 Governing equations for a monodisperse aerosol The original SIR model equations for a monodisperse aerosol, in the absence of protective equipment, are (Stilianakis & Drossinos, 2010):   dSdt=−β˜DSN, (2.1a)  dIdt=−dSdt−μII, (2.1b)  dDdt=κI−1νD, (2.1c) where   β˜=cBVclτpqNp (2.1d) and   v−1=(1+cτ)BVclq+μd+θ. (2.1e) 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, μI the infection recovery rate, κ the droplet production rate, c the contact rate between a susceptible person and an infected person, B the breathing rate, Vcl the volume of the personal cloud of an infected person, τ the characteristic breathing time, p the probability of infection by an inhaled pathogen, q the inhaled-droplet deposition probability, Np the number of pathogens per droplet, μd the inactivation rate of airborne pathogens, and θ the gravitational settling rate. The boundary conditions are   S(0)=S0, (2.1f)  I(0)=I0, (2.1g)  D(0)=0. (2.1h) Equations (2.1) describe the evolution of the susceptible and infected populations, as a function of the relevant parameter groups. The processes described by the evolution equations are illustrated in Fig. 1a and b. The protected population will be introduced in the next section; for now this population can be ignored, as can minor notational differences between the figure and (2.1). (Notation in Fig. 1 corresponds to that used in the RPD model below.) The susceptible population converts to the infected population with a time rate of change proportional to β˜DS. The infected persons convert to the recovered population at a rate proportional to μII.The total rate of change of infected individuals (2.1b) is equal to the rate at which susceptible persons convert to infected minus the rate at which infected persons recover. The total number of droplets, described by (2.1c) and shown in Fig. 1b, increases by one mechanism and decreases by four mechanisms. The generation mechanism is production by infected persons, with a rate proportional to κwI. The attenuation mechanisms are: (1) inhalation of an infected person's own droplets (rate proportional to BqD/Vcl), with the droplet ultimately residing with the infected person generating the droplet; (2) inhalation of someone else's droplet (rate proportional to cτ BqD/Vcl), with the destination of the droplet being potentially any of the three populations; (3) loss of infectivity (rate proportional to μDD), with the droplet remaining in the global environment but not being infectious and (4) droplets settling onto a surface under the influence of gravity (rate proportional to θD). Fig. 1. View largeDownload slide (a) Flux of people between the susceptible, infected and recovered populations, along with parameters characterizing the rate at which people convert from one population to the other. Dashed line symbolizes new features of the dynamics introduced by protective equipment. (b) Schematic of droplet dynamics. D(t) is total number of droplets. Droplet generation mechanism is shown to the left of D(t); annihilation mechanisms are shown to the right. Box on the far left denotes origin of the drops; boxes on the far right denote destinations. Next to the mechanism arrows are parameters characterizing the rate of droplet creation or removal Dashed arrow indicates feature that is affected by protective equipment. Fig. 1. View largeDownload slide (a) Flux of people between the susceptible, infected and recovered populations, along with parameters characterizing the rate at which people convert from one population to the other. Dashed line symbolizes new features of the dynamics introduced by protective equipment. (b) Schematic of droplet dynamics. D(t) is total number of droplets. Droplet generation mechanism is shown to the left of D(t); annihilation mechanisms are shown to the right. Box on the far left denotes origin of the drops; boxes on the far right denote destinations. Next to the mechanism arrows are parameters characterizing the rate of droplet creation or removal Dashed arrow indicates feature that is affected by protective equipment. Disease transmission is initiated by a small number of infected persons ( I0≪S0) that have not yet generated a significant number of droplets at the time the disease progression begins to be tracked ( D(0)=0). At any time during the analysis of the disease progression, the recovered population can be derived from the equations S+I+R=N. That is, the total population is assumed constant. The SIR equations can be modified to account for the presence of RPDs in the following manner. We consider a susceptible population that contains a group Sr deploying RPDs, and a group Snr not deploying RPDs. The two groups comprise the total susceptible population, i.e.   S=Sr+Snr. (2.2) The transmission rate β˜nr for the no-RPD group Snris identical to β˜ in (2.1d). For Sr, the reduction in pathogen transmission introduced by the RPD is accounted for by the introduction of the factor Tin representing the fraction of incoming (to the respiratory system of the susceptible population) pathogens transmitted by the RPD when the susceptible person breathes in. For example, for an N95 respirator capable of filtering out 95% of the pathogens, Tin=5%. The transmission rate in the presence of respiratory protection can be written   β˜r=TincBVclτpqNp=Tinβ˜nr. (2.3) This equation can be thought of as a pathogen transmission rate proportional to a reduced breathing rate TinB or a reduced probability of infection Tinp. While the protected group of susceptible individuals ( Sr) converts to the next phase (infected) at a different rate than the unprotected ( Snr) group, a similar dynamic does not hold true for the infected population. The infected population not deploying RPDs ( Inr) converts to the next phase (recovered) at the same rate as those deploying RPDs ( Ir). The two groups produce pathogens at different rates, however, and the droplet product term can be written   κrIr+κnrInr. (2.4a) Analogously with the reduction in incoming pathogens associated with a barrier having an incoming transmission rate Tin, we can model the reduction in outgoing airborne pathogens introduced (through coughing, sneezing, breathing …) by an infected person deploying an RPD using an outward transmission rate Tout. The effective respirable droplet production rate for an infected person wearing an RPD is   κr=Toutκnr. (2.4b) The infected population changes because susceptible people convert to the infected group, and infected persons convert to the recovered group. To prescribe what proportion of the susceptible population deploys RPDs after becoming part of the infected population, at least two possibilities arise. The first is that the susceptible people continue to wear RPDs in the same proportion after becoming infected as they did before they were infected. A second possibility is that they deploy RPDs in the same proportion as the rest of the infected population. We adopt the second choice, for the following mathematical simplification. Since the Ir and Inr groups convert to the recovered population at the same rate, and the newly converted (from the susceptible population) infected persons mimic the behaviour of previously infected persons, the initial fraction of infected people deploying RPDs maintains throughout time. We denote this initial fraction of the infected population by fi and write the infected population as   I=Ir+Inr=fiI+(1−fi)I. (2.5) The production rate in (2.4) can similarly be written   [κrfi+κnr(1−fi)]I=κwI, (2.6a) where κw represents a ‘weighted’ droplet production rate. Using (2.4b), the weighted droplet production rate can also be expressed as:   κw=κnr[Toutfi+(1−fi)]. (2.6b) By employing the modified variables and parameters in (2.2)–(2.6) in the original SIR equations (2.1a)–2.1h, we can compose the SIR equations for the case where a fraction of the population deploys RPDs. After additionally normalizing S, I and D by N to work in terms of fractional populations, we obtain   ddtS^r=−β˜rD^S^r, (2.7a)  ddtS^nr=−β˜nrD^S^nr, (2.7b)  ddtI^=−ddt(S^r+S^nr)−μII^, (2.7c)  ddtD^=κwI^−1νD^, (2.7d) where   S^r=Sr/N, (2.7e)  S^nr=Snr/N, (2.7f)  I^=I/N, (2.7g)  D^=D/N. (2.7h) The boundary conditions are:   S^r(0)=S^r,0=Sr,0/N, (2.7i)  S^nr(0)=S^nr,0=Snr,0/N, (2.7j)  I^(0)=I^0=I0/N, (2.7k)  D^(0)=0. (2.7l) Here S^r,0 and S^nr,0 are the normalized susceptible populations initially deploying RPDs and not deploying RPDs, respectively. It is important to note that it is not precisely true that S^r,0+S^nr,0=1, since S^r,0 and S^nr,0 are normalized by the total population N, not the total initial susceptible population S^0. (The precise relation is S^r,0+S^nr,0+I^0+R^0=1, where R^0 can assumed to be 0 and I^0 to be small.) For convenience in later expressions, we introduce the fraction of the total initial susceptible population, fs, that deploys RPDs. Then (2.7i) and (2.7j) can also be written   S^r(0)=fsS^0, (2.7m)  S^nr(0)=(1−fs)S^0. (2.7n) The additional features of disease-transmission dynamics introduced by the presence of protective equipment are denoted with the dashed lines in Fig. 1a and b. As shown in Fig. 1a, protective measures introduce a new population, the protected susceptible population. This group converts to the infected population at a slower rate ( β˜r<β˜nr) than the unprotected susceptible population. The dashed lines in Fig. 1b signify that, with the availability of protective equipment, the infected population will generate droplets at a slower rate than in the absence of protective equipment (κw<κnr). 2.2 Approximate solution for a monodisperse aerosol We now analyse (2.7a), with the goal of identifying conditions under which the infected population decays with time due to the presence of the RPDs. We begin by noting that (2.7a) may be rewritten   ddt¯(lnS^r)=−β˜nrD^ (2.8) Integration of (2.8) and application of boundary condition (2.7e) yield   S^r(t)=S^r,0exp−(β˜r∫0tD^(t′)dt′) (2.9) Similarly,   S^nr(t)=S^nr,0exp−(β˜nr∫0tD^(t′)dt′) (2.10) Differentiating (2.7d) with respect to t and using (2.7c) for dI^/dt, with ddt(S^r+S^nr) provided by differentiating (2.9) and (2.10), results in the second-order differential equation   d2D^dt¯2+1νdD^dt¯=κwD^[β˜rS^r,0 exp(−β˜nr∫0tD^(t′)dt′)+β˜nrS^nr,0 exp(−β˜nr∫0tD^(t′)dt′)]−κwμII^ (2.11) Using (2.7d) to recast the final term of (2.11) in terms of D^ results in   d2D^dt¯2(μI+1ν)dD^dt¯+μIνD^=κwD^[β˜rS^r,0 exp(−β˜nr∫0tD^(t′)dt′)+β˜nrS^nr,0 exp(−β˜nr∫0tD^(t′)dt′)] (2.12) Because D^ is initiated by the fractional infected population (2.7d) and (2.7h), which is typically small, the exponents in (2.12) will initially be small, and the exponentials can be approximated by 1 for small times. Additionally, if conditions can be found where the solution for D^ exhibits long-term exponential decay, approximating the exponentials by 1 may be valid for longer time periods. This approximation will be checked a posteriori, as well as by full numerical solution to the full equation (2.12). For now, we set the exponentials to 1.0 and obtain the equation   d2D^dt¯2+(μI+1ν)dD^dt¯+μIνD^−κwβ˜wD^=0 (2.13a) where   β˜w=[β˜rS^r,0+β˜nrS^nr,0]= β˜nrS^0[Tinfs+(1−fs)] (2.13b) Assuming solutions to (2.13a) of the form   D^(t)=exp[Mt] (2.14a) yields   M1,2=−(μI+1/ν)±(μI+1/ν)2−4(μI/ν−κwβ˜w)2 (2.14b) or   M1,2=−(μI+1/ν)±(μI−1/ν)2+4κwβ˜w2. (2.14c) A solution to (2.13a) satisfying (2.7l) is   D^(t)=C{exp[M1t]−exp[M2t]}, (2.15) where C is a constant to be determined. The solution for I^ can be derived directly from (2.7d). After performing the operations in (2.7d) and using (2.7k) to prescribe C, we have   I^(t)=I^0M1−M2{M1exp[M1t]−M2exp[M2t]+(1/ν)exp[M1t]−(1/ν)exp[M2t]} (2.16a) or   I^(t)=I^0(μI−1/ν)2+4κwβ˜w{M1exp[M1t]−M2exp[M2t]+1νexp[M1t]−1νexp[M2t]} . (2.16b) The expression (2.16b) for the infected population will have exponentially decaying solutions if both M1 and M2 are negative. From (2.14c), this will occur if   (μI+1ν)>(μI−1/ν)2+4κwβ˜w (2.17a) or   μIν>κwβ˜w. (2.17b) In terms of original variables, (2.17b) can be written (using (2.6b) and (2.13b))   [Tout fi+(1−fi)][Tin fs+(1−fs)]<μI(1/v)S^0β˜nrκnr. (2.18) This threshold equation, which provides the criterion for determining whether the infected population will grow with time, is analysed in Section 3. The assumption that the exponents in the full governing equation (2.12) are small (allowing (2.13) to apply) can be checked using the solution (2.15) for the normalized droplet number. Using (2.15) in the integral   −β˜r∫0tD^(t′)dt′ (2.19a) occurring in the exponents in (2.12), with   C=κwI^0M1−M2 (2.19b) derived from (2.16a), we obtain   −β˜r∫0tD^(t′)dt′=−β˜rκwI^0M1−M2[exp[M1t]−1M1−exp[M2t]−1M2] (2.19c) Since M1 and M2 are both negative, the exponential terms are small for times greater than about 1/(μI+1ν) (2.14c). For the conditions considered by Stilianakis & Drossinos (2010), this time scale is on the order of 10 days. Beyond this time, the approximation   −β˜r∫0tD^(t′)dt′≈−β˜rκwI^0M1−M2[M2−M1M2M1]=β˜rκwI^0M1M2 (2.20) can be made. The product M1M2 is given by (2.14c)   4μI⋅1ν−4κwβ˜w (2.21a) and thus   −β˜r∫0tD^(t′)dt′≈β˜rκwI^04(μI⋅1ν−κwβ˜w) (2.21b) The denominator is assumed positive from (2.17b). The entire expression will be small if   μI⋅1ν−κwβ˜wβ˜rκw≫ I^0/4 (2.21c) This result is analysed further in Section 4. 2.3 Governing equations for a polydisperse aerosol While the monodisperse aerosol model provides valuable insight into the role of RPDs in reducing the spread of infection, quantitative predictions should be based upon a polydisperse model, as most pathogen size distributions span the range from a few nanometers to a few micrometers. In the polydisperse model, a more realistic aerosol distribution consisting of multiple sizes is considered, and the size dependence of the parameters is incorporated. For example, the deposition rate within the lungs, the settling rate of the droplet, and the level of protection provided by the RPD are all a function of droplet size. In the case of a polydisperse aerosol, we follow Stilianakis & Drossinos (2010) and identify M discrete droplet size bins. The total number of droplets varies for each size bin; the number in the jth bin is designated Dj(t). The transmission rate β˜ is a function of droplet diameter, because the number of pathogens per droplet Npand the deposition probability q are functions of the aerosol size, as are the fractions of pathogens Tin and Tout transmitted by the RPDs. The production rate κ and droplet removal rate 1/ ν also depend upon droplet size. Writing the model equations, given by (2.7a)–(2.7d) for a single size, for all of the size bins yields the following system:   ddtS^r=−(β˜r,1D^1+β˜r,2D^2+⋯)S^r, (2.22a)  ddtS^nr=−(β˜nr,1D^1+β˜nr,2D^2+⋯)S^nr, (2.22b)  ddtI^=−ddt(S^r+S^nr)−μII^, (2.22c)  ddtD^1=κw,1I^−1ν1D^1,ddtD^2=κw,2I^−1ν2D^2,… (2.22d) where   β˜r,j=Tin,jcBVclτpqiNp,j (2.22e)  β˜nr,j=cBVclτpqjNp,j (2.22f)  κw,j=κr,jfi+κnr,j(1−fi) (2.22g)  vj−1=(1+cτ)BVclqj+μd+θj (2.22h) We consider boundary conditions analogous to those for the monodisperse case:   S^r(0)=S^r0, (2.22i)  S^nr(0)=S^nr,0, (2.22j)  I^(0)=I^0 (2.22k)  D^j(0)=0,j=1,2,⋯,M. (2.22l) 2.4 Approximate solution for a polydisperse aerosol We next determine whether (2.22a), like (2.7) for monodisperse droplet distribution, admit solutions where the number of infected individuals is exponentially decaying due to the presence of RPDs. To that end, we multiply the governing equation for D^j (2.22d) by β˜r,j, and add all the equations:   ddt(β˜r,1D^1+ β˜r,2D^2+⋯)=(κw,1β˜r,1+κw,2β˜r,2+⋯)I^−(β˜r,1ν1D^1+β˜r,2ν2D^2+⋯) (2.23a) An identical operation using the parameter β˜nr,j yields   ddt( β˜nr,1D^1+ β˜nr,2D^2+⋯)=(κw,1β˜nr,1+κw,2β˜nr,2+⋯)I^−(β˜nr,1ν1D^1+β˜nr,2ν2D^2+⋯) (2.23b) We now assume that the droplet removal rate 1/νj is the same for all size bins. This can be an unrealistic assumption, particularly because the gravitational settling rate θis different for different droplet diameters. Several strategies can be employed to overcome this shortcoming. The first is to use in the computations the value of 1/νj yielding the most stringent requirements on the RPDs necessary to avoid an epidemic. This conservative value would be the smallest of the 1/νj values, since 1/νj characterizes droplet removal. Another strategy for making use of a model with identical 1/νj values in all size bins is to use the smallest and largest values of 1/νj in the calculations, in order to bracket the actual infected rate. Assuming 1/νj=1/ν for all j, and setting   D˙r=β˜r,1D^1+β˜r,2D^2+⋯ (2.24a)  D˙nr=β˜nr,1D^1+β˜nr,2D^2+⋯ (2.24b) we obtain the following modification of (2.22a)–(2.22d) and (2.23a) and (2.23b):   ddtS^r=−D˙rS^r (2.25a)  ddtS^nr=−D˙nrS^nr (2.25b)  ddtI^=−ddt(S^r+S^nr)−μII^ (2.25c)  ddtD˙r=δrI^−1νD˙r (2.25d)  ddtD˙nr=δnrI^−1νD˙nr (2.25e) with   δr=κw,1β˜r,1+κw,2β˜r,2+⋯ (2.25f)  δnr=κw,1β˜nr,1+κw,2β˜nr,2+⋯ (2.25g) The relevant boundary conditions are:   S^r(0)=S^r0, (2.25h)  S^nr(0)=S^nr,0 (2.25i)  I^(0)=I^0 (2.25j)  D˙r(0)=0, (2.25k)  D˙nr(0)=0. (2.25l) Equations (2.25a) and (2.25b) can be integrated into forms similar to (2.9) and (2.10). Using those results in (2.25d) and incorporating this modified form of (2.25d) into the time derivative of (2.25e) yields (see steps leading to (2.12)):   d2D˙rdt2+(μI+1ν)dD˙rdt+μIνD˙r=δr[D˙rS^r,0exp(−∫0tD˙r(t′)dt′)+D˙nrS^nr,0exp(−∫0tD˙nr(t′)dt′)] (2.26a)  d2D˙nrdt2+(μI+1ν)dD˙nrdt+μIνD˙nr=δnr[D˙rS^r,0exp(−∫0tD˙r(t′)dt′)+D˙nrS^nr,0exp(−∫0tD˙nr(t′)dt′)]. (2.26b) As with the monodisperse case, we assume that the exponents in (2.26a) and (2.26b) are small and the exponential can be approximated by 1.0. The condition for the validity of this assumption is discussed below. Here we set the exponentials to 1.0 and obtain   d2D˙rdt2+(μI+1ν)dD˙rdt+μIνD˙r−δr[D˙rS^r,0+D˙nrS^nr,0]=0, (2.27a)  d2D˙nrdt2+(μI+1ν)dD˙nrdt+μIνD˙nr−δnr[D˙rS^r,0+D˙nrS^nr,0]=0. (2.27b) Equations (2.27a) and (2.27b) are coupled, precluding a simple solution similar to (2.15). However, a relatively simple solution can be obtained for three important special cases: (i) All susceptible persons initially deploy RPDs, i.e. fs=1. In this case, (2.27a) reduces to   d2D˙rdt2+(μI+1ν)dD˙rdt+μIνD˙r−δrD˙r=0. (2.28a) This equation is similar to (2.13a) and will have exponentially decaying solutions ( D˙r~exp(Mt) ,M<0) when   μIν>δr, (2.28b) i.e.   κw,1β˜r,1+κw,2β˜r,2+⋯<μI⋅1ν (2.28c) Since the simplification was made prior to (2.24a) that the 1/νj values were the same for all bins, in order for (2.28c) to be accurate the smallest (over all bins) 1/νj value, (1/ν)min, should be used. By making this change, and explicitly entering the features of the RPD (analogous to (2.18) in the monodisperse case), we can rewrite (2.28c) as   ∑jκnr,jβ˜nr,jTin,j[Tout,jfi+(1−fi)]<μI(1/v)min (2.28d) When this criterion for exponentially decaying solutions holds, (2.27b) will have solutions of the form D˙nr~texp(Mt) ,M<0, i.e. solutions manifesting some algebraic growth but also exponentially decaying. The number of infected people is (from (2.25d))   I^(t)=I^0(μI−1/ν)2+4δr{M1exp[M1t]−M2exp[M2t]+1νexp[M1t]−1νexp[M2t]} , (2.28e) where   M1,2=−(μI+1/ν)±(μI−1/ν)2+4δr2. (2.28f) (ii) No susceptible individuals initially deploy RPDs, i.e. S^r,0=0,S^nr,0=1 .In this case, (2.27b) reduces to   d2D˙nrdt2+(μI+1ν)dD˙nrdt+μIνD˙nr−δnrD˙nr=0. (2.29a) The criterion for exponentially decaying solutions, analogous to (2.28b), is   μIν>δnr (2.29b) i.e.   κw,1β˜nr,1+κw,2β˜nr,2+⋯<μI⋅1ν (2.29c) Again setting 1/ν equal to the minimum value of 1/νj from all of the bins, and entering the RPD characteristics explicitly, we obtain   ∑jκnr,jβ˜nr,j[Tout,jfi+(1−fi)]<μI(1/v)min. (2.29d) The solution for the number of infected persons is:   I^(t)=I^0(μI−1/ν)2+4δnr{M1exp[M1t]−M2exp[M2t]+1νexp[M1t]−1νexp[M2t]} , (2.29e) where   M1,2=−(μI+1/ν)±}(μI−1/ν)2+4δnr2. (2.29f) (iii) Half of the susceptible population initially deploys RPDs, i.e. S^r,0=1/2,S^nr,0=1/2. Adding (2.27a) and (2.27b) gives,   d2dt2(D˙r+D˙nr)+(μI+1ν)ddt(D˙r+D˙nr)+μIν(D˙r+D˙nr)−12(δr+δnr) (D˙r+D˙nr)=0. (2.30a) The combination field D˙r+D˙nr will have exponentially decaying solutions if   μIν>12(δr+δnr), (2.30b) i.e.   {κw,1(β˜r,1+β˜nr,12)+κw,2(β˜r,2+β˜nr,22)+⋯}<μI⋅1ν (2.30c) or   ∑jκnr,jβ˜nr,j[12Tin,j+12][Tout,jfi+(1−fi)]<μI(1/v)min. (2.30d) The evolution of the infected population is described by:   I^(t)=I^0(μI−1/ν)2+2(δr+δnr){M1exp[M1t]−M2exp[M2t]+1νexp[M1t]−1νexp[M2t]} , (2.30e) where   M1,2=−(μI+1/ν)±}(μI−1/ν)2+2(δr+δnr)2 (2.30f) The criterion for validity of the small-exponent approximation leading to (2.28), derived in the same manner as (2.21c), is   μI⋅1ν−δrδr≫I^0/4 (2.31a) or, writing δr explicitly,   μI•1ν−(κw,1β˜r,1+κw,2β˜r,2+⋯)κw,1β˜r,1+κw,2β˜r,2+⋯≫I^0/4. (2.31b) Similar criteria involving δnr and (δr+δnr)/2 apply for (2.29) and (2.30), respectively. 3. Results The threshold inequality, (2.18), is the essential result for the case of a monodisperse aerosol. For convenience in describing the results we recast it as:   Ωm=[Toutfi+(1−fi)][Tinfs+(1−fs)]{S^0β˜nrκnrμI(1/v)}<1, (3.1) where Ω is the ‘reproduction number’ (Stilianakis & Drossinos, 2010; Diekmann et al., 1990) and the subscript ‘m’ denotes monodisperse. The term is curly brackets is the reproduction number in the absence of protective equipment. Roughly speaking, it is the ratio of the mechanisms promoting infection to those preventing it. For high values of this number, conditions are favourable for the growth of infection. To prevent the growth of infection, the characteristics of the protective strategy, given by the two factors in square brackets, must be such that the overall reproduction number is less than 1. Examining the factors in square brackets, we see that the prevention strategy involves a high degree of compliance in deploying RPDs by the susceptible population ( fs≈1) and a low-inward-transmission-rate barrier for the susceptible population ( Tin≈0), or a high degree of compliance in deploying RPDs by the infected population ( fi≈1) combined with a low-outward-transmission-rate barrier for the infected population ( Tout≈0), or both. The threshold inequality (3.1) was analysed numerically in the computations described below. For the polydisperse case, generalized reproduction numbers and threshold inequalities, corresponding to (2.28d)–(2.30d), can be defined by   Ωp=∑jκnr,jβ˜nr,jTin,j[Tout,jfi+(1−fi)]μI(1/v)min<1(fs=1) , (3.2a)  Ωp=∑jκnr,jβ˜nr,j[Tout,jfi+(1−fi)]μI(1/v)min<1(fs=0) , (3.2b)  Ωp=∑jκnr,jβ˜nr,j[12Tin,j+12][Tout,jfi+(1−fi)]μI(1/v)min<1(fs=12) . (3.2c) These threshold inequalities were also analysed in the computations. Numerical solutions to the governing equations (2.7a)–(2.7h) for the monodisperse case, and (2.25a)–(2.25l) for the polydisperse case, were obtained to validate the approximate analytic solution. The differential equations were solved using a variable time-step Runge–Kutta method, as implemented in the Matlab (MathWorks, Inc., Natick, MA) mathematical software. The analytical approximations for the infected population are given by (2.16b) in the monodisperse case and (2.28d)–(2.30d) for a polydisperse aerosol. The threshold criteria involving the reproduction number (3.1) for the monodisperse case and (3.2a)–(3.2c) for polydisperse were used to prescribe the quality of protection required to avoid a growing infected population, under various conditions. Unless otherwise stated, the parameters used in the calculations are those provided by Stilianakis & Drossinos (2010) for an influenza epidemic. The values are provided in Table 1. For the polydisperse case, we considered two bins, one having a characteristic droplet diameter of 4 microns, and the other a diameter of 8 microns. In all computations, the initial infected population was 0.1% of the total population. Parameters characterizing the protective equipment were varied in the calculations. Table 1. Parameter values used in calculations Parameters  Values  Breathing rate, B ( m3/day)  24  Volume of the personal cloud of an infected person, Vcl ( m3)  8  Characteristic breathing time, τ(days)  0.0139  Probability of infection by an inhaled pathogen, p  0.052  Droplet diameters, d1 and d2(μm)  4 and 8  Inhaled-droplet deposition probability, q1 and q2  0.88 and 0.99  Number of pathogens per droplet, Np1 and Np2  9.15×10−4 and 8.2×10−3  c (per day)  13  Inactivation rate of airborne pathogens, μd (per day)  8.64  Gravitational settling rate, θ1 and θ2 (per day)  28.8 and 113.2  Infection recovery rate, μI (per day)  0.2  Droplet production rate for two particle sizes without respirator, κnr,1κnr,2  4.1×105, 1.92×104  Fraction of initial number of infecteds, I0  0.001  Fraction of incoming pathogens transmitted by the RPD, Tin  0.25  Outward transmission rate, Tout  0.1  Parameters  Values  Breathing rate, B ( m3/day)  24  Volume of the personal cloud of an infected person, Vcl ( m3)  8  Characteristic breathing time, τ(days)  0.0139  Probability of infection by an inhaled pathogen, p  0.052  Droplet diameters, d1 and d2(μm)  4 and 8  Inhaled-droplet deposition probability, q1 and q2  0.88 and 0.99  Number of pathogens per droplet, Np1 and Np2  9.15×10−4 and 8.2×10−3  c (per day)  13  Inactivation rate of airborne pathogens, μd (per day)  8.64  Gravitational settling rate, θ1 and θ2 (per day)  28.8 and 113.2  Infection recovery rate, μI (per day)  0.2  Droplet production rate for two particle sizes without respirator, κnr,1κnr,2  4.1×105, 1.92×104  Fraction of initial number of infecteds, I0  0.001  Fraction of incoming pathogens transmitted by the RPD, Tin  0.25  Outward transmission rate, Tout  0.1  View Large 3.1 Numerical results for monodisperse case In Fig. 2, the size infected population is plotted as a function of time, for a monodisperse aerosol. Different values of the reproduction number (3.1), generated by considering different values of the droplet production rate κnr, are considered. In the calculations it was assumed that 50% of the susceptible population deploys RPDs ( fs=50%), and that the RPDs allow 50% of the aerosols incoming into the RPD during breathing to be transmitted ( Tin=50%). A Tin value of 50% is representative of a respirator that has not been fit-tested (Coffey et al. 2004). The fraction of infected persons deploying RPDs was 10% ( fi=10%) and the outward transmission rate was 10% ( Tout=10%, Patel et al., 2015). Decaying infected populations can be seen for the reproduction numbers less than 1, while growing infected populations (at least initially) can be seen for reproduction numbers greater than 1. For reproduction numbers less than 1, the analytical approximation closely matches the numerical one for the entire 10-day time period. For Ωm=3, the analytical approximation also closely approximates the numerical one for the 10-day duration, even though the infected field is growing. For the higher reproduction number ( Ωm=6), the analytical approximation tracks the numerical one closely for about half the 10-day period, but then fails to predict the turnover in the infected-population count. Fig. 2. View largeDownload slide Comparison of the analytical and numerical solutions for the infected population, in situations where the monodisperse reproduction number is less than 1 (bottom curves) and greater than 1 (upper curves). Fig. 2. View largeDownload slide Comparison of the analytical and numerical solutions for the infected population, in situations where the monodisperse reproduction number is less than 1 (bottom curves) and greater than 1 (upper curves). A subsequent set of calculations was performed in which it was assumed that 50% of the susceptible population deploys RPDs ( fs=50%), and that the fraction of infected individuals deploying RPDs was zero ( fi=0). It was determined that for an inward transmission rate of slightly greater than 16%, the threshold inequality was just satisfied ( Ωm=1) and the infected population should not grow. In other words, as long as the RPDs deployed by the susceptible population were able to contain 84% of the pathogens, even if the infected individuals were not protected, a pandemic would not spread. This prediction was checked by full numerical solution, using Tin=16% as well as values slightly lower and higher, namely 15, 17 and 18%. The results are shown in Fig. 3. The dynamics of the infected population follow the prediction of the threshold inequality, with Tin values of 15 and 16% leading to decay, and 17 and 18% leading to initial growth. The initial growth rate, as well as the turn-over time, increased significantly with increasing barrier transmission rate. Fig. 3. View largeDownload slide Infected population as a function of time, when the susceptible population deploys RPDs having different inward transmission rates. Curves were computed from full numerical model for a monodisperse droplet. Analytical threshold inequality predicts that the transition between a growing infected population and a decaying one occurs for Tin slightly above 16. Fig. 3. View largeDownload slide Infected population as a function of time, when the susceptible population deploys RPDs having different inward transmission rates. Curves were computed from full numerical model for a monodisperse droplet. Analytical threshold inequality predicts that the transition between a growing infected population and a decaying one occurs for Tin slightly above 16. The threshold inequality (3.1) was used to identify the level of integrity required of the barriers worn by the susceptible population, in order to avoid a growing infection rate, as a function of the level of compliance in wearing RPDs by the infected population. The outward transmission rate was taken to be Tout=10%. Different levels of compliance by the susceptible population ( fs=0.1, 0.3, 0.5, 0.7 and 0.9) were considered. The results are shown in Fig. 4. As expected, when the compliance rate of the infected population increases (higher fi), lower-integrity barriers (higher Tinvalues) worn by the susceptible population can still prevent growth in the infection population. When the compliance of the infected population exceeds approximately 0.14 (for the parameter values considered), the number of pathogens produced is low enough that infection growth does not occur regardless of the behaviour of the susceptible population. That is, for any level of compliance by the susceptible population (all fs values), and the use of barriers that transmit all pathogens ( Tin=100%), infection growth does not occur. This is manifested in Fig. 4 in the convergence of all the curves to the point ( fi=0.14, Tin=100%). For a compliance rate by the susceptible group of 10% ( fs=0.1), the curve originating at ( fi=0.03, Tin=0) implies that for compliance rates by the infected group of less than 3%, a growth in the infected population cannot be prevented except by perfect protection ( Tin=0) of the susceptible population. Fig. 4. View largeDownload slide Level of barrier integrity (designated by the transmission rate Tin) required in order to avoid a growing infection rate, as a function of the fraction of infected persons deploying RPDs. A monodisperse droplet distribution was assumed. Curves are plotted for different fractions of the susceptible population initially wearing RPDs. Outward transmission rate for the infected population is 10% ( Tout=0.1). Fig. 4. View largeDownload slide Level of barrier integrity (designated by the transmission rate Tin) required in order to avoid a growing infection rate, as a function of the fraction of infected persons deploying RPDs. A monodisperse droplet distribution was assumed. Curves are plotted for different fractions of the susceptible population initially wearing RPDs. Outward transmission rate for the infected population is 10% ( Tout=0.1). 3.2 Numerical results for polydisperse case Figure 5 for the polydisperse case is analogous to Fig. 2 for the monodisperse. The barrier parameters were the same ( fi=10%, Tout=10%, fs=50%, Tin=50%) as in Fig. 2, and again the reproduction number Ωp was varied by adjusting the droplet production rate κnr The infected population decreases as a function of time for reproduction numbers less than 1 (bottom two sets of curves in Fig. 5), and increases for reproduction numbers greater than 1 (top two sets of curves). As with the monodisperse droplets, the analytical approximation matches the numerical one closely for the decaying infected populations. For the increasing populations, the analytical approximation matches the numerical one well for about half of the 10-day duration, but does not predict the maximum and eventual decrease in the population. Fig. 5. View largeDownload slide Comparison of the analytical and numerical solutions for the infected population, in situations where the polydisperse reproduction number is less than 1 (bottom curves) and greater than 1 (upper curves). Fig. 5. View largeDownload slide Comparison of the analytical and numerical solutions for the infected population, in situations where the polydisperse reproduction number is less than 1 (bottom curves) and greater than 1 (upper curves). Figure 6 is the polydisperse equivalent to Fig. 4 for the monodisperse case. As indicated above, for the polydisperse distribution, a second size bin (8-micron) was added to the first bin (4-micron). There are thus additional inhalable droplets for the polydisperse case, making the challenge against RPDs deployed by the susceptible population greater. The allowable transmission rate ( Tin) to prevent the growth of infection is therefore lower for the polydisperse case. When no susceptibles deploy RPDs ( fs=0), the transmission rate is immaterial, and hence the result is a vertical line. The fi value (roughly 0.45) at which the vertical line is plotted is the threshold value for compliance by the infected population (assuming a Tout of 0.1), above which there will be decaying infected populations. Fig. 6. View largeDownload slide Level of barrier integrity (designated by the transmission rate Tin) that is required in order to avoid a growing infection rate, as a function of the fraction of infected people deploying RPDs, for a polydisperse distribution. Curves are plotted for different fractions of the susceptible population initially wearing RPDs. Outward transmission rate for the infected population is 10% ( Tout=0.1). Fig. 6. View largeDownload slide Level of barrier integrity (designated by the transmission rate Tin) that is required in order to avoid a growing infection rate, as a function of the fraction of infected people deploying RPDs, for a polydisperse distribution. Curves are plotted for different fractions of the susceptible population initially wearing RPDs. Outward transmission rate for the infected population is 10% ( Tout=0.1). 4. Discussion This work extends the SIR model of Stilianakis & Drossinos (2010) for predicting the spread of infectious disease by inhalable droplets to include effect of protective measures. The most important use made of the extended equations was to derive an analytical threshold inequality for predicting the level of protection required from RPDs in order to avoid a growing infection rate. This inequality (3.1) for the monodisperse case and (3.2a,b,c) for polydisperse aerosols, is written in terms of a generalized reproduction number containing the properties of the RPDs and the level of compliance of the susceptible and infected populations. The predictions of threshold for growth from the analytical model agreed well with full numerical solutions to the SIR equations. Another useful result emerging from the analytical approximation of the SIR equations is the ability of the approximate model to track the initial dynamics of the infected population (Figs 2 and 5), even in the case of reproduction numbers substantially bigger than 1. This ability could prove useful, e.g. in identifying the time interval available for medical assistance to arrive before the infected population reaches a certain size, as a function of the properties of the pathogen, the level of protection, and the characteristics of the susceptible population. An initial normalized infected population of I^0=0.001 was used in our calculations. To the extent that the linearized equation (2.13a) is valid, the results (e.g. threshold equation) are independent of the initial infected population (solutions for different initial conditions are the same to within a scaling factor). For decaying solutions, the criterion for validity of the small-exponent approximation, and hence of the linear equation (2.13a), was shown to be (2.21c). For the parameters considered in the calculations (including Tin=Tout=50%, fs=fi=50%), this condition is roughly I^0≪1. For growing solutions, the analysis leading to (2.21c) is not applicable, as the exponential terms grow large with time. However, for sufficiently small times the small-exponent approximation is still valid with growing solutions. If we consider the exponent, provided in (2.19c), and take the largest of the terms exp[ M1] and exp[ M2], we can say that the exponent will be small if   β˜rκwI^0M2−M1exp[M1t]M1≪1 (4.1a) or   t≪1M1ln[M1(M2−M1)β˜rκwI^0] . (4.1b) For the parameters used in the calculations, this inequality becomes, roughly, t≪30 days. This criterion is consistent with the results of Figs 2 and 5, where the analytical approximation correctly predicted the infection dynamics for the growing population over a period of at least a few days. As noted above, this period is useful for identifying the initial growth rate of an epidemic, and devising an initial intervention strategy. The symmetry in (3.1) between the parameters ( Tin,fs) characterizing the level of protection for the susceptible population and the parameters ( Tout,fi) characterizing the level of protection for the infected population indicate that, formally, deployment of protective equipment by one population is as effective in reducing the likelihood of an epidemic as deployment by the other. However, as reflected in the parameter selection in the computations ( Tout=10%, Tin=50%) protective equipment may be more effective in reducing the outward flux of pathogens than the inward flux (Coffey et al., 2004; Patel et al., 2015). In such a case, preferential attention to the infected population (e.g. to bring about a higher level of compliance) may be a worthwhile strategy. In the derivation of the analytical approximation for the polydisperse aerosol, including the threshold equations (3.2a)–(3.2c), it was assumed that the droplet removal rate 1/ν was the same for all droplet sizes. This assumption led to a set of equations, (2.25), that was reduced in size and more mathematically tractable than the full set of polydisperse equations, (2.22). Choosing the smallest removal rate over all size bins as the one removal rate, as we did in our calculations, results in a conservative estimate of the infection growth rate. That is, the number of infected persons determined through a numerical or analytical approximation to (2.25) will be greater than or equal to the number of infected persons determined by solution to (2.22). For the bimodal set of data used in our polydisperse calculations, the difference between assuming a uniform removal rate and using different removal rate for the two size bins can be gleaned from Fig. 7. The conditions match those of Fig. 5, for a reproduction number of 4. The predicted number of infected individuals is substantially lower (about a factor of 4) for the variable removal rate, owing to the very large settling rate for the 8-micron size. In practice, if highly conservative values are to be avoided in favour of more accurate predictions, average droplet removal rate can be used instead of the smallest, or calculations using the lowest and highest removal rates can be performed to establish a range of infected-population values. In any event, we highlight the desirable property that the analytical approximation leads to a conservative estimate of the infected population. Fig. 7. View largeDownload slide Infected population as a function of time determined three ways: (1) numerical solutions the full polydisperse equations (25); (2) numerical solution of the polydisperse equations assuming uniform droplet removal rate (22) and (3) approximate analytical solution to the polydisperse equations for uniform droplet removal rate (22). Reproduction number is equal to 4. Fig. 7. View largeDownload slide Infected population as a function of time determined three ways: (1) numerical solutions the full polydisperse equations (25); (2) numerical solution of the polydisperse equations assuming uniform droplet removal rate (22) and (3) approximate analytical solution to the polydisperse equations for uniform droplet removal rate (22). Reproduction number is equal to 4. The majority of the effort in this study was devoted to the development of the equations and approximations for the polydisperse case. The value of the polydisperse formulation appeared in the example just given where a significant (factor of 4) difference in the infected population relative to the monodisperse model was found. Additionally, the polydisperse formulation was developed in order to handle more complicated distributions than the present bimodal one. Nicas et al. (2005) list more than 10 size bins in the 1−200 μm range where the number of cough droplets produced in each bin varies by more than 10% from the previous bin. There are 4 bins containing more than 50% of the maximum number of droplets. Even if the size dependence of some parameters is unknown and an average (effectively monodisperse) value must be used, the ability to leverage the known size dependence of important quantities such as the droplet production rate in a model provides more semblance of reality to the calculations. Many of the calculations performed were limited to a short time (e.g. 10 days in Figs 2, 5, 7) relative to the actual duration of a pandemic in which the protective equipment might be deployed. The duration was kept short in order to focus on the applicability of the analytical model in the challenging case of a reproduction number significantly larger than 1. There is no time limit at which the accuracy of the full model ((2.7) for the monodisperse case and (2.22) for polydisperse), or the approximate model for reproduction number less than 1, degrades. The calculations featured in Fig. 3, e.g. cover the full range of infection dynamics: growth, peak and decay. These aspects of disease transmission occur over a period of 100 days. Ideally, the results generated by the present model should be validated in experiments involving protective devices. It should be recognized that validating even standard SIR models, which involves deliberate infection of a susceptible population, is extremely difficult. Experiments involving controlled use of protective equipment are even more complicated. The feasibility of attaching protective devices to animals is likewise low. We do note that the SIR equations were modified in a very systematic manner to account for the effect of RPDs, by adjusting the flux of droplets input to and output from of the respiratory system. Thus, other methods devised to test SIR models that can also modify the inflow or outflow of pathogens, e.g. with pharmacological interventions or by controlling the air quality in an environment, can be used to indirectly validate the present model. In the absence of rigorous validation, the accuracy of the model's absolute predictions for the number of persons infected is not known. However, the accuracy of relative predictions is likely to be more reliable, and such relative predictions represent an important use for the model. In the purchasing and stockpiling of protective equipment for a possible spread of infection, example, a valuable piece of information is an estimate of the relative risk reduction achieved when one possible barrier is replaced by another, more expensive one with a lower transmission rate. Another useful estimate that the model can provide is the degree of compliance by the target population that would be required for an RPD to produce a decay (rather than growth) in the infected population. In such cases, sociological data may suggest that the level of compliance is not achievable, and resources would not be well spent pursuing a strategy involving the given RPD for the given population. The computations in this paper were performed primarily to validate the analytical approximation and explore the dynamics of the infected population near the threshold for infection growth. Much remains to be done using the generalized equations derived in this paper in order to assess the impact of protective measures during an epidemic (i.e. well above threshold). The full differential equations should be numerically solved with the relevant parameters informed by data from actual barrier tests (e.g. Guha et al., 2016). The uncertainty in the predictions must also be quantified. In this regard, the analytical approximation (e.g. (2.16)) is useful. The variability in a prediction, say for the infected field I, with respect to a given parameter p may be written   δI=∂I∂pδp (4.2) where δp is the uncertainty in the parameter p. The sensitivity ∂I∂p may be obtained through direct differentiation of the analytical approximation. Calculations of the infected and susceptible populations under epidemic conditions, as well as the uncertainties in the predictions, are currently underway. Funding This work was supported by the FDA Medical Countermeasures Initiative, Grant Numbers MCM2D20511 and MCM2J270HT. References Chao C. Y. H. Wan M. P. Morawska L. Johnson G. R. Ristovski Z. D. Hargreaves M. Mengersen K. Corbett S. Li Y. Xie X. Katoshevski D. ( 2009) Characterization of expiration air jets and droplet size distributions immediately at the mouth opening. J. Aerosol. Sci. , 40, 122– 133. Google Scholar CrossRef Search ADS   Chen S.-C. Liao C.-M. ( 2008) Modelling control measures to reduce the impact of pandemic influenza among schoolchildren. Epidemiol. Infect. , 136, 1035– 1045. Google Scholar CrossRef Search ADS PubMed  Coburn B. J. Wagner B. G. Blower S. ( 2009) Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1). BMC Med.,  7, 30. Google Scholar CrossRef Search ADS   Coffey C. C. Lawrence R. B. Campbell D. L. Zhuang Z. Calvert C. A. Jensen P. A. ( 2004) Fitting characteristics of eighteen N95 filtering-facepiece respirators. J. Occup. Environ. Hyg. , 1, 262– 271. Google Scholar CrossRef Search ADS PubMed  Committee on the Certification of Personal Protective Technologies. ( 2011) Certifying Personal Protective Technologies: Improving Worker Safety  Washington DC: Institute of Medicine, National Academy of Sciences. Diekmann O. Heesterbeek J. A. P. Metz J. A. J. ( 1990) On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J. Math. Biol. , 28, 365– 382. Google Scholar CrossRef Search ADS PubMed  Furuya H. ( 2007) Risk of transmission of airborne infection during train commute based on mathematical model. Environ. Health. Prev. Med. , 12, 78– 83. Google Scholar CrossRef Search ADS PubMed  Guha S. McCaffrey B. Hariharan P. Myers M. ( 2016) Quantifying the leakage through gaps in surgical masks and pediatric face masks. J. Occup. Hyg. (to appear). Huang S.-H. Chen C.-W. Chang C.-P. Lai C.-Y. Chen C.-C. ( 2007) Penetration of 4.5 nm to 10 μm aerosol particles through fibrous filters J. Aerosol Sci. , 38 719– 727. Google Scholar CrossRef Search ADS   Lee S. A. Grinshpun S. A. Adhikari A. Li W. McKay R. Maynard A. Reponen T. ( 2005) Laboratory and field evaluation of a new personal sampling system for assessing the protection provided by the N95 filtering facepiece respirators against particles. Ann. Occup. Hyg. , 49 245– 257. Google Scholar PubMed  Mniszewski S. M. Del Valle S. Y. Priedhorsky R. Hyman J. M. Hickman K. S. ( 2014) Understanding the impact of face mask usage through epidemic similation of large social networks. Theories and Simulations of Complex Social System (D. Vahid and M. Vijay Kumar eds), 97–115. Nicas M. Nazaroff W. W. Hubbard A. ( 2005) Towards understanding the risk of secondary airborne infection: emission of respirable pathogens. J. Occup. Environ. Hyg. , 2, 143– 154. Google Scholar CrossRef Search ADS PubMed  Patel R. B. Skaria S. D. Mansour M. M. Smaldone G. C. ( 2015) Respiratory source control using a surgical mask: an in vitro study. J Occup. Environ. Hyg. , 15, 1– 29. Stilianakis N. I. Drossinos Y. ( 2010) Dynamics of infectious disease transmission by inhalable respiratory droplets. J. R. Soc. Interface , 7, 1355– 1366. Google Scholar CrossRef Search ADS PubMed  Sze To G. N. Chao C. Y. H. ( 2010) Review and comparison between the Wells–Riley and dose-response approaches to risk assessment of infectious respiratory diseases. Indoor Air , 20, 2– 16. Google Scholar CrossRef Search ADS PubMed  Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications 2016. This work is written by (a) US Government employee(s) and is in the public domain in the US.

Journal

Mathematical Medicine And Biology: A Journal Of The ImaOxford University Press

Published: Mar 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off