General

If the pathogen concentration in an infectious individual’s respiratory tract fluid ρ_p is low enough, almost all exhaled pathogen copies will be the only pathogen in their aerosols, i.e. mono-multiplicity aerosols, and poly-multiplicity aerosols can reliably be ignored. We will use the tailing subscript k to denote aerosols with k pathogen copies inside them. An aerosol cannot contain more pathogen copies than will fit in its volume, and there is a limit to how large an aerosol/droplet a person can exhale. Let M be the maximum number of pathogen copies that can fit in the largest aerosol/droplet that can possibly be exhaled. This is the hard cutoff/limit on k. There also exists a soft cutoff/limit M_c ≤ M for which contributions of aerosols with k > M_c is negligible. In a worst case M_c = M, but in practice it is much lower since the pathogen volume fraction of respiratory tract fluid is quite low even at the upper pathogen load for some diseases and the largest droplets don’t stay airborne and ballistically fall to the ground. For example, SARS-CoV-2 at the very upper end of its concentration range at 10¹¹ cm^-3 [35, 36] would give a volume fraction of approximately 5 × 10⁻⁵, if we treat the virus as a 100 nm sphere (approximate size of the SARS-CoV-2 virus [37]). This is important because an aerosol with a diameter of 1 μm could contain up to approximately 740 spherical pathogen copies with diameter 100 nm, if we assume hard-sphere packing (packing fraction of 74%). An aerosol with a diameter of 10 μm could contain up to approximately 7.4 × 10⁵ of the same pathogen copies for the same packing fraction.

To properly account for higher multiplicities, we must consider the separate doses for each multiplicity. Let Δ_k be the number of pathogen copies absorbed from aerosols with multiplicity k, and let m_k be the number of aerosols absorbed with multiplicity k. The aerosol and pathogen doses are related by Δ_k = km_k. The total pathogen dose from all aerosols is just the sum of the doses for each multiplicity, which is $$ \Delta ={\sum }_{k=1}^{\infty }{\Delta }_k$$. Let μ_k = 〈m_k〉 = 〈Δ_k〉/k be the average number of absorbed aerosols with multiplicity k.

As long as the aerosols are randomly distributed in space (well-mixed with no clustering nor avoidance), then the PDF of each m_k follows a Poisson distribution with mean μ_k. Since Δ_k = km_k, the PDF of Δ_k is not a Poisson distribution for k > 1. It is instead a scaled-Poisson distribution of the form

The deviation from the Poisson distribution is most visible in the fact that this distribution has holes. For example with k = 2, P_k = 0 for all odd Δ_k. Since Δ is the sum of a Poisson distribution for k = 1 and some number of possibly non-negligible scaled-Poisson distributions, the PDF of Δ will not be a Poisson distribution unless the contributions from k > 1 are negligible compared to k = 1. So we can’t just naively put the expected average dose into dose-response models expecting a Poisson distribution.

Instead, we must change the summation in Eq (2) to get the infection risk $$ \mathfrak{R}$$. Let us consider the p’th moment, $$ {\mathcal{M}}_p$$, of the infection probabilities as a function of the average aerosol doses μ_k (note, we use p in later sections of this manuscript as a summation index). To determine $$ {\mathcal{M}}_p$$, we must sum over all possible combinations of exact aerosol doses m_k of each multiplicity for k ∈ [1, ∞) of the product of the Poisson probabilities of each m_k and the infection risk for the dose raised to the power of p. This is

Exponential model corrections

Then, putting R_E from Eq (3) into Eq (11), the exponential model mean infection risk is

Beta-poisson model corrections

The integral over r commutes with the sums in Eq (10). So as was with the case when multiplicity is not considered in Eq (7), we can get the moments by taking the result for the exponential model and integrating it times the beta distribution PDF over r. This is

Unfortunately, as is the case for when the dose is Poisson distributed [13], the integral cannot be solved analytically and must be solved numerically or approximated, though now it is harder with the extra terms for M_c > 1.

Looking ahead

Now that we have dose-response models corrected for the multiplicity via Eq (11), we must determine the average aerosol doses μ_k for each multiplicity before the infection risk can be calculated. We now generalize the Wells-Riley formulation for multi-pathogen aerosols to get this. In the following sections, we will describe the environment, people, aerosols, sources, sinks, etc. to get the model equations. Let n_k(d₀, t) be the concentration density of aerosols with original diameter d₀ (diameter at production) and k pathogen copies in them over time, which has units of [L]^-4 where [L] is the unit of length since n_k(d₀, t)dd₀ is the concentration of infectious aerosols with original diameters between d₀ and d₀+ dd₀. To get a concentration, n_k(d₀, t) must be integrated with respect to d₀.

In the end, we will get the following system of ODEs (Ordinary Differential Equations) in time t and the original diameter at production d₀ for the n_k, which is

Table 1

Summary of all the source (the β) and sink (the α) terms considered in this manuscript. “Individuals” is abbreviated as “ind.” See their relevant sections for details on where they come from and the meanings of their terms.

Source and sink term summary.

Term	Meaning	Form
βr,k(d0, t)	transport from other rooms	qr(t)nr,k(d0, t)
βI,k(d0, t)	production by infectious individuals	$$
αo(t)	air exchange with outside	qo(t)
αr(t)	air exchange with other rooms	qr(t)
αv(d0, t)	filtering by ventilation	qv(t)Ev(w(d0, t)d0)
αg(d0, t)	gravitational settling	$$
αd(d0, t)	deposition on surfaces	found elsewhere
αI,f(d0, t)	filtering by infectious ind. inhaling	$$
αS,f(d0, t)	filtering by susceptible ind. inhaling	$$
αO,f(d0, t)	filtering by other ind. inhaling	$$

Environment

Like most Wells-Riley formulations, we consider the infection risk in one sufficiently well-mixed indoor environment such as a room or set of rooms sufficiently coupled together with respect to their air that they have the same infectious aerosol concentration densities. And we assume that sources, sinks, and individuals are far enough apart from each other that the local concentration densities at their locations are approximately equal to the average concentration density in the whole environment. Note that the particular kind of ventilation has an impact on the validity of this assumption [30]. See the Discussion for when this assumption is not valid. The environment could also be split into coupled well-mixed zones with weaker mixing between them [7, 33], but that shall not be considered here.

Let the volume of the environment be V. Air is exchanged with outside, with other rooms, and circulated internally through the ventilation system. Let Q_o, Q_r, and Q_v be the volumetric rate of air exchange with outdoors, other rooms, and the circulating ventilation of the environment (ventilation system that pulls air out of the environment and puts it back in). These will be normalized by the environment volume; yielding q_o ≡ Q_o/V, q_r ≡ Q_r/V, and q_v ≡ Q_v/V since target values of these parameters are often the design goals for HVAC systems.

Aerosols

Consider the concentration of infectious aerosols over time. To be completely accurate, we need to consider the concentration density for each multiplicity k as a function of time, current diameter d while in the environment, and the solute content (including inactivated pathogen copies). We have to consider both d and the solute content because an exhaled aerosol’s equilibrium diameter is a function of its solute content, the humidity, and the temperature [27]. Higher solute concentrations decrease the vapor pressure of the aerosol, which allows equilibrium to be reached as long as the environment isn’t super-saturated or too close to saturated [26, 27]. For higher humidities, an aerosol will continue to grow by condensation indefinitely, though the growth rate slows towards a crawl for d > 20 μm [26, 38]. But such super-saturated conditions can cause clouds/fog, which rarely occur in indoor environments. So we will assume the environment is sub-saturated. If the environment is dry, the aerosols can evaporate at most to the point where they are purely precipitated solid with no water left. Note that as a drop (whether large or a small aerosol) dries, the solute fraction increases, until at some point the solute makes the shape non-spherical (not enough water to spherically encapsulate the insoluable components, solute causing anisotropy and/or inhomogeneity in the surface tension, etc.). This will occur at a humidity no lower than the efflorescence relative humidity of the solute mix, where the soluble solutes will homogeneously nucleate and the water completely evaporates away. Infectious aerosols always have at least two components of the solute (whatever is in the respiratory tract fluid plus the pathogen/s), so there is the possibility of heterogeneous nucleation causing the water to completely evaporate away at a higher humidity.

This means that we have four different diameters to consider, which are

dcurrent diameter in the environment (spherical equivalent diameter if it is completely dry or almost dry and the solute causes a non-spherical shape)

d_eequilibrium diameter in the environment

d₀wet diameter at production (original diameter), which determines the distribution of initial multiplicities

d_Dspherical equivalent dry diameter when all water is evaporated away and just solute remains (note that the aerosol may no longer be spherical, so the spherical equivalent diameter for the same volume must be used)

For any aerosol; d₀ and d_D are fixed and never change as long as collisional-coalescence and shattering don’t occur (can be treated as fixed if these processes are negligible), d_e is dynamic in time if the environment’s temperature and/or humidity changes, and d is dynamic in time unless the environment’s temperature and humidity exactly match those inside the respiratory tract at the point of production.

Small wet/nucleated aerosols respond very quickly to the humidity and temperature, evaporating/condensing to their equilibrium diameter in a very short period of time due to their high surface area to volume ratio [9, 25, 26]. Assuming the environment is well-mixed enough that the time between exhalation from an infectious individual and inhalation by any person is long compared to the evaporation/condensing time scale, we can make the approximation that all aerosols are at their equilibrium diameter when in the environment (d ≈ d_e). This means that when d_e increases from d_e = d_N (completely evaporated) to d_e > d_N (wet/nucleated), we are assuming that the time the aerosols require to nucleate and grow to d_e is short compared to other time scales in the model and therefore also make the approximation d ≈ d_e even when d_e increases from d_e = d_N to d_e > d_N. This means that we just need to worry about the equilibrium diameter and its changes, and not the non-immediate response to shifting equilibrium diameters. There is one complication, however. Aerosols will initially stay in the exhaled plume where the humidity is higher, so they won’t reach the well-mixed equilibrium diameter till they leave the plume or the plume is diluted and mixed with the environment, which brings us back to the well-mixed environment assumption.

We will also make the assumption that the temperatures and humidities in different individuals’ respiratory tracts (and the volume under their facemasks if they are wearing any) are similar enough and change negligibly enough over time that the equilibrium diameter in people’s respiratory tracts is d₀. If the aerosols have not completely dried out in the environment (d_e > d_D), the aerosols will start to grow inside people’s respiratory tracts back towards d₀ and thus d ≈ d₀ inside the respiratory tract. But the time scale of breathing is short and for completely dried out aerosols it takes time to nucleate and grow back to d₀, which means that some fraction of dry aerosols might not reach d = d₀ while in the respiratory tract. However, we will make the assumption/approximation that dry aerosols have returned to their original diameter by the time they are exhaled back out if they were not absorbed in respiratory tract. This last approximation only affects the sink from individuals inhaling aerosols α_C,f(d₀, t) from Eq (42) if they are wearing masks, which is usually small compared to other sinks. When the individuals in the environment are wearing masks and the α_C,f(d₀, t) sinks dominate, then a better approximation or an explicit treatment of the diameter when exhaled should be used. Combined, our assumption/approximation is

Let us define ratios between the remaining diameters: the evaporation ratio w, the dilution ratio δ, and the initial solute ratio ζ as

Note that w and δ are potentially functions of time, as well as diameter due to the effect of surface curvature (through surface tension) on equilibrium vapor pressure [26, 27]. Also, different solutes have different molar densities, different practical osmotic coefficients, and maximum concentrations before they precipitate; and therefore different functional relationships between the saturation vapor pressure and the concentration [27]. So different solute compositions will cause w and δ to be different even for aerosols with the same ζ.

But, we will make the assumption that the value of ζ and the solute composition (except for the pathogen copies) is approximately constant from each infectious individual to the next and over time with each infectious individual, and we will ignore the contribution of the pathogen copies (both active and inactivated) to the equilibrium vapor pressure and therefore d_e. We will also assume that ζ has no diameter dependence (i.e. attraction and repulsion of solutes from the liquid surface at production has a negligible effect on solute fraction and composition). With these approximations, we have a single constant value of ζ and single functions for w and δ, possibly over time and d₀ (or equivalently d_D), for all infectious aerosols in the environment.

This means we can choose to track one of d_e, d₀, or d_D and always know the other two through the ratios that are the same for all infectious aerosols at the same moment of time with the same value of the chosen diameter parameter. Thus we have two independent variables, t and one diameter parameter.

Processes such as gravitational settling, deposition, filtering or exchange by the ventilation, filtering by facemasks when inhaling are all functions of the current diameter, which is approximately d_e, making d_e convenient. Additionally, any non-drying aerosol instrument can be used in the environment to measure d_e. But, because d_e can change over time for a fixed d_D or d₀, the equations for the aerosol concentration density in terms of t and d_e have a flux term (from evaporation/growth) with a partial derivative with respect to d_e; making the equations PDEs (Partial Differential Equations) which adds complications in the analysis. This can be seen by considering the total time derivative of the aerosol concentration density $$ \tilde{n}$$ expressed in terms of t and d_e, which is

Since d_D and d₀ are fixed for a given aerosol over time regardless of how the temperature or humidity in the environment might be changing, the equivalent flux term is zero and thus the equivalent functions are ODEs, which are much easier to solve. Thus, we eliminate d_e as a choice for the diameter parameter.

The model in this manuscript can be constructed with either choice of d₀ or d_D, with w appearing in places if d₀ is chosen, and both δ and ζ appearing in places if d_D is chosen. We choose d₀ because then we only need one of the ratios (w only), the diameter limits are easier to express in it, and the literature on the diameter distributions of exhaled aerosols generally work hard to convert their measurements (vary between whether they are d_e or d_D) into expressions in terms of d₀ rather than d_D.

Now, n_k(d₀, t) is the concentration density of aerosols in terms of t and the original diameter d₀. Let $$ {\tilde{n}}_k$$ be the concentration density in terms of t and d_e, and $$ {\breve{n}}_k$$ be the concentration density in terms of t and d_D. To make conversions between them; consider the original diameter interval d₀ to d₀+ dd₀, and its corresponding intervals d_e to d_e + dd_e and d_D to d_D + dd_D. The number of aerosols in each interval must all be equal: n_k dd₀, $$ {\tilde{n}}_{k}d{d}_e$$, and $$ {\breve{n}}_{k}d{d}_D$$. Thus, the conversions are

Let n_0,k(d₀) be the initial concentration density in the room for a multiplicity k at the initial time t = t₀ and n_r,k(d₀, t) be the volume averaged concentration density of the air coming in from other rooms. We are assuming that the concentration density outdoors is negligible.

Diameter limits

For the model, we will limit ourselves for each multiplicity to the range d₀ ∈ [d_m,k, d_M] where d_m,k is the minimum aerosol diameter required to hold k pathogen copies, and d_M is a diameter cutoff separating larger aerosols that are more ballistic and gravitationally settle to the ground too quickly to become well mixed and smaller aerosols that more closely follow the flow and mix. Let K_m(d₀) be the largest number of pathogen copies that can fit in an aerosol at production. We will consider

For a spherical pathogen with diameter d_p, we can use the crude approximation of just considering the total pathogen volume and a packing efficiency e = 0.74 (hard pack spheres) with a minimum of 1 and completely neglect the aerosol shape that small number of pathogen copies would force (two pathogen copies, for example, can’t be arranged into a configuration that even vaguely resembles a sphere). We can use the same idea to get K_m(d₀). Both of them are

At the lower limit near d_m,k, the pathogen/s take up a disproportionate amount of the space in the aerosol compared to other solutes and the assumption of approximately equal solute concentrations at production is violated and the evaporation ratio has a strong dependence on d₀ and the initial multiplicity, the latter of which we aren’t tracking at all. However, as long as the total liquid volume of exhaled aerosols with diameters close to d_m,k (say, those whose diameters are small enough that their volume is only a few times larger) is small compared to total liquid volume of the rest of the range in d₀, this problem will have a negligible effect. Additionally, the diameter dependence of many of the sink terms may be much smaller close to d_m,k for submicron pathogen copies which means that the effect of assuming the wrong evaporation ratio may be small. The smaller the pathogen, the less issues this will pose. It will be least important for small viruses, and possibly quite important for large bacteria and eukaryotic pathogens.

The upper limit is rather imprecise since there is no single hard separation scale that could be chosen unless the air is completely still in which case one can use a so called “Wells curve” (same Wells as of the Wells-Riley model) for the environment’s humidity to determine the largest size that won’t settle to the ground before evaporating to their equilibrium diameter, such as the original one [24] or newer ones [25]. But mixing of any sort complicates this. One might think that one could just rely on the fact that the gravitational settling sink term keeps growing with diameter and not bother with the problem. But, the well-mixed assumption breaks down and the lifetime of the aerosols converges towards depending solely on the initial diameter and the height of the infectious individual’s mouth and nose from the ground. Additionally, the time to evaporate to the equilibrium diameter increases with increasing size. And from a practical standpoint, it is necessary in order to keep M_c from getting too large since $$ M_c\sim \mathcal{O}\left(d_M^3\right)$$ for sufficiently large d_M and pathogen concentration in the infectious individual’s respiratory tract fluid ρ_p. If we assume that the aerosols are approximately spherical (reasonably true except potentially when completely dried out) and their density is approximately equal to that of water ρ_w, the aerosols’ inertial response times τ_p to fluid motions from Stokes drag (we are assuming they are small enough that contributions beyond Stokes drag are negligible) and gravitational settling terminal velocity u_g are

Both grow quadratically with diameter, which does not lend itself to a well defined cutoff scale. And additionally one must consider that once exhaled, the aerosols will tend to evaporate (relative humidity in the environment is typically lower than in the respiratory tract where it is close to 100%) thereby reducing their inertia and terminal velocities. For 10 μm, 20 μm, and 50 μm diameter aerosols; the terminal velocities at 20°C and atmospheric pressure are 3.0 mm s^-1, 1.2 cm s^-1, and 7.5 cm s^-1 respectively. However, larger aerosols take longer to evaporate/grow to their equilibrium diameter and therefore will settle at a faster rate initially than their final equilibrium diameter suggests, which makes them even more likely to be lost due to settling than smaller aerosols.

The simulations of Chong et al. [21] indicate that 100 μm aerosols are quite ballistic and quickly fall out of the exhaled plume, but 10 μm aerosols are carried along with the plume and stay in the air despite their evaporation being greatly slowed. This suggests that d_M should be chosen somewhere in the 10–100 μm range, which is further supported by the Wells curves found by Xie et al. [25]. For lack of a better suggestion; we suggest the use of d_M = 50 μm, which will be explored in the Discussion. Before evaporating, the terminal velocity is 7.5 cm s^-1. If the evaporation ratio is a typical value in the $$ \frac{1}{2}-\frac{1}{5}$$ range, the final evaporated diameter would be in the 10–25 μm range and have terminal velocities in the 3–19 mm s^-1 range which is still in the range that indoor environment air flow can keep suspended (though with a high loss rate).

People and infectious aerosol production

We will denote infectious individuals by the subscript I, susceptible individuals by the subscript S, and other individuals by the subscript O. The Other category is all the individuals who are non-infectious non-susceptible. This includes individuals that are immune before they enter the environment (following Jimenez [10]), all of the Removed group in SIR and SEIR models except for the individuals who died or leave the environment, and all of the Exposed group in SEIR models. If one wants to make a full SEIR model from the model presented in this manuscript, the two subgroups (Exposed, and the part of Removed that is still within the environment and breathing plus the previously immune individuals) within this group will have to be treated explicitly. Let the number of individuals in category C be N_C. The total number of individuals is N = N_I + N_S + N_O. The subscript A will be used to refer to all individuals in all categories. Each count is potentially a function of time as individuals can come in and out of the environment. Let 〈⋅〉_C denote taking the average over all individuals in category C.

Let λ_C,j(t) be the volumetric breathing rate of the j’th person in category C. Let E_C,m,in,j(d) and E_C,m,out,j(d) be the filtering efficiency of the mask (if any) of the j’th person in category C for inhalation and exhalation respectively.

The filter efficiencies of most masks vary significantly with aerosol diameter. Note that it is important that the leak rate of the mask be included in its filtering efficiency. These two filtering efficiencies are generally not equal because masks tend to leak more during exhalation than inhalation and aerosols have higher velocities on exhalation than inhalation. We will assume that all infectious aerosols caught by the mask aren’t later re-aerosolized.

Let E_C,r,j(d₀, w, λ_C,j) be the filtering/absorption efficiency of the respiratory tract of the j’th person in category C. This term is non-zero, but it is also not equal to one since the respiratory tract does not absorb all infectious aerosols that pass through it [5, 7, 9, 12]. The best example of this is the observation that individuals can inhale smoke (which is composed of many aerosols) and then exhale some of it back out. The filtering efficiency depends on the original diameter of the aerosols and the evaporation ratio in the environment since d₀ and w give both the initial diameter on inhalation (d ≈ d_e = wd₀), the diameter the aerosols grow towards (d₀) if they are wet on inhalation or nucleate inside the respiratory tract if they are completely dry on inhalation, as well as the time they spend inside the respiratory tract which is inversely proportional to λ_C,j. It must capture the time it takes for the aerosols to nucleate and grow if they are dried out, the growth process inside the respiratory tract, and the absorption probability as they pass through the respiratory tract. A useful reference for the nucleation and the growth processes would be Pruppacher & Klett [27], and a useful reference for the absorption processes for particles in the respiratory tract would be ICRP [39].

The diameter will be d_e = wd when passing through the mask on inhalation, and d₀ when passing through the mask on exhalation since the humidity between the mouth and nose and the mask is high and the distance is short, so there is little time for evaporation. It is often easier to work with the survival efficiencies rather than the filtering efficiencies, defined as

We will assume that the number of infectious pathogen copies in each exhaled droplet/aerosol follow a Poisson distribution where the mean count is equal to the droplet/aerosol’s initial volume times the pathogen load in respiratory tract fluid at the point of production. This excludes diseases where pathogenic agents stick together and clump. Note that this implicitly means we are assuming that the pathogen volume fraction in the respiratory tract fluid is small. Otherwise, the non-Poissonity caused by there being a maximum number of pathogen copies that can fit in a finite sized drop will NOT be negligible.

Let ρ_j(d₀, t)dd₀ be the number density in exhaled air of the aerosols with diameters between d₀ and d₀ + dd₀ exhaled by the j’th infectious individual at time t. Let ρ_p,j(t) be the pathogen concentration in the j’th person’s respiratory tract fluid where the aerosols are being produced. The mean/expected multiplicity for infectious aerosols produced by the j’th infectious individual for any d₀ is

If the pathogen copies are Poisson distributed in the fluid that makes up the aerosols (no clumping, etc.), then

Sources

We will denote sources by the symbol β with a subscript denoting the individual source. All of them are normalized by the volume of the environment, V.

First, ventilation with other rooms brings infectious aerosols inside at a rate, normalized by the environment volume, of

The other source is the infectious individuals exhaling aerosols with pathogen copies in them. The total production from the infectious individuals normalized by the environment volume is the sum of the products of the breathing rate, the exhaled aerosol concentration density, and the survival efficiency of the mask [7, 10]; which is

Other than not replacing the average of the product with the product of the averages, following aerosols by multiplicity and diamater, and not using quanta; this term is identical to the equivalent term by Nazaroff, Nicas & Miller [7] and Jimenez [10] and, if masks are removed, that of the original formulation [4].

Now, it may be the case that an infectious person has different respiratory tract pathogen concentrations at different locations where exhaled aerosols are produced (e.g. different concentrations in the lungs and mouth). In this case, one would split the term in Eq (37) for the particular infectious person into separate terms for each location of production and use different ρ_j(d₀, t) and 〈k〉(d₀, t)_j in n_I,j,k(d₀, t) from Eq (35).

Sinks

Sinks are proportional to the concentration density n_k. We will denote all sinks divided the concentration density by the symbol α with a subscript denoting the individual source. All of them are normalized by the volume of the environment, V. Unlike the sources, none of the sinks (except inactivation, considered separately) depend on the multiplicity and therefore the subscript k is dropped. Note that inactivation is treated separately later since it is a flux term when considering each multiplicity separately, unlike in the traditional formulation where it is a sink.

The volume normalized loss rate coefficients of infectious aerosols due to exchange of clean air with outdoors and other rooms are just the volume normalized flow rates [9, 33] and are

Let E_v(d) be the filtering efficiency of the circulating ventilation system for aerosols with diameter d. The diameter when an aerosol reaches this filter is d ≈ d_e = w(d₀, t)d₀. Then the volume normalized loss rate coefficient from the circulating ventilation system [4] is

Aerosols also gravitationally settle and deposit onto surfaces. We will treat these processes as simple loss rates proportional to their concentration densities just as one does with radioactive decay. The volume normalized loss rates divided by the concentration density, of gravitational settling and deposition are defined to be α_g(w(d₀, t)d₀) and α_d(w(d₀, t)d₀) respectively; which depend on the room geometry, aerosol diameter, and air flow in the room. A possible approximate expression for the settling loss term [9] would be

Sinks from individuals inhaling aerosols

Unfortunately, when individuals inhale infectious aerosols, some are absorbed thereby causing a risk of infection. While this phenomena is not desired for susceptible individuals, we must consider the loss rate from this process by the susceptible individuals as well as the infectious individuals and the non-infectious non-susceptible individuals. There are three steps to the filtering process for the j’th person of category C: passing through the mask on inhalation, passing through the respiratory tract, and then passing through the mask on exhalation.

The total survival probability of an aerosol going through all three steps is the product of the individual survival rates. The total filtering efficiency is then one minus the total survival rate. But, there is a time delay between when the aerosols are removed from the environment on inhalation and when the survivors are exhaled back out. As long as this time is short compared to all other time scales such as mixing times in the room, the time scales of all other sinks, the time scale of inactivation, etc.; we can ignore this time delay and consider the re-exhalation to occur at the same time. This assumption implies that we can neglect possible changes in multiplicity by inactivation while the aerosols are in the respiratory tract. In most situations, this is a reasonably good assumption. But, at a swimming pool where people regularly hold their breath for long periods of time, this assumption could be violated for the highest multiplicities since the inactivation rate from k to k − 1 is proportional to k.

We assume that the individuals are far enough away from sources and that the environment is well-mixed enough that the concentration density in the air inhaled by each individual is approximately the average concentration density n_k(d₀, t). See the Discussion for a brief qualitative discussion of what the required corrections would look like when this assumption is not valid. Note that we will make the assumption that the self-proximity correction for infectious individuals is negligible (each infectious individual is by definition in close proximity to an infectious individual, themself), though this could pose an issue when the transport of infectious aerosols in the environment to an individual is weak [29]. Then the number of aerosols that are inhaled by a person is equal to λ_C,j(t)n_k(d₀, t). The volume normalized sink coefficient from this filtering is then

Flux: Inactivation

When a pathogen in an aerosol with multiplicity k inactivates, the aerosol’s multiplicity changes to k − 1. We will model inactivation of pathogen copies as exponential decay with inactivation rate γ(t), which might depend on time (e.g. dependence on UV light intensity, humidity, etc. that could be fluctuating in time). For aerosols with a multiplicity of k, the volume normalized loss rate to multiplicity k − 1 is just

Two pathogen copies will never inactivate at exactly the same time; so we don’t have to consider flux terms beyond the two neighboring multiplicities.

General concentration density equations

All of the sources, sinks, and flux terms can be collected to make the system of differential equations describing the infectious aerosol concentration density, which is

We have assumed that shattering and collisional coalescence of infectious aerosols, whether from turbulent induced collisions or differential gravitational settling, is negligible. Collisional coalescence could begin to be important if there are a significant number of very large aerosols and/or n_k is very large. Particularly, d > 100 μm aerosols/droplets, even though they will generally settle to the ground/floor before evaporating to their equilibrium diameter [24, 25], can capture smaller aerosols on their way to the ground/floor [26, 27, 38]. This will generally be negligible unless individuals are situated in the environment such that the large aerosols exhaled by one person (who need not be infectious) will fall through the exhaled aerosol plume of an infectious individual, and potentially negligible even then. If the aerosol concentration, including non-infectious aerosols, reach the levels seen in atmospheric clouds, collisional coalescence might also have to be considered along with keeping track of k = 0 aerosols; though this is very unlikely in indoor environments except when there is a lot of smoke or artificial fog machines are in use, like in a discotheque or theater.

Then, putting the flux terms into Eq (44), we have the following system of ODEs to get the concentration density

Luckily this is a system of ODEs rather than PDEs with flux terms in diameter (involving derivatives with respect to diameters). This is the advantage of choosing d₀ or d_D instead of d_e. For practical applications, this also means that we can also split the diameter range into bins and solve it for each bin separately since there are no flux terms between bins. (See S3 Appendix for how to bin the model with respect to diameter.).

This is a linear inhomogeneous finite system of coupled ODEs at each d₀. The number of equations in the system is finite since k is non-negative and there is the maximum theoretical multiplicity M. Moreover, we don’t even need to care about k = 0 since those aerosols are no longer an infection hazard. Additionally, the system that needs to be solved is smaller if M_c < M. If M_c = 1, then we have only one ODE. This situation occurs if the pathogen load of respiratory tract fluid is low enough that very few aerosols have 2 or more pathogen copies in them.

Note that this model demonstrates superposition with respect to sources since it is linear, as expected intuitively—each aerosol is independent of all others, therefore the response (concentration density and expected dose) from each individual source is independent of all other sources. If n_k,1 and n_k,2 are solutions for the same α and γ but different sources β_k,1 and β_k,2 respectively, then the solution for β_k = β_k,1 + β_k,2 is n_k = n_k,1 + n_k,2.

Infection risk

Let μ_j,k be the average number of aerosols with multiplicity k absorbed by the j’th susceptible individual from time t₀ to time t. At any particular instant of time, the average number of such aerosols of each original diameter d₀ entering the person’s mask if they are wearing a mask or their mouth and nose if they aren’t is λ_S,j(t)n_k(d₀, t). Note that we have assumed that the j’th susceptible individual is not close enough to any sources or filtering sinks that the concentration density of the air they are inhaling deviates significantly from n_k(d₀, t). For susceptible individuals in close proximity to infectious individuals, close to the output of ventilation, etc.; corrections must be applied. See the Discussion for a qualitative discussion on what the required corrections would look like.

A fraction S_S,m,in,j(d₀, t) will survive the mask to enter the respiratory tract [5–8, 10, 12]. A fraction E_S,r,j(d₀) of those survivors will be absorbed by the respiratory tract [5, 7, 9, 12], which contributes to the dose. The expected average aerosol dose is then the double integral of this over the d₀ and the time between t₀ and t, which is

In order to use the μ_j,k in the multiplicity-corrected dose-response model for the particular disease of interest $$ \mathfrak{R}$$, we need to first assume that the aerosol dose for each multiplicity follows a Poisson distribution with μ_j,k as the means and that each is independent of each other (no correlations). This requires the well-mixed assumption like many other parts of the model.

But it also requires that the effect of turbulent inertial clustering is negligible. We will now show that it is negligible except possibly at extremely high aerosol concentrations. It will be negligible if the aerosol Stokes numbers St = τ_p/τ_η are very small (St ≪ 1) [41, 42] where τ_p is the aerosol inertial response time scale from Eq (29) and τ_η is the Kolmogorov time scale of the turbulence in the environment, which is $$ {\tau }_{\eta }=\sqrt{{\nu }_a/ϵ}$$ where ϵ is the turbulent dissipation rate. It will also be small if the typical inter-aerosol distance $$ {\overline{d}}_a\sim {\mathcal{N}}^{-1/3}$$, where $$ \mathcal{N}$$ is the total infectious aerosol concentration for all d₀ and k, is much larger than the typical scale of turbulent inertial clustering (i.e. the fraction of aerosols with a neighbor in the clustering range is low). The typical scale of turbulent inertial clustering is about 10η [41, 42] where $$ \eta ={({\nu }_a^3/ϵ)}^{1/4}$$ is the Kolmogorov length scale of the turbulence. This means that as long as St ≪ 1 and/or $$ {\mathcal{N}}^{-1/3}\gg 10\eta $$, the deviations of the aerosol doses from independent Poisson distributions will be negligible. The situation will be worst for the largest w(d_M, t)d_M sized aerosols in high enough humidity that w(d_M, t) ≈ 1. For a low dissipation rate of ϵ = 1 mW kg^-1; St = 0.06 for a d_M sized aerosol and the number density limit is $$ \mathcal{N}⪡4\times {10}^5$$ m^-3. The Stokes number is small, so the turbulent inertial clustering’s effect will be small even if $$ \mathcal{N}$$ exceeded that limit. For a higher dissipation rate of ϵ = 1 W kg^-1; St = 2.0 for a d_M sized aerosol and the number density limit is $$ \mathcal{N}⪡7\times {10}^7$$ m^-3. While the Stokes number is large, the number density limit is very high so turbulent inertial clustering’s effect will generally be small. For a high for indoors dissipation rate of ϵ = 10 W kg^-1; St = 6.3 for a d_M sized aerosol and the number density limit is $$ \mathcal{N}⪡4\times {10}^8$$ m^-3. While the Stokes number is large, the number density limit is very high so turbulent inertial clustering’s effect will generally be small. Thus, turbulent inertial clustering will have a negligible effect on the Poissonity and independence of the aerosol dose distributions except possibly at extraordinarily high aerosol concentrations.

General

There is an analytical solution to Eq (45), though it is not closed form unless the time dependence of α, β, and γ allow it. Eq (45) can be rewritten in matrix-vector form as

The general solution in matrix-vector form, shown in S1 Appendix, is

Working this out using the structure of the diagonalization of A in S1 Appendix, the general solution for each k is

Coefficients constant in time

We cannot go further in simplifying the general solution from Eq (50) without knowing the time dependence of α, $$ \overrightarrow{\beta }$$, and γ. In many situations; α, $$ \overrightarrow{\beta }$$, and γ are approximately constant with respect to time. If this is so; the general solution from Eq (50) and its time integral from t₀ to t (needed for the dose) for the trivial case that γ = 0 but α ≠ 0 is

It is possible for λ_S,j to be a function of t but α not be (i.e. there is cancelation). But if λ_S,j and w are constant, the expected average aerosol dose of multiplicity k for the j’th susceptible individual in Eq (46) becomes

Calculation of $$ {\overrightarrow{n}}_k(d_0,t)$$, $$ {\overrightarrow{n}}_{\infty ,k}\left(d_0\right)$$, $$ {\int }_{t_0}^t{n}_k(d_0,v)dv$$ scales as $$ \mathcal{O}\left(M_c^3\right)$$ due to there being M_c multiplicities and double sums in U_k and W_k that scale as M_c. There is a recursive solution for $$ {\overrightarrow{n}}_{\infty ,k}\left(d_0\right)$$ which is linear in M_c, and recursive solutions for all the U_k and W_k which are quadratic in M_c. Additionally, the recursive formulas don’t require as much numerical precision in the intermediate steps to get a desired final precision as shown in S5 Appendix. From S1 Appendix, the recursive solutions start at k = M_c and proceed downwards to k = 1. They are

This recursive analytical solution for $$ \overrightarrow{n}$$ is checked against a numerical solution of Eq (47) for a simple case and a very small time step in S2 Appendix. The relative differences for the simple case are very small at less than 10⁻¹². See S5 Appendix for numerical considerations for evaluating the analytical solutions on a computer or solving Eq (47) with a numerical ODE solver. The number of terms for both are discussed, as well as the required precision and maximum magnitude required for floating point numbers used to calculate the analytical solution formulas.

Room, people, and filter efficiencies

We consider a hypothetical example based on the ongoing SARS-CoV-2 pandemic—a poorly ventilated seminar room with two infectious individuals with SARS-CoV-2 at the very upper end of viral concentrations (viral load) and one of them continuously coughing. We assume that the room is well-mixed and that the individuals are far enough apart from each other and the ventilation that no corrections to n_k(d₀, t) need to be applied at any source or sink, nor in the calculated absorbed doses. Let the room have volume V = 200 m^-3 with a height of h = 4 m, with ventilation q_r = 0, q_v = 0, and q_o = 0.5 hr^-1. We will ignore surface tension’s effects on w. Let the humidity be such that the evaporation ratio is $$ w=\frac{1}{3}$$, which is a constant with respect to both t and d₀. We ignore deposition (α_d = 0). Let there be N_S = 15 susceptible individuals in groups of 5 wearing no mask, a simple1 mask, and a simple2 mask (defined later); and no non-infectious non-susceptible individuals (N_O = 0). The susceptible individuals will be assumed to be sedentary/passive adults with a breathing rate of λ_S,j = 0.3 m^-3 hr^-1, which is in the range of mean breathing rates for this activity from the U.S. EPA’s Exposure Factors Handbook Table 6.2 [14]. The pathogen concentration for SARS-CoV-2 varies widely across individuals, location in the body, and stage of the disease [35, 36, 43, 44], and can sometimes get as high as the 10¹⁰–10¹¹ cm^-3 range [35, 36]. We will use this upper range because it makes the model more challenging to solve due to the larger M_c and due to the interest in so called “super-spreading events”. The situation is composed of two stages (Stages 1 and 2) that each start when an infectious individual enters the room. Initially, there are no infectious aerosols in the room, meaning n_0,k(d₀) = 0. Stage 1; at t = t₀ = 0, one infectious individual enters the room who is speaking, wearing no mask, breathing at a rate λ_I,j = 0.5 m^-3 hr^-1 (just below an 0.54 m^-3 hr^-1 average value for reading out loud [15]), and has a high respiratory tract fluid pathogen concentration of ρ_p,j = 10¹⁰ cm^-3. Stage 2; then at t = 3 hr, one more infectious individual enters the room who is continuously coughing while wearing a simple2 mask, breathing at a higher rate of λ_I,j = 2.0 m^-3 hr^-1, and has a higher respiratory tract fluid pathogen concentration of ρ_p,j = 10¹¹ cm^-3 at the very upper range for SARS-CoV-2. We chose this estimated continuous coughing breathing rate by deducing a breathing rate range from Hegland, Troche & Davenport [17] for continuous 3 cough cycles (heavily using their Fig 1), getting a breathing rate range of 1.9–2.3 m^-3 hr^-1 from which we chose 2.0 m^-3 hr^-1.

We use mask filter efficiencies of the functional form

none (no mask) E₀ = E_∞ = 0.

mask simple1 E₀ = 0.2 and E_∞ = 0.8.

mask simple2 E₀ = 0.95 and E_∞ = 0.99.

The filtering efficiencies of both the simple1 and simple2 masks are shown in S1 Fig. The mask parameters were chosen such that they are more efficient at filtering large aerosols/droplets than small ones, with the simple2 mask being better than the simple1 mask. The simple1 and simple2 masks could reasonably correspond to a reasonably well fitted home-made cloth mask and an excellently fitted FFP2 mask, though here we have treated their leak rate to be the same during inhalation as exhalation (not true with most real masks). At the largest sizes, leakage doesn’t matter as much since the aerosols are more ballistic. Let us assume that $$ E_{C,r,j}(d_0,w(d_0,t),{\lambda }_{C,j}\left(t\right))=\frac{1}{2}$$ for everyone.

Disease and infectious aerosol production

We assume that an exponential-dose response model is the correct model to use for SARS-CoV-2 since the exponential model works better than the beta-Poisson model for two other human infecting corona viruses (SARS-CoV-1 and HCoV-229E) [34]. In absence of a good value to use for r, we use the same value of r as found for SARS-CoV-1 in mice which is r = 2.45 × 10⁻³ and the same value of r as found for HCoV-229E in humans which is r = 5.39 × 10⁻² [34]. We use γ = 0.64 hr^-1 as the inactivation rate for SARS-CoV-2 [45].

We approximate the SARS-CoV-2 pathogen as a sphere with a diameter of 100 nm, which is close to the correct size and the rough shape with the surface proteins removed (actually an ellipsoid) [37]. We use the aerosol size distributions for speaking and coughing from Johnson et al. [22], but extrapolate them to smaller diameters (from 800 nm to 100 nm). This is used with Eqs (26) and (27) to get the β_I,k. They are shown in the top-right panel of Fig 3. The aerosol size distributions have two peaks at approximately 2 μm and 100 μm. This puts d_M between the trough (between the two peaks) and the second larger diameter peak.

Fig 3

Model solution for example.

Solution to the example case. (Top-Left) The total pathogen and infectious aerosol concentrations over time. (Top-Right) The infectious aerosol concentration densities in the room as a function of d₀ at t = 6 hr compared to the aerosol concentration densities being exhaled by speaking and coughing individuals from Johnson et al. [22] scaled by 10⁻⁴ to make them have comparable magnitudes. (Bottom-Left, Bottom-Right) The mean infection risk $$ {\mathfrak{R}}_E$$ for the susceptible individuals based on the mask they are wearing (none, simple1, or simple2) using (Bottom-Left) r = 2.45 × 10⁻³ (Bottom-Right) r = 5.39 × 10⁻².

Concentration densities and infection risk

We now find the infectious aerosol concentration densities and doses, and mean infection risks $$ {\mathfrak{R}}_E$$. First, we split the diameter range between d_m,1 = 0.1 μm and d_M = 50 μm into 20 logarithmically spaced bins; and determine the bin average values for the coefficients over each bin by integration following the scheme in S3 Appendix. The infectious individuals source parameters for the i’th bin, β_I,k|_i, are calculated numerically via Simpson’s rule for integration with 1000 equal linear width sub-bins in each bin. The particular choice of the mask survival efficiency in Eq (81) and w being constant lets the other binning integrals be calculated analytically.

The model is solved for Stage 1 and then the final values used as initial values for Stage 2 because this makes it so that α and β_k are constant in time when solving the model (all changes are between stages). For M_c, we used the maximum value of M_c,I,j for each infectious individual present at each Stage with T = 10⁻³. Note that M_c stayed the same or increased for each bin going from Stage 1 to Stage 2 with the addition of one more infectious individual.

For the i’th bin, the n_k|_i(t) and μ_j,k|_i(t) are solved analytically if M_c ≤ 500 using the recursive solution and numerically if M_c > 500, both in IEEE-754 binary64 floating point (also known as double precision and float64). This threshold between analytical and numerical solving was chosen to use the analytical solution as much as possible without overflow in V_k (see S5 Appendix). As shown in S5 Appendix, binary64 numbers provide sufficient precision and allowed maximum magnitude. Note that overflow is easy to spot as infinities, which were not seen so this number format was sufficient to prevent overflow. When doing it numerically, Eq (47) along with $$ {\int }_0^{t}n(d_0,v)dv$$ were solved using Runge-Kutta 4 with a time step of 10⁻⁴ hr, which is required for stability and an accurate solution with the large α|_i+ M_c γ values in the largest bin. After determining α|_i and $$ \overrightarrow{\beta }{|}_i$$, the solutions were calculated with the help of the PMADRA (Poly-Multiplicity Airborne Disease Risk Assessment) software suite we wrote for the purpose (https://gitlab.gwdg.de/mpids-lfpn-public/pmadra), specifically the Python 3.5 or newer implementation pypmadra version 0.2.1 (https://gitlab.gwdg.de/mpids-lfpn-public/pmadra/pypmadra) using the Fortran 2008 accelerator library libpmadra version 0.2.1 (https://gitlab.gwdg.de/mpids-lfpn-public/pmadra/libpmadra). The main results are shown in Fig 3.

The total pathogen concentration is slightly less than double the infectious aerosol concentration in Stage 1, and slightly higher than double in Stage 2. This means that the average multiplicity in both stages is approximately two, and it increases slightly from Stage 1 to Stage 2 which is expected with the higher viral load in the second infectious individual. Also, as expected, increasing r (infection risk of each individual pathogen) increases the infection risk. As expected, susceptible individuals wearing masks decrease their infection risk and increasing exposure increases their infection risk.

Comparing the infectious aerosol concentration density in the room with the aerosol concentration densities exhaled by the infectious individuals as a function of d₀ (see top-right panel of Fig 3); we can see how as d₀ increases, the probability of an aerosol being infectious increases (infectious aerosol concentration density decreases slower after the first peak than the exhaled aerosol concentration densities) but at the largest d₀ > 15 μm the increasing α due to stronger gravitational settling causes the infectious aerosol concentration density to grow slower after the trough than the exhaled aerosol concentration densities from the infectious individuals (including the speaking individual who is not wearing a mask). To see the latter, the strengths of the sinks α and total sinks α+ kγ are shown in Fig 4 and we can see that settling causes α to increase by over a factor of 10 from 100 nm to 50 nm. Fig 4 additionally shows the increase in the total sink strength for the largest multiplicities M_c being considered due to inactivation. The large difference between the total sink strength between k = M_c and k = 1 makes the system of ODEs stiff.

Fig 4

Sink strength by bin.

The strength of the sink terms for each bin with 80 bins, which is α without inactivation, α + γ for k = 1, and α + M_c γ for k = M_c (different values for Stage 1 and 2).

The pathogen concentrations as a function of d₀ and k right after the beginning and at the end of each Stage are shown in S2 Fig. For large diameters, the concentrations at the beginning of each Stage are initially in a narrow band around the expected multiplicity in each diameter bin but by the end of each Stage the distributions have widened downward as inactivation fills in the lower multiplicities.

The results of choosing different numbers of bins (5, 20, and 80) is shown in S3 Fig. The difference in the concentration densities between 5 bins and 20 bins is substantial, but the difference between 20 and 80 is small. This means that in our example; for concentration densities, 20 bins is sufficient to capture the variation in α(d₀) and β_k(d₀) with respect to diameter, but 5 is too few and 80 is a lot more effort for little gain. But for the $$ {\mathfrak{R}}_E$$, the difference between the solutions for different number of bins is very small for the smaller r = 2.45 × 10⁻³, but more noticeable but still small for the larger r = 5.39 × 10⁻².

Effect of multiplicity on dose-response

We consider a few hypothetical examples to ellucidate the importance of multiplicity in the dose-response using the corrected exponential model in Eq (12). Another dose-response model could be chosen and the resulting values would differ, but the general pattern would be the same.

First, let’s reconsider the example case but with all pathogen production forced to be mono-multiplicity. We set the new $$ {\beta }_{1,new}={\sum }_{k=1}^{M_c}k{\beta }_k$$ and all other β_k,new = 0 ∀ k ≠ 1 and then set M_c = 1 for all bins. This is equivalent to going to each bin, taking the total aerosol volume production, finding the expected number of pathogen copies in that volume, and redistributing the volume so that each pathogen is alone in an aerosol but not changing d₀ anywhere. Or put equivalently, making Eq (47) track pathogen copies instead of aerosols and ignoring multiplicity. To quantify the difference, we took a simplified version of the example where the second coughing infectious individual was removed, the ρ_p of the first speaking infectious individual was adjusted, and we took the steady state case where $$ {\overrightarrow{n}}_0={\overrightarrow{n}}_{\infty }$$ and calculated the constant dμ_j,k/dt for each susceptible individual. Then using the constant dμ_j,k/dt and an initial dose of zero, we found the time, τ₅₀, required for $$ {\mathfrak{R}}_E$$ to be 50% (note that the particular choice does not matter, the curve is identical for any chosen risk). This was calculated for the 80 diameter bins example to keep errors from finite bin width small, and a range of r values up to the maximum value r = 1. Ignoring multiplicity causes τ₅₀ to be underestimated (overestimation of risk). The underestimate of τ₅₀ is shown in Fig 5.

Fig 5

Effect of ignoring multiplicity, full version.

Full version of Fig 1 with more ρ_p and the effect of masks. Plot of the ratio of the time required to reach a 50% infection risk when multiplicity is ignored τ_50,ignore to when it is fully accounted for τ_50,full for different respiratory tract fluid pathogen concentrations ρ_p. We are considering the same situation as in the worked example, but at steady-state with just the speaking mask-less infectious individual and the risk to a susceptible individual whose exposure starts after steady state is reached. The ratio is shown for different combinations of mask on the susceptible individual (none and simple2) and for different r. The legend lists the r, mask combinations in the same order as the lines from top to bottom. We assumed a 100 nm diameter spherical pathogen and used 80 diameter bins and chose the M_c (maximum multiplicity considered) heuristic threshold to be T = 0.01 (include 99% of pathogen production).

The underestimation increases with increasing ρ_p and r, and decreases when wearing a mask that is more efficient at filtering large aerosols than small aerosols. The largest aerosols have the greatest multiplicities, which means that a mask that filters them out better than small aerosols reduces the effect of ignoring multiplicity. As ρ_p increases, the expected multiplicity range for each d₀ increases which makes ignoring multiplicity underestimate τ₅₀ more. For the r values considered here, ρ_p ≤ 10⁹ cm^-3 underestimates τ₅₀ by at most 20% and ρ_p ≤ 10⁸ cm^-3 underestimates it by at most 12%. But for ρ_p = 10¹¹ cm^-3, the underestimation is up to 67%. To better understand these patterns, we need to consider two more hypothetical situations.

Let the average pathogen dose be 〈Δ〉 = r⁻¹ and all infectious aerosols have the exact same multiplicity k. Then, the μ for all other multiplicities is zero and μ_k = 〈Δ〉/k. Essentially, we are dividing the same number of pathogen copies among fewer and fewer aerosols as we increase the number of pathogen copies in each one. The mean infection risk for this constant average dose is shown on the left side of Fig 6 as a function of k for four different r. As the multiplicity increases, the mean infection risk decreases even though the average dose is the same. For k ≪ r⁻¹, the effect of multiplicity on $$ {\mathfrak{R}}_E$$ is small. It starts to rapidly decrease near k ∼ r⁻¹ and converges towards zero, because the number of pathogen copies in each aerosol is large enough that each aerosol has a high probability of causing infection by itself but the aerosols are decreasing in number faster than the risk can increase. The risk per aerosol can’t exceed 100% no matter how many pathogen copies are in an aerosol.

Fig 6

Multiplicity’s impact on infection risk.

Plots of mean infection risk ($$ {\mathfrak{R}}_E$$) using the modified exponential dose-response model when all infectious aerosols have the same number of pathogen copies in them k. (Left) The infection risk as a function of k for fixed average dose 〈Δ〉 = r⁻¹ for different single pathogen infection probabilities r. (Right) The infection risk as a function of the dose scaled by r (〈Δ〉r) for different k and the same fixed r = 10⁻² (r⁻¹ = 100).

Another way to see this is to consider another hypothetical. Let’s consider the mean infection risk if all aerosols have multiplicity k as we vary r〈Δ〉 for fixed r. This is shown on the right side of Fig 6 for r = 10⁻². For low k ≪ r⁻¹, the infection risk curves are nearly identical. For k ≥ r⁻¹, the infection risk decreases for increasing k.

Overall, this means that if the typical infectious aerosol multiplicity is on the order of or greater than r⁻¹, there can be a significant decrease in the infection probability for the same average dose. This has implications for large aerosols when the respiratory tract fluid pathogen concentration ρ_p,j is large. Large aerosols where 〈k〉 ≳ r⁻¹ will contribute less to the infection risk than would otherwise be expected from their resulting average pathogen dose 〈Δ_k〉. While we must have M_c > 〈k〉, M_c is usable as a proxy for which diameters the multiplicity causes a substantial correction to the dose-response. If we were to consider r = 2.45 × 10⁻³ as was done in the example, Fig 2 shows that this would be important for d₀ > 15 μm for a high viral concentration of ρ_p,j = 10¹¹ cm^-3 and d₀ > 30 for the lower but still high viral concentration of ρ_p,j = 10¹⁰ cm^-3. If we were to consider r = 5.39 × 10⁻² as was also done in the example, Fig 2 shows that this would be important for d₀ > 5 μm for a high viral concentration of ρ_p,j = 10¹¹ cm^-3 and d₀ > 10 for the lower but still high viral concentration of ρ_p,j = 10¹⁰ cm^-3.

Going back to the risk overestimation from ignoring multiplicity in Fig 5, decreasing r decreases the underestimation in τ₅₀ because the ratio of the average multiplicity in the larger diameter bins to r⁻¹ is smaller. A mask that filters large aerosols better than small aerosols reduces the effect of ignoring multiplicity because larger aerosols have higher multiplicities.

Filtering by the people

We introduced the sink terms α_C,f for filtering by the individuals in the environment as they inhale aerosols with many being absorbed by their mask or respiratory tract rather than being exhaled back out into the environment. To determine when this sink matters, we need to consider the total volume of air that is filtered, ignore the filtering efficiencies, and compare it to the ventilation. The volumetric rate of air filtration by the individuals normalized by the volume of the environment is

The mean adult breathing rates from sedentary/passive to high intensity activity ranges between 0.25 m³ hr^-1 and 3.2 m³ hr^-1 [14]. For sitting, it would be hard to get σ_A to be more than 1 m^-2 but it would be possible while standing (some public events) though the well-mixed assumption would be breaking down in either case. For a typical room height of 〈h〉 = 4 m, this density limit would yield max(q_p)∈[0.063, 0.8] hr^-1. If the environment is poorly ventilated (total ventilation rate q_v + q_o + q_r less than 1 hr^-1), this high people density would mean the filtering effect of the people would not be negligible compared to the ventilation. But with even moderate ventilation, the contribution of α_C,f would be negligible unless all the ventilation is circulating ventilation (q_o = q_r = 0) with no filter or a very poor filter. For 1.5 and 2 m social distancing, the maximum σ_A are 0.14 and 0.080 m^-2 respectively. For a typical room height of 〈h〉 = 4 m, this density limit would yield max(q_p) ∈ [0.005, 0.11] hr^-1 which would be negligible in almost all circumstances. For taller rooms, the contribution would be smaller if the total ventilation rate is held constant.

If the fraction of individuals who are infectious is held constant, then N_I ∼ σ_A. Since β_k ∼ N_I and α_C,f ∼ N_C but the non α_C,f terms of α stay constant, the source increases faster than the sinks meaning that n_k increases and therefore $$ \mathfrak{R}$$ increases. So, increasing σ_A with everything else held constant increases the risk for the susceptible individuals. Thus, deliberately making α_C,f non-negligible is not a viable strategy to decrease risk. If the α_C,f dominate over the ventilation, the situation is actually quite hazardous from an infection transmission perspective. It is just that if one ignores the terms, one would overestimate the risk in such a crowded and poorly ventilated space.

Effect of masks

The filtering effects of masks show up in the source β_I,k, the sinks α_C,f, and the total dose over time μ_j,k. Masks can substantially improve the total filtering efficiency of the people in α_C,f since aerosols have to pass through the mask twice, once on inhalation and again on exhalation at a larger diameter (many masks are better at filtering larger diameters than small diameters). But unless the ventilation is poor and there are a lot of people, this increase in α_C,f will have only a small effect on the total sink α. Instead, the main contribution is to reducing β_I,k and μ_j,k which are both linearly proportional to the mask survival efficiency, which can be seen in the example situation.

In the example during Stage 1, there is one infectious individual in the room who is not wearing a mask and the total pathogen concentration reaches about 40 m^-3 after 3 hr (Fig 3). During Stage 2, an addition infectious individual has entered the room. The second infectious individual’s ρ_p is 10 times greater than the first person’s and they are breathing at 4 times the rate; which would mean 40 times the pathogen exhalation rate by itself. Additionally, they are coughing rather than speaking, with the resulting larger exhaled aerosol concentration density ρ_j (top-right panel of Fig 3); which increases the number of exhaled pathogen copies further. But, they are wearing a mask which reduces the number of infectious aerosols that survive to reach the environment by a factor of 20–100 depending on the diameter. Due to this, the total pathogen concentration doesn’t increase by a factor of over 40 but instead approximately triples, reaching approximately 140 m^-3.

The reduction in the average dose μ_j,k and therefore infection risk $$ \mathfrak{R}$$ when susceptible individuals wear masks can also be seen in Fig 3. Even the simple1 mask gives some improvement, and the simple2 mask reduces the infection risk by over an order of magnitude.

Let’s consider the case where all infectious individuals have the same mask survival efficiency and all susceptible individuals have the same mask survival efficiency. If the effects of masks on α is negligible (α_C,f is generally small compared to the other sinks) and β_r,k is negligible; the combined effect of both infectious and susceptible individuals wearing masks on the dose is quadratic in the survival efficiencies, which has shown up in other Wells-Riley formulations in the past [7, 10]. Due to superposition of sources, n_k ∼ S_I,m,out since β_I,k ∼ S_I,m,out. Then, μ_j,k ∼ S_S,m,in n_k ∼ S_S,m,in S_I,m,out, which is a quadratic term. Now α_C,f ∼ S_C,m,in S_C,m,out makes the effect stronger (usually only slightly stronger) than quadratic since it only serves to increase α and therefore decrease n_k further. If everyone wears masks with the exact same survival efficiency S for both inhalation and exhalation that is constant with respect to d₀, then if exposure starts at steady state, μ_j,k ∼ Sn_k,∞ ∼ Sβ_k/α ∼ S²/(1 − cS²) where c ∈ [0, 1) is a constant that depends on the relative importance of the α_C,f in the total α. In this form, it is easier to see how μ_j,k scales super-quadratically in the mask survival efficiency. If just the susceptible or just the infectious individuals wear masks, the reduction drops to being stronger than linear (direct contribution of the mask on reducing β_I,k or reducing μ_j,k plus the effect on α_C,f). If only non-susceptible non-infectious individuals wear masks, there is still a reduction in the dose but it is small since α_O,f is generally small compared to the other sinks, giving a sublinear reduction.

Well-mixed limitation and corrections

The biggest limitation to the model presented here, like all Wells-Riley formulations, is the well-mixed environment assumption. In almost all indoor environments, the assumption breaks down to varying degrees—the infectious aerosol concentration densities at the locations of susceptible individuals and all sinks (except possibly inactivation) depend on their locations in the environment relative to the sources and the air flow. Social distancing helps with this assumption (reduces direct inhalation of undiluted exhaled puffs of aerosols from infectious individuals), but the assumption is still often dubious.

In situations where people, other sources, and localized sinks (or their outputs) are located close to each other; corrections to n_k(d₀, t) must be applied at the location of the individual, other source, or sink. Here, we will qualitatively discuss what simple partial corrections that don’t depend on the history of n_k(d₀, t) would look like. For proximity to the output of filtering sinks, a multiplicative correction would need to be applied with a factor between the sink’s filtering efficiency and one, inclusive, that depends on the location and the properties of the sink such as the flow rate. For proximity to the output of ventilation, the respective filtering efficiency is E_v(w(d₀,t)d₀). For proximity to individuals, the respective filtering efficiency is 1−S_C,m,in,j S_C,r,j,k S_C,m,out,j,k. For proximity to sources, the correction would be to use a weighted average of n_k(d₀, t) and the concentration of the air coming from the source/s with the weights depending on the location and the nature of the source flows and mixing, such as flow rates. For close proximity to ventilation coming from other rooms, this would mean a weighted average with n_r,k(d₀, t) (if there is more than one room, it would be the concentration coming from the room/s whose air is not yet diluted at the location). For close proximity to infectious individuals, this would mean a weighted average with n_I,j,k(d₀, t)[1 − E_I,m,out,j(d₀)]. These partial corrections could be done for specific cases (e.g. susceptible individual 2 is 1 m directly in front of infectious individual 5) or in a statistical way if the pair correlation functions between individuals of each two categories (including in-category) as well as the equivalent correlation functions for relative angles of orientation by distance. More extensive corrections could depend on the history of n_k(d₀, t) and would turn the system of ODEs into a system of Delay Differential Equations (DDEs) or Integro-Differential Equations (IDEs), which would most likely be much harder to solve. At some point, however, it could be easier to do a full fluid and aerosol dynamics treatment.

Any corrections developed for mono-multiplicity Wells-Riley formulations could either be used as is or could be adapted to the poly-multiplicity model presented in this manuscript. Full fluid dynamics simulations with infectious aerosols simulated as passive scalars or as discrete aerosols such as those done by Löhner et al. [40] are the common way to address this limitation entirely and can be used to develop corrections, which are considerably more difficult. Further investigation is needed to find simple approximate ways to generalize the Wells-Riley formulation presented in this manuscript for non-well-mixed environments that are easier than full fluid dynamics with suspended aerosols simulations.

Other model limitations

Another limitation of the model presented here is that it assumes that all infectious aerosols have the same ζ and solute composition, and therefore the same w(d₀, t). This is more easily circumvented in one case. If the solute concentration and composition is constant over time for each individual source (reasonable assumption over small time spans), the model can be solved for each source individually and then the resulting n_k and μ_k,j summed over the individual solutions. This would also be the solution if ζ varies in different locations in the respiratory tract where infectious aerosols are produced for an infectious person. If ζ changes over time for the sources but the solute composition is constant, then one could generalize the model to additionally track ζ (or equivalently d_D) and initial diameter at production d₀ separately.

Another problem is the choice of diameter limits d₀ ∈ [d_m,k, d_M] for each multiplicity. We have neglected the fact that the solute concentration is much greater for d₀ near the lower limit d_m,k as pathogen copies are taking up a very large fraction of the volume and that surface effects may cause additional deviations in the number of pathogen copies in the aerosol from a Poisson distribution. Further work is needed to lift this limitation; though for small pathogens, the total fluid volume and therefore pathogen content in the smallest aerosols where this matters is much less than that of the larger aerosols (see top-right panel of Fig 3) meaning that the effect could be small for small pathogens.

The upper limit d_M is the cutoff where aerosols are so large that they are more ballistic and either settle to the ground before evaporating to equilibrium or still settle too quickly to be mixed even after evaporating to their equilibrium diameter. Based on Xie et al. [25] and Chong et al. [21], we suggested a value d_M = 50 μm. To look at it, we took the example case and re-calculated it for 23 equal log-width bins between 100 μm and 100 μm and considered the concentration densities and mean infection risks if the top 0, 2, and 4 bins were discarded, thereby setting decreasing d_M to 100 μm, 54.8 μm, and 30.1 μm. The time step for the numerical solution had to be reduced to 5 × 10⁻⁶ hr due to the increase in M_c at the larger d_M. This is shown in Fig 7. Increasing d_M increases the total pathogen concentration being tracked since a lot of exhaled respiratory tract fluid volume is contained in the large diameter aerosols, but the total number concentration does not increase much since these big aerosols are few in number. For the larger r = 5.39 × 10⁻², the effect on $$ {\mathfrak{R}}_E$$ is very small as d_M is increased by a factor of approximately three. But for the smaller r = 2.45 × 10⁻³, there is a larger fractional difference in the mean infection risk but the additive difference is no more than 5% for the worst case (no mask). The masks as we have defined them in the example, are better at filtering large particles than small, so they attenuate the effect of increasing d_M on $$ {\mathfrak{R}}_E$$. More investigation is required on this upper diameter limit. Generalizing the model to track d and d₀ and treating evaporation/growth explicitly over time would help alleviate this problem as the high settling rates and the slower evaporation of the largest aerosols could be treated explicitly.

Fig 7

Effect of upper diameter limit d_M.

The example situation was calculated for different values of the upper diameter limit d_M (technically, calculated at the largest and then truncated down as needed). (Left) The total pathogen and infectious aerosol concentration densities over time for each d_M. Note that the differences in the total infectious aerosol concentration density are so small that the lines are right on top of each other. The mean infection risk for each combination of masks on a susceptible individual (none, simple1, simple2) for (Middle) r = 2.45 × 10⁻³ and (Right) r = 5.39 × 10⁻².