Risk assessment for airborne disease transmission by poly-pathogen aerosols

Freja Nordsiek1, Eberhard Bodenschatz1,2,3*, Gholamhossein Bagheri1,

1 Max Planck Institute for Dynamics and Self-Organization (MPIDS), Göttingen,
Niedersachsen, Germany

2 Institute for Dynamics of Complex Systems, University of Göttingen,
Göttingen, Niedersachsen, Germany

3 Laboratory of Atomic and Solid State Physics and Sibley School of Mechanical
and Aerospace Engineering, Cornell University, Ithaca, New York, USA

*<EMAIL_ADDRESS>
## Abstract

In the case of airborne diseases, pathogen copies are transmitted by droplets
of respiratory tract fluid that are exhaled by the infectious that stay
suspended in the air for some time and, after partial or full drying, inhaled
as aerosols by the susceptible. The risk of infection in indoor environments
is typically modelled using the Wells-Riley model or a Wells-Riley-like
formulation, usually assuming the pathogen dose follows a Poisson distribution
(mono-pathogen assumption). Aerosols that hold more than one pathogen copy,
i.e. poly-pathogen aerosols, break this assumption even if the aerosol dose
itself follows a Poisson distribution. For the largest aerosols where the
number of pathogen in each aerosol can sometimes be several hundred or several
thousand, the effect is non-negligible, especially in diseases where the risk
of infection per pathogen is high. Here we report on a generalization of the
Wells-Riley model and dose-response models for poly-pathogen aerosols by
separately modeling each number of pathogen copies per aerosol, while the
aerosol dose itself follows a Poisson distribution. This results in a model
for computational risk assessment suitable for mono-/poly-pathogen aerosols.
We show that the mono-pathogen assumption significantly overestimates the risk
of infection for high pathogen concentrations in the respiratory tract fluid.
The model also includes the aerosol removal due to filtering by the
individuals which becomes significant for poorly ventilated environments with
a high density of individuals, and systematically includes the effects of
facemasks in the infectious aerosol source and sink terms and dose
calculations.

## Introduction

It is well known that some diseases such as influenza, the common cold,
Mycobacterium tuberculosis, measles, and Severe Acute Respiratory Syndrome
Coronavirus 1 (SARS-CoV-1) are airborne; meaning they can be transmitted by
particles (also called liquid droplets, aerosols, or, if completely dried,
droplet nuclei) exhaled by infected individuals that stay suspended in the air
for some time rather than immediately falling to the ground. These particles
come from the fluid of the lungs, vocal chords, mouth, and nose; which
hereafter are all noted as “respiratory tract”. While these particles that
stay airborne as well as larger ones that tend to fall on the ground and
surfaces are all drops/droplets unless they have completely dried out to solid
solute and they are technically aerosols (albeit, sometimes large), the
literature usually refers to small airborne ones as aerosols and the larger
ones that don’t get suspended in the air as drops/droplets, which we shall do
here as well. Note that these diseases can have additional transmission
pathways, which can be more or less significant depending on the
circumstances. Whether Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-
CoV-2) is an airborne disease or not and the relative importance of the
airborne pathway to the pathway of exhaled droplets too large to stay airborne
ballistically getting on susceptible individuals and surfaces have been topics
of ongoing discussion and debate throughout the pandemic [1, 2, 3]. Due to the
possibility that SARS-CoV-2 might be an airborne disease among other
transmission pathways, the SARS-CoV-2 pandemic has brought an increased
interest in airborne disease transmission dynamics and models.

The risk of getting infected from such airborne particles for an individual or
a population has been the subject of numerous studies and analyzes [4, 5, 5,
6, 7, 8, 9, 1, 10]. Many of the transmission mitigation strategies rely on
results obtained by models that take into account a variety of factors to
assess the likelihood of transmission, a good example of which is the World
Health Organization’s 2009 guidelines Natural Ventilation for Infection
Control in Health-Care Settings [11]. Two well-known families of models are
dose-response and Wells-Riley models, which have been extensively used to
model the spread of airborne diseases [12].

There are several dose-response models for various diseases in existence which
consider the risk of infection for an average dose of pathogen copies, taking
full account of the counting statistics [13]. Two common models are the
exponential and beta-Poisson models, which are described in great detail by
Haas, Rose & Gerba [13]. Many diseases follow the exponential model, which has
the added simplicity of having only a single adjustable parameter. Both the
exponential and beta-Poisson models assume that the minimum number of pathogen
copies required for infection, the threshold, is one; but other models exist
for non-unity thresholds. Both models, along with many others, assume that the
number of pathogen copies absorbed follows a Poisson distribution; though
modification of the exponential model for doses following a beta or gamma
distribution has been conducted[5].

The Wells-Riley model, in its original form, takes the steady state balance of
sources and sinks of airborne infectious pathogen copies (in units of quanta)
over a period of time in a well-mixed indoor environment such as a room or
several rooms connected via ventilation (homogeneous concentration assumption)
to calculate the average dose received by susceptible individuals over a time
period, which is then run through an exponential dose-response model [4]. The
original model measures pathogen copies in units of quanta, which is defined
as ID63.21 pathogen copies [12]. Sources such as exhalation by infectious
individuals in the environment and air exchange with other environments with
infectious aerosols and sinks due to fluxes with outside, filtering by the
ventilation, filtering by masks, inactivation, settling, and deposition have
all been considered as well as full temporal modelling of the infectious
aerosol concentration rather than assuming steady-state [4, 5, 5, 6, 7, 8, 9,
1, 10]. At the model’s heart, it is essentially a conservation of infectious
aerosols model, choosing some sources and sinks to explicitly include and
considering others to be negligible, to get the pathogen concentration and
then the average inhaled dose, before using a dose-response model (usually the
exponential model) for the infection risk. Note, in the literature the term
“Wells-Riley model” is sometimes used to refer only to when this formulation
is used with an exponential model, and the terms “Wells-Riley equation” and
“dose-response model” used if other dose-response models are used instead
(e.g. [12]). We will use the term “Wells-Riley formulation” to refer to both.

Wells-Riley formulations are a statistical treatment of airborne disease
transmission. Underneath its source, sink, and respiratory tract absorption
parameters (as well as the choices of which to include and exclude) and its
well-mixed assumption and their caveats/limitations are a mix of fluid
dynamics with inertial particles (aerosols), the biological processes of the
respiratory tract and diseases, thermodynamics, aerosol chemistry, human
behavior and safety interventions (e.g. wearing masks), etc. This includes
breathing rates for different activities [14, 15, 16, 17]; the dynamics of
exhaled puffs and the particles within them by breathing, speech, coughing,
etc. [18, 19, 20, 21]; the generation and ejection of aerosols and larger
droplets by breathing, speech, coughing, etc. [16, 22, 19, 23];
aerosol/droplet growth/evaporation in response to temperature and humidity
[24, 25, 26, 27, 19, 28]; the dynamics of inertial particles in turbulence;
mixing and transport [18, 29, 30, 19, 20, 21]; ventilation and convection in
indoor environments [30]; etc. There have been a number of recent papers that
each go into several of these topics written during the course of the ongoing
SARS-CoV-2 pandemic [31, 29, 30, 19, 28], which while focused on SARS-CoV-2
are also applicable to other airborne diseases. In this manuscript, we will
mostly focus on a statistical treatment.

In the past, various generalizations and improvements have been applied to the
Wells-Riley formulation for situations beyond its original design and to
address its limitations [12]. For example; Nicas, Nazaroff & Hubbard [9]
included sink terms for pathogen inactivation, aerosol settling, and
deposition as well as less than unitary efficiency of the respiratory tract
absorbing infectious aerosols. Wells-Riley formulations have also been
combined with SIR (Susceptible-Infectious-Removed) and SEIR (Susceptible-
Exposed-Infectious-Removed) models [6, 32]. Noakes & Sleigh [33] made a
stochastic model with compartmentalization of the environment into well-mixed
subregions that have less mixing with other regions that can work for periods
of time longer than the incubation period. Recent Wells-Riley based analyzes
during the ongoing SARS-CoV-2 pandemic also include the effects of masks (such
as [10]) unless they are investigating scenarios in which individuals are not
wearing any mask [1], though including the effect of masks predates the
pandemic by decades [5, 6, 7, 8].

One of the biggest assumptions of the Wells-Riley formulation is that the
indoor environment is sufficiently well-mixed [4, 7, 8, 9, 33, 12, 1, 10].
Essentially, it assumes that the infectious aerosol concentration is
homogeneous enough that the concentration inhaled by susceptible individuals
and at all sinks is approximately equal to the volume average concentration
[4, 7, 8, 9, 33, 12, 1, 10]. In reality, there can be concentration gradients
on both large and small length scales in the environment. For example, the
infectious aerosol concentration at close range directly in front of an
infectious individual will usually be larger than the volume average of the
whole environment since exhaled puffs from the infectious individual will not
have dispersed much before inhalation. This means that a susceptible
individual located where they can inhale such puffs would be at greater risk
of infection than if they were not directly in the exhaled puffs of any
infectious individuals. The nature of the ventilation plays a significant role
in the validity of the assumption [30] The practice of social distancing,
using fans to better mix the room, etc. all improve the quality of this
assumption, but room conditions and people’s proximity to each other in real-
world situations can be far away from the well-mixed state with everyone
inhaling well-mixed air. We will make this same exact assumption in the model
presented in this manuscript, and will neither be using nor developing
corrections for close proximity between individuals and localized sinks and
other sources, though the nature of partial corrections will be briefly
discussed.

Besides the well-mixed assumption, there are several other assumptions
associated with Wells-Riley formulations, which are not necessarily always
true. As an example, there is an additional loss term that has not been fully
considered yet that is the loss of the infectious aerosols absorbed by the
individuals themselves, though the self-proximity depletion of infectious
aerosols in the vicinity of susceptible people has previously been mentioned
as an effect to consider [29]. This is despite the fact that this is exactly
the reason that susceptible individuals get infected. In some cases this can
be safely neglected, e.g. if the combined breathing volume exchange rate of
all individuals in the environment is negligible compared to that of
ventilation. But in a poorly ventilated room with many individuals inside,
this sink term must be taken into account – not incorporating it leads to
false risk predictions.

Another large assumption is that the absorbed doses follow a Poisson
distribution, which is implicit in the use of the exponential dose-response
model even if not stated explicitly [4, 6, 9, 1, 10], though there has been
work on doses following beta and gamma distributions [5]. The Poisson
distribution assumption requires that the pathogen-carrying aerosols have at
most one pathogen inside, i.e. a mono-pathogen assumption. However, this
assumption is violated if the pathogen concentration in an infectious
individual’s respiratory tract is high. For this poly-pathogen situation the
Wells-Riley formulation and the dose-response models must be generalized to
consider a larger number of pathogen in an individual aerosol explicitly. We
will use the term multiplicity to refer to the number of pathogen copies in an
aerosol.

Ignoring multiplicity causes the infection risk to be overestimated even
though the expected average pathogen dose does not change. Using a modified
version of the worked example later in this manuscript, Fig 1 shows this
effect on the time required to reach a 50% infection risk for different
pathogen concentrations in the respiratory tract fluid with and without
considering multiplicity. For low pathogen concentrations and small infection
probabilities per pathogen, ignoring multiplicity has only a small effect. But
for high pathogen concentrations and/or pathogen copies with a high infection
probability per pathogen, ignoring multiplicity has a significant impact. For
a respiratory tract pathogen concentration of $10^{11}$ ${\mathrm{cm}}^{-3}$
where the average number of pathogen copies per aerosol is approximately 6500
for a 50 $\mathrm{\SIUnitSymbolMicro m}$ in diameter at production, if the
single pathogen infection probability ($r$) is large enough that multiplicity
matters, this means taking into account multiplicities up to approximately
7000.

Fig 1: Effect of Ignoring Multiplicity. Ratio of the time required to reach a
50% infection risk when multiplicity is ignored $\tau_{50,ignore}$ to when it
is fully accounted for $\tau_{50,full}$ for single pathogen infection
probabilities $r$ (an average dose of $r^{-1}$ Poisson distributed pathogen
copies gives a mean infection risk of $63.21$%) and different pathogen
concentrations $\rho_{p}$ in the respiratory tract fluid of the infectious
individual as in the worked example later in the manuscript with a disease
following the exponential model, but at steady-state with just the speaking
mask-less infectious individual and the risk to a mask-less susceptible
individual whose exposure starts after steady state is reached. This is a
simplified version of Fig 5.

In this manuscript, we will consider the following generalizations and
modifications to the Wells-Riley formulation:

  * •

Fully accounting for the multiplicity of pathogen copies in aerosols and the
effect on the dose-response models.

  * •

Additional sink terms due to the filtering of air by people inhaling and then
exhaling it back out, including the effects of masks.

  * •

Working exclusively in units of pathogen copies and aerosols instead of quanta
(note, quantum is undefinable when accounting for multiplicity).

We will first generalize dose-response models that assume Poisson distributed
doses for the distribution that results from poly-pathogen aerosols being
present. Then we will develop the general pathogen concentration model that is
a generalization of the Wells-Riley formulation. This results in a linear
inhomogeneous coupled system of ODEs (Ordinary Differential Equations) for
each initial aerosol diameter at production (diameter when exhaled), with one
equation for each multiplicity that must be considered. We then derive the
general solution, and then simplify the general solution for coefficients that
are constant in time. Requirements and heuristics are developed for finding
the appropriate cutoff in the multiplicity, $M_{c}$. This is important because
the number of ODEs to solve is equal to $M_{c}$; and the computational effort
scales as $\mathcal{O}\left(M_{c}^{2}\right)$ for the numerical solution, or
worse than for $\mathcal{O}\left(M_{c}^{2}\right)$ or
$\mathcal{O}\left(M_{c}^{3}\right)$ for the different analytical solutions for
coefficients constant in time. Some circumstances allow small $M_{c}=1$ or
close to one. We consider a full hypothetical example situation for SARS-CoV-2
with very high viral loads to apply the generalized Wells-Riley formulation
developed in this manuscript. Finally, we discuss the effects of poly-pathogen
aerosols, the filtering by the people in the environment, the effects of face-
masks, and the model limitations and possible corrections. As a tool to aid
solving the model presented in this paper, we wrote the PMADRA (Poly-
Multiplicity Airborne Disease Risk Assessment) software suite
(https://gitlab.gwdg.de/mpids-lfpn-public/pmadra).

## Fundamentals

Throughout this manuscript, we will use the Poisson distribution, which
describes the probability of counting some number, $m$, of independent
events/objects/etc. as a function of the ensemble mean of the number counted,
$\mu$. The Probability Distribution Function (PDF) of the Poisson distribution
is

$P_{P}\left(\mu,m\right)=e^{-\mu}\left(\frac{\mu^{m}}{m!}\right)\quad.$ (1)

Most dose-response models assume that the number of pathogen copies absorbed
follows a Poisson distribution. For the case of a dose-response model, the
average number of pathogen copies absorbed over some period of time would be
the $\mu$ and then $P_{P}$ would give the probability that a person absorbed
exactly $m$ pathogen copies. For clarity in the rest of this manuscript, we
will now define $\Delta$ to be the number of pathogen copies absorbed (instead
of $m$) and the average number of pathogen copies absorbed is
$\left<\Delta\right>$, where we have used $\left<\cdot\right>$ to denote the
average. The use of a Poisson distribution for the doses requires that the
pathogen copies are independent (i.e. no clumping); and as we will later show,
that the number of pathogen copies in aerosols is assumed to be one or zero,
which is generally assumed by existing models but won’t be in the model
presented in this manuscript.

Let $R(\Delta)$ denote the infection probability when exactly $\Delta$
pathogen copies are absorbed, and $\mathfrak{R}(\left<\Delta\right>)$ denote
the average infection probability when the average number of pathogen copies
absorbed is $\left<\Delta\right>$. For a disease where the threshold (minimum
number of pathogen copies required for infection) is greater than one, the
threshold must be included into the definition of $R(\Delta)$ such that it is
zero for $\Delta$ less than the threshold, which makes $R(\Delta)$ be a piece-
wise function.

There are two ways to construct $\mathfrak{R}(\left<\Delta\right>)$ from
$R(\Delta)$. We use the method of taking the sum over all possible
$\Delta\in[1,\infty)$ of the product of the probability of absorbing each
particular $\Delta$ and the resulting infection risk $R(\Delta)$ [12]. If the
number of pathogen copies absorbed follows a Poisson distribution, then

$\mathfrak{R}\left(\left<\Delta\right>\right)=\sum_{\Delta=1}^{\infty}P_{P}\left(\left<\Delta\right>,\Delta\right)R(\Delta)\quad.$
(2)

The other method instead considers the number of pathogen copies that survive
to try to infect, $\Delta_{i}$, and does a double sum over $\Delta_{i}$
(starting from the threshold) and $\Delta$ of the product of the probability
of the dose $\Delta$ and the probability of exactly $\Delta_{i}$ out of
$\Delta$ surviving to try to infect [13] (this is NOT $R(\Delta)$). The two
methods are equivalent, with this extra sum being implicitly included in the
definition of $R(\Delta)$. This is why $R(\Delta)$ is a piece-wise function
when the threshold is not one. For some models it may be easier to do this
other method explicitly rather than try to construct $R(\Delta)$.

The exponential model assumes that all pathogen copies are identical, all
people are equally vulnerable to infection, the pathogen copies are acting
independently of each other, and that each pathogen has an equal probability
of causing infection $r$ [13]. These assumptions implicitly mean that the
threshold is one. Each pathogen has a probability $1-r$ to not infect. Then
the exponential model’s infection risk for an exact dose $\Delta$ is just one
minus the probability that all $\Delta$ pathogen copies did not infect.

$R_{E}(\Delta)=1-(1-r)^{\Delta}\quad.$ (3)

If the dose follows a Poisson distribution, then Eq (2) can be calculated for
the exponential model [12], yielding

$\mathfrak{R}_{E}\left(\left<\Delta\right>\right)=1-e^{-r\left<\Delta\right>}\quad.$
(4)

Note that often, the parameter $D\equiv 1/r$ is used instead of $r$ (the
symbol $k$ is also used [34]), which is the ID63.21 (Infective Dose required
for $63.21$% chance of infection). We will be making non-Poissonity
corrections to this later.

The beta-Poisson model is essentially the exponential model but instead of
considering everyone to be equally vulnerable, each person has their own value
for $r$ which comes from the beta distribution [12, 13]. The beta distribution
PDF [13] is

$P_{B}(r)=\frac{\Gamma(\varepsilon+\theta)}{\Gamma(\varepsilon)\Gamma(\theta)}r^{\varepsilon-1}(1-r)^{\theta-1}\quad,$
(5)

where $r\in[0,1]$ and the symbols $\varepsilon$ and $\theta$ have been used in
place of the conventional alpha and beta parameters respectively to avoid
clashing with symbols used later in this manuscript. This means that to get
the mean infection risk for a beta-Poisson model $R_{BP}(\Delta)$, we must
include an integral over all $r\in[0,1]$. Specifically,

$R_{BP}(\Delta)=\int_{0}^{1}{P_{B}(r)R_{E}(\Delta)dr}\quad.$ (6)

Since the integral commutes with the sums used to calculate
$\mathfrak{R}(\left<\Delta\right>)$, the integral can be calculated as an
outer integral rather than an inner integral yielding [13]

$\mathfrak{R}_{BP}(\left<\Delta\right>)=\int_{0}^{1}{P_{B}(r)\mathfrak{R}_{E}(\left<\Delta\right>)dr}\quad.$
(7)

Wells-Riley formulations, both the original model and many subsequent uses,
measure pathogen copies in units of quanta [4, 8, 9, 33, 12, 1, 10]. A quanta
is defined as ID63.21 pathogen copies [12]. This means that one quantum is
equal to $D=1/r$ pathogen copies. For the case of $r=1$ such as Mycobacterium
tuberculosis, one quantum is one pathogen [9, 12]. Using these units, the
exponential model from Eq (4) becomes

$\mathfrak{R}_{E}\left(\left<\mathbb{Q}\right>\right)=1-e^{-\left<\mathbb{Q}\right>}\quad,$
(8)

where $\mathbb{Q}$ is the number of absorbed quanta [4, 8, 9, 33, 12, 1, 10].

Let $N_{I}$ be the number of infectious individuals, $\sigma$ be the average
production rate of infectious quanta per infectious individual, $\lambda$ be
the volumetric breathing rate of susceptible individuals, $Q$ be the
volumetric rate that clean air is brought into the particular indoor
environment, and $\tau$ be the time period of exposure of susceptible
individuals. Then, in its simplest form, the Wells-Riley Model’s infection
probability for time periods smaller than the incubation period of the disease
[4] is

$\mathfrak{R}_{WR}\left(\tau\right)=1-\exp{\left[-\left(\frac{N_{I}\sigma}{Q}\right)\lambda\tau\right]}\quad.$
(9)

For time periods longer than the incubation period of the disease, one must
either break the time period into subintervals smaller than the incubation
period [4] or model both $\mathfrak{R}_{WR}$ and the number of infectious and
susceptible individuals over time with a SIR or SEIR model [6, 32].

## Dose-Response Models for poly-Pathogen Aerosols

### General

If the pathogen concentration in an infectious individual’s respiratory tract
fluid $\rho_{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}\leq 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^{11}$ ${\mathrm{cm}}^{-3}$ [35, 36] would give a
volume fraction of approximately $5\times 10^{-5}$, if we treat the virus as a
100 $\mathrm{nm}$ sphere (approximate size of the SARS-CoV-2 virus [37]). This
is important because an aerosol with a diameter of 1
$\mathrm{\SIUnitSymbolMicro m}$ could contain up to approximately 740
spherical pathogen copies with diameter 100 $\mathrm{nm}$, if we assume hard-
sphere packing (packing fraction of 74%). An aerosol with a diameter of 10
$\mathrm{\SIUnitSymbolMicro m}$ could contain up to approximately $7.4\times
10^{5}$ 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 $\Delta_{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 $\Delta_{k}=km_{k}$. The total pathogen dose from all aerosols is
just the sum of the doses for each multiplicity, which is
$\Delta=\sideset{}{{}_{k=1}^{\infty}}{\sum}\Delta_{k}$. Let
$\mu_{k}=\left<m_{k}\right>=\left<\Delta_{k}\right>/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 $\mu_{k}$. Since $\Delta_{k}=km_{k}$, the PDF of
$\Delta_{k}$ is not a Poisson distribution for $k>1$. It is instead a scaled-
Poisson distribution of the form

$P_{k}\left(\mu_{k},\Delta_{k}\right)=\begin{cases}P_{P}\left(\mu_{k},\frac{\Delta_{k}}{k}\right)&\text{if
$\quad\Delta_{k}\mod k=0$}\quad,\\\ 0&\text{otherwise}\quad.\end{cases}$

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
$\Delta_{k}$. Since $\Delta$ is the sum of a Poisson distribution for $k=1$
and some number of possibly non-negligible scaled-Poisson distributions, the
PDF of $\Delta$ 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 $\mu_{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\in[1,\infty)$ 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

$\mathcal{M}_{p}\left(\mu_{1},\dots,\mu_{\infty}\right)=\overbrace{\sum_{m_{1}=0}^{\infty}\dots\sum_{m_{\infty}=0}^{\infty}}^{\text{all
combinations}}\overbrace{\left[\prod_{k=1}^{\infty}{P_{P}\left(\mu_{k},m_{k}\right)}\right]}^{\text{probability
of
dose}}\bigg{[}\overbrace{R\Big{(}\underbrace{\sum_{k=1}^{\infty}km_{k}}_{\text{pathogen
dose}}\Big{)}}^{\text{infection probability}}\bigg{]}^{p}\quad,$ (10)

where we have written out the dose $\Delta$ inside $R$. The mean infection
risk is the first moment ($p=1$), which is

$\mathfrak{R}\left(\mu_{1},\dots,\mu_{\infty}\right)=\sum_{m_{1}=0}^{\infty}\dots\sum_{m_{\infty}=0}^{\infty}\left[\prod_{k=1}^{\infty}{P_{P}\left(\mu_{k},m_{k}\right)}\right]R\left(\sum_{k=1}^{\infty}km_{k}\right)\quad.$
(11)

### Exponential Model Corrections

Then, putting $R_{E}$ from Eq (3) into Eq (11), the exponential model mean
infection risk is

$\displaystyle\mathfrak{R}_{E}\left(\mu_{1},\dots,\mu_{\infty}\right)$
$\displaystyle=$
$\displaystyle\sum_{m_{1}=0}^{\infty}\dots\sum_{m_{\infty}=0}^{\infty}\left[\prod_{k=1}^{\infty}{P_{P}\left(\mu_{k},m_{k}\right)}\right]\left[1-\left(1-r\right)^{\sum_{k=1}^{\infty}km_{k}}\right]$
(12) $\displaystyle=$ $\displaystyle
1-\sum_{m_{1}=0}^{\infty}\dots\sum_{m_{\infty}=0}^{\infty}\prod_{k=1}^{\infty}{e^{-\mu_{k}}e^{(1-r)^{k}\mu_{k}}e^{-(1-r)^{k}\mu_{k}}\frac{\left[\left(1-r\right)^{k}\mu_{k}\right]^{m_{k}}}{m_{k}!}}$
$\displaystyle=$ $\displaystyle
1-\exp{\left[-\sum_{k=1}^{\infty}{\left(1-(1-r)^{k}\right)\mu_{k}}\right]}\quad,$

where the fact that the sum of all probabilities over the Poisson distribution
is equal to one has been used extensively. The final sum has a finite number
of terms due to the cutoff $M$ as long as the $\mu_{k}$ are finite for $k\leq
M$. For small $M_{c}$, we can truncate the risk probability and get an easier
to calculate approximation. Except for $M_{c}=1$, this is different from Eq
(4) due to the non-Poissonity in $\Delta$. The expression for the first few
values of $M_{c}$ are

$\displaystyle\mathfrak{R}_{E}\approx\begin{cases}1-e^{-r\mu_{1}}&\text{if
$M_{c}=1$}\quad,\\\ 1-e^{-r\mu_{1}}e^{-r\left(2-r\right)\mu_{2}}&\text{if
$M_{c}=2$}\quad,\\\
1-e^{-r\mu_{1}}e^{-r\left(2-r\right)\mu_{2}}e^{-r\left(3-3r+r^{2}\right)\mu_{3}}&\text{if
$M_{c}=3$}\quad.\end{cases}$ (13)

### 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

$\mathcal{M}_{BP,p}\left(\mu_{1},\dots,\mu_{\infty}\right)=\int_{0}^{1}{P_{B}(r)\mathcal{M}_{E,p}\left(\mu_{1},\dots,\mu_{\infty}\right)dr}\quad.$
(14)

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$.

## General Pathogen Concentration Model

### Looking Ahead

Now that we have dose-response models corrected for the multiplicity via Eq
(11), we must determine the average aerosol doses $\mu_{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_{0},t)$ be the
concentration density of aerosols with original diameter $d_{0}$ (diameter at
production) and $k$ pathogen copies in them over time, which has units of
${\mathrm{[L]}}^{-4}$ where $\mathrm{[L]}$ is the unit of length since
$n_{k}(d_{0},t)dd_{0}$ is the concentration of infectious aerosols with
original diameters between $d_{0}$ and $d_{0}+dd_{0}$. To get a concentration,
$n_{k}(d_{0},t)$ must be integrated with respect to $d_{0}$.

In the end, we will get the following system of ODEs (Ordinary Differential
Equations) in time $t$ and the original diameter at production $d_{0}$ for the
$n_{k}$, which is

$\frac{dn_{k}}{dt}=\overbrace{-\alpha(d_{0},t)n_{k}}^{\text{sinks}}+\overbrace{\left(k+1\right)\gamma(t)n_{k+1}-k\gamma(t)n_{k}}^{\text{flux
from inactivation}}+\overbrace{\beta_{k}(d_{0},t)}^{\text{sources}}\quad,$
(15)

where $\alpha(d_{0},t)$ is the sum of all sink term coefficients,
$\beta_{k}(d_{0},t)$ is the sum of all sources for each $k$, $\gamma(t)$ is
the pathogen inactivation rate, and we have assumed that the time period
considered is shorter than the incubation time of the disease. Then the
combined source and sink terms are

$\displaystyle\beta_{k}(d_{0},t)$ $\displaystyle=$
$\displaystyle\beta_{r,k}+\beta_{I,k}\quad,$ (16)
$\displaystyle\alpha(d_{0},t)$ $\displaystyle=$
$\displaystyle\alpha_{o}+\alpha_{r}+\alpha_{v}+\alpha_{g}+\alpha_{d}+\alpha_{I,f}+\alpha_{S,f}+\alpha_{O,f}\quad,$
(17)

which don’t depend on $n_{k}(d_{0},t)$ (i.e. no quadratic or higher order
terms), though they may depend on $t$. The different sources and sinks are
summarized in Table 1. See their relevant sections for the meanings of their
terms, their assumptions, and where they come from.

Table 1: Source And Sink Term Summary Summary of all the source (the $\beta$) and sink (the $\alpha$) 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. Term | Meaning | Form
---|---|---
$\beta_{r,k}\left(d_{0},t\right)$ | transport from other rooms | $q_{r}(t)n_{r,k}\left(d_{0},t\right)$
$\beta_{I,k}\left(d_{0},t\right)$ | production by infectious individuals | $\frac{N_{I}}{V}\left<\lambda_{I}(t)n_{I,k}(d_{0},t)\left[1-E_{I,m,out}(d_{0})\right]\right>_{I}$
$\alpha_{o}(t)$ | air exchange with outside | $q_{o}(t)$
$\alpha_{r}(t)$ | air exchange with other rooms | $q_{r}(t)$
$\alpha_{v}\left(d_{0},t\right)$ | filtering by ventilation | $q_{v}(t)E_{v}\left(w(d_{0},t)d_{0}\right)$
$\alpha_{g}\left(d_{0},t\right)$ | gravitational settling | $\approx\frac{1}{h}u_{g}\left(w(d_{0},t)d_{0}\right)$
$\alpha_{d}\left(d_{0},t\right)$ | deposition on surfaces | found elsewhere
$\alpha_{I,f}\left(d_{0},t\right)$ | filtering by infectious ind. inhaling | $\frac{1}{V}\sum_{j=1}^{N_{I}}\lambda_{I,j}(t)\left[1-S_{I,m,in,j}(d_{0},t)S_{I,r,j,k}(d_{0},w,\lambda_{I,j})S_{I,m,out,j,k}(d_{0})\right]$
$\alpha_{S,f}\left(d_{0},t\right)$ | filtering by susceptible ind. inhaling | $\frac{1}{V}\sum_{j=1}^{N_{S}}\lambda_{S,j}(t)\left[1-S_{S,m,in,j}(d_{0},t)S_{S,r,j,k}(d_{0},w,\lambda_{S,j})S_{S,m,out,j,k}(d_{0})\right]$
$\alpha_{O,f}\left(d_{0},t\right)$ | filtering by other ind. inhaling | $\frac{1}{V}\sum_{j=1}^{N_{O}}\lambda_{O,j}(t)\left[1-S_{O,m,in,j}(d_{0},t)S_{O,r,j,k}(d_{0},w,\lambda_{O,j})S_{O,m,out,j,k}(d_{0})\right]$

### 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}\equiv
Q_{o}/V$, $q_{r}\equiv Q_{r}/V$, and $q_{v}\equiv 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$
$\mathrm{\SIUnitSymbolMicro 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

d

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

de

equilibrium diameter in the environment

d0

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

dD

spherical 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_{0}$ 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 [26, 9,
25]. 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\approx 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\approx 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_{0}$. 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_{0}$ and thus $d\approx d_{0}$ 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_{0}$, which means that
some fraction of dry aerosols might not reach $d=d_{0}$ 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
$\alpha_{C,f}(d_{0},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 $\alpha_{C,f}(d_{0},t)$ sinks dominate, then a
better approximation or an explicit treatment of the diameter when exhaled
should be used. Combined, our assumption/approximation is

$d(t)\approx\begin{cases}d_{e}(t)&\text{if in the environment outside of the
respiratory tract}\quad,\\\ d_{0}&\text{at re-exalation after
inhalation}\quad.\end{cases}$ (18)

Let us define ratios between the remaining diameters: the evaporation ratio
$w$, the dilution ratio $\delta$, and the initial solute ratio $\zeta$ as

$\displaystyle w$ $\displaystyle\equiv$
$\displaystyle\frac{d_{e}}{d_{0}}\quad,$ (19) $\displaystyle\delta$
$\displaystyle\equiv$ $\displaystyle\frac{d_{e}}{d_{D}}\quad,$ (20)
$\displaystyle\zeta$ $\displaystyle\equiv$
$\displaystyle\frac{d_{D}}{d_{0}}\quad.$ (21)

Note that $w$ and $\delta$ are potentially functions of time, as well as
diameter due to the effect of surface curvature (through surface tension) on
equilibrium vapor pressure [27, 26]. 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 $\delta$ to be
different even for aerosols with the same $\zeta$.

But, we will make the assumption that the value of $\zeta$ 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 $\zeta$ 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 $\zeta$ and single
functions for $w$ and $\delta$, possibly over time and $d_{0}$ (or
equivalently $d_{D}$), for all infectious aerosols in the environment.

This means we can choose to track one of $d_{e}$, $d_{0}$, 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_{0}$, 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

$\frac{d\tilde{n}(d_{e},t)}{dt}=\frac{\partial\tilde{n}}{\partial
t}+\frac{\partial\tilde{n}}{\partial d_{e}}\,\frac{d_{e}}{dt}\quad.$ (22)

Since $d_{D}$ and $d_{0}$ 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_{0}$
or $d_{D}$, with $w$ appearing in places if $d_{0}$ is chosen, and both
$\delta$ and $\zeta$ appearing in places if $d_{D}$ is chosen. We choose
$d_{0}$ 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_{0}$ rather than $d_{D}$.

Now, $n_{k}(d_{0},t)$ is the concentration density of aerosols in terms of $t$
and the original diameter $d_{0}$. 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_{0}$ to $d_{0}+dd_{0}$, 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_{0}$, $\tilde{n}_{k}dd_{e}$, and $\breve{n}_{k}dd_{D}$. Thus, the
conversions are

$\displaystyle\tilde{n}_{k}$ $\displaystyle=$
$\displaystyle\frac{n_{k}}{w}\quad,$ (23) $\displaystyle\breve{n}_{k}$
$\displaystyle=$ $\displaystyle\frac{n_{k}}{\zeta}\quad,$ (24)
$\displaystyle\tilde{n}_{k}$ $\displaystyle=$
$\displaystyle\frac{\breve{n}_{k}}{\delta}\quad.$ (25)

Let $n_{0,k}(d_{0})$ be the initial concentration density in the room for a
multiplicity $k$ at the initial time $t=t_{0}$ and $n_{r,k}(d_{0},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_{0}\in[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_{0})$ be the largest number of
pathogen copies that can fit in an aerosol at production.. We will consider

$n_{k}(d_{0},t)=0\;\;\forall\;\;d\notin[d_{m,k},d_{M}],\;k>K_{m}(d_{0})\quad.$
(26)

All of these limits have problems, but there is no obvious better choice
without adding a lot more complexity to the model.

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_{0})$. Both of them are

$\displaystyle d_{m,k}$ $\displaystyle\approx$
$\displaystyle\begin{cases}d_{p}&\text{if $k=1$}\quad,\\\
\left(\frac{k}{e}\right)^{1/3}d_{p}&\text{if $k>1$}\quad,\end{cases}$ (27)
$\displaystyle K_{m}(d_{0})$ $\displaystyle\approx$
$\displaystyle\max{\left[1,e\left(\frac{d_{0}}{d_{p}}\right)^{3}\right]}\quad.$
(28)

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_{0}$ 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_{0}$, 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}(d_{M}^{3})$
for sufficiently large $d_{M}$ and pathogen concentration in the infectious
individual’s respiratory tract fluid $\rho_{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 $\rho_{w}$, the aerosols’ inertial response times $\tau_{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

$\displaystyle\tau_{p}$ $\displaystyle=$
$\displaystyle\frac{\rho_{w}d^{2}}{18\rho_{a}\nu_{a}}\quad,$ (29)
$\displaystyle u_{g}$ $\displaystyle=$
$\displaystyle\frac{\left(\rho_{w}-\rho_{a}\right)gd^{2}}{18\rho_{a}\nu_{a}}\approx
g\tau_{p}\quad,$ (30)

where $\rho_{a}$ is the density of air, $\nu_{a}$ is the kinematic viscosity
of air, and $g$ is the acceleration due to gravity.

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
$\mathrm{\SIUnitSymbolMicro m}$, 20 $\mathrm{\SIUnitSymbolMicro m}$, and 50
$\mathrm{\SIUnitSymbolMicro m}$ diameter aerosols; the terminal velocities at
20 $\mathrm{\SIUnitSymbolCelsius}$ and atmospheric pressure are $3.0$
$\mathrm{mm}\text{\,}{\mathrm{s}}^{-1}$, $1.2$
$\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$, and $7.5$
$\mathrm{cm}\text{\,}{\mathrm{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
$\mathrm{\SIUnitSymbolMicro m}$ aerosols are quite ballistic and quickly fall
out of the exhaled plume, but 10 $\mathrm{\SIUnitSymbolMicro 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 $\mathrm{\SIUnitSymbolMicro 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$ $\mathrm{\SIUnitSymbolMicro m}$, which will
be explored in the Discussion. Before evaporating, the terminal velocity is
$7.5$ $\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$. If the evaporation ratio is a
typical value in the $\frac{1}{2}\text{--}\frac{1}{5}$ range, the final
evaporated diameter would be in the 10–25 $\mathrm{\SIUnitSymbolMicro m}$
range and have terminal velocities in the 3–19
$\mathrm{mm}\text{\,}{\mathrm{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 $\left<\cdot\right>_{C}$ denote taking the average over all
individuals in category $C$.

Let $\lambda_{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_{0},w,\lambda_{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_{0}$ and $w$ give both the
initial diameter on inhalation ($d\approx d_{e}=wd_{0}$), the diameter the
aerosols grow towards ($d_{0}$) 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 $\lambda_{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_{0}$ 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

$\displaystyle S_{C,m,in,j}(d_{0},t)$ $\displaystyle=$ $\displaystyle
1-E_{C,m,in,j}\left(w(d_{0},t)d_{0}\right)\quad,$ (31) $\displaystyle
S_{C,r,j,k}(d_{0},w,\lambda_{C,j})$ $\displaystyle=$ $\displaystyle
1-E_{C,r,j}(d_{0},w,\lambda_{C,j})\quad,$ (32) $\displaystyle
S_{C,m,out,j,k}(d_{0})$ $\displaystyle=$ $\displaystyle
1-E_{C,m,out,j}(d_{0})\quad.$ (33)

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 $\rho_{j}\left(d_{0},t\right)dd_{0}$ be the number density in exhaled air
of the aerosols with diameters between $d_{0}$ and $d_{0}+dd_{0}$ exhaled by
the $j$’th infectious individual at time $t$. Let $\rho_{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_{0}$ is

$\left<k\right>\\!\left(d_{0},t\right)_{j}=\frac{\pi}{6}d_{0}^{3}\rho_{p,j}(t)\quad.$
(34)

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

$n_{I,j,k}(d_{0},t)=\begin{cases}\rho_{j}\left(d_{0},t\right)P_{P}\left(\left<k\right>\\!\left(d_{0},t\right)_{j},k\right)&\text{if
$d_{0}\geq d_{m,k}$}\quad,\\\ 0&\text{if $d_{0}<d_{m,k}$}\quad.\end{cases}$
(35)

Note that no infectious aerosols with multiplicity $k$ can be generated with
diameters too small to contain them (i.e. no $d_{0}<d_{m,k}$ aerosols).

### Sources

We will denote sources by the symbol $\beta$ 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

$\beta_{r,k}(d_{0},t)=q_{r}(t)n_{r,k}\left(d_{0},t\right)\quad.$ (36)

where we have lumped all other rooms that might be exchanging air with the
room of interest together rather than summing over them as done by Noakes &
Sleigh [33]. A coupled model for multiple rooms would have to split this into
a sum and model the whole system. Note that we are assuming, like elsewhere,
the aerosols brought in from other rooms reach their equilibrium diameter
quickly compared to other processes.

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

$\displaystyle\beta_{I,k}(d_{0},t)$ $\displaystyle=$
$\displaystyle\frac{1}{V}\sum_{j=1}^{N_{I}}{\overbrace{\lambda_{I,j}(t)n_{I,j,k}(d_{0},t)}^{\text{production
rate}}\overbrace{\left[1-E_{I,m,out,j}(d_{0})\right]}^{\text{mask survival}}}$
(37) $\displaystyle=$
$\displaystyle\frac{N_{I}}{V}\left<\lambda_{I}(t)n_{I,k}(d_{0},t)\left[1-E_{I,m,out}(d_{0})\right]\right>_{I}\quad,$

where the $j$ subscript has been dropped in the average. Any terms in the
average of a product ($\lambda_{I,j}$, $n_{I,j,k}$, and $1-E_{I,m,out,j,k}$)
that have no correlation with the others can be pulled out to make a product
of averages. But any correlated terms cannot be separated, which means it must
be kept as an average of a product. As an example, if there are two infectious
individuals in a room and one is singing and the other is listening in
silence; they will be strongly correlated. The singing person will on average
be breathing at a higher rate, could have a higher concentration density of
infectious aerosols in their exhaled air, and probably won’t be wearing a mask
while the listener might be wearing a mask. Now, if all individuals are
wearing the same mask, the mask term could be pulled out but the other two
terms would remain since they could still be correlated.

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
$\rho_{j}(d_{0},t)$ and $\left<k\right>\\!\left(d_{0},t\right)_{j}$ in
$n_{I,j,k}(d_{0},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 $\alpha$ 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 [33, 9] and are

$\displaystyle\alpha_{o}(t)$ $\displaystyle=$ $\displaystyle q_{o}(t)\quad,$
(38) $\displaystyle\alpha_{r}(t)$ $\displaystyle=$ $\displaystyle
q_{r}(t)\quad,$ (39)

respectively.

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\approx d_{e}=w(d_{0},t)d_{0}$. Then the volume normalized
loss rate coefficient from the circulating ventilation system [4] is

$\alpha_{v}(d_{0},t)=q_{v}(t)E_{v}\left(w(d_{0},t)d_{0}\right)\quad.$ (40)

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 $\alpha_{g}(w(d_{0},t)d_{0})$ and
$\alpha_{d}(w(d_{0},t)d_{0})$ 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

$\alpha_{g}\left(w(d_{0},t)d_{0}\right)\approx\frac{1}{h}u_{g}\left(w(d_{0},t)d_{0}\right)\quad,$
(41)

where $h$ is the characteristic height of the indoor environment and
$u_{g}(d)$ is the terminal velocity. For small spherical aerosols, Eq (30)
provides $u_{g}(d)$. Larger aerosols need additional diameter corrections [9,
25, 40].

### 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_{0},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
$\lambda_{C,j}(t)n_{k}(d_{0},t)$. The volume normalized sink coefficient from
this filtering is then

$\displaystyle\alpha_{C,f}(d_{0},t)$ $\displaystyle=$
$\displaystyle\frac{1}{V}\sum_{j=1}^{N_{C}}\overbrace{\lambda_{C,j}(t)}^{\text{volume
rate}}\overbrace{\Big{[}1-\overbrace{S_{C,m,in,j}(d_{0},t)}^{\text{mask
in}}\overbrace{S_{C,r,j,k}(d_{0},w,\lambda_{C,j})}^{\text{resp.
tract}}\overbrace{S_{C,m,out,j,k}(d_{0})}^{\text{mask
out}}\Big{]}}^{\text{total filtering efficiency}}$ (42) $\displaystyle=$
$\displaystyle\frac{N_{C}}{V}\left<\lambda_{C}(t)\left\\{1-\left[1-E_{C,m,in,j}\left(w(d_{0},t)d_{0}\right)\right]\right.\right.$
$\displaystyle\bullet\left.\left.\left[1-E_{C,r}(d_{0},w(d_{0},t),\lambda_{C}(t))\right]\left[1-E_{C,m,out}(d_{0})\right]\right\\}\right>_{C}\quad,$

where the $j$ subscript has been dropped in the average over category $C$. As
was the case before with the average of a product, only terms that are
uncorrelated with the others can be pulled out or be replaced by their average
value inside. Note that if aerosols completely dry out in the environment, we
have made the assumption that their diameters have approximately returned to
$d_{0}$ upon leaving the respiratory tract at re-exhalation. This assumption
only effects the value of $\alpha_{C,f}(d_{0},t)$ if an individual is wearing
a mask.

### 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 $\gamma(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

$f_{k,k-1}(t)n_{k}(d_{0},t)=k\gamma(t)n_{k}(d_{0},t)\quad.$ (43)

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

$\frac{dn_{k}}{dt}=-\alpha(d_{0},t)n_{k}+f_{k+1,k}(d_{0},t)n_{k+1}-f_{k,k-1}(d_{0},t)n_{k}+\beta_{k}(d_{0},t)\quad.$
(44)

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$ $\mathrm{\SIUnitSymbolMicro 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

$\frac{dn_{k}}{dt}=-\alpha(d_{0},t)n_{k}+\left(k+1\right)\gamma(t)n_{k+1}-k\gamma(t)n_{k}+\beta_{k}(d_{0},t)\quad.$
(45)

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_{0}$ 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_{0}$.
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 $\alpha$ and $\gamma$ but different
sources $\beta_{k,1}$ and $\beta_{k,2}$ respectively, then the solution for
$\beta_{k}=\beta_{k,1}+\beta_{k,2}$ is $n_{k}=n_{k,1}+n_{k,2}$.

### Infection Risk

Let $\mu_{j,k}$ be the average number of aerosols with multiplicity $k$
absorbed by the $j$’th susceptible individual from time $t_{0}$ to time $t$.
At any particular instant of time, the average number of such aerosols of each
original diameter $d_{0}$ entering the person’s mask if they are wearing a
mask or their mouth and nose if they aren’t is
$\lambda_{S,j}(t)n_{k}(d_{0},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_{0},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_{0},t)$ will survive the mask to enter the
respiratory tract [5, 6, 7, 8, 12, 10]. A fraction $E_{S,r,j}(d_{0})$ 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_{0}$ and the time between $t_{0}$ and $t$, which
is

$\displaystyle\mu_{j,k}(t)$ $\displaystyle=$
$\displaystyle\int_{d_{m,k}}^{d_{M}}{d\phi\,\int_{t_{0}}^{t}{dv\,\overbrace{E_{S,r,j}(\phi,w(\phi,v),\lambda_{S,j}(v))}^{\text{absorption
efficiency}}\overbrace{S_{S,m,in,j}(\phi,t)}^{\text{survive
mask}}\overbrace{\lambda_{S,j}(v)n_{k}(\phi,v)}^{\text{inf. aerosol inhalation
rate}}}}$ (46) $\displaystyle=$
$\displaystyle\int_{d_{m,k}}^{d_{M}}\\!\\!\\!\\!d\phi\\!\int_{t_{0}}^{t}dv\,\Biggl{\\{}E_{S,r,j}(\phi,w(\phi,v),\lambda_{S,j}(v))$
$\displaystyle\quad\bullet\left[1-E_{S,m,in,j}\left(w(\phi,v)\phi\right)\right]\lambda_{S,j}(v)n_{k}(\phi,v)\Biggr{\\}},$

where we have $\phi$ as the integration variable over $d_{0}$ and $v$ as the
integration variable over time. We will continue to use $\phi$ and $v$
exclusively for this purpose in the rest of the manuscript.

In order to use the $\mu_{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 $\mu_{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=\tau_{p}/\tau_{\eta}$ are very small ($St\ll 1$) [41, 42]
where $\tau_{p}$ is the aerosol inertial response time scale from Eq (29) and
$\tau_{\eta}$ is the Kolmogorov time scale of the turbulence in the
environment, which is $\tau_{\eta}=\sqrt{\nu_{a}/\epsilon}$ where $\epsilon$
is the turbulent dissipation rate. It will also be small if the typical inter-
aerosol distance $\bar{d}_{a}\sim\mathcal{N}^{-1/3}$, where $\mathcal{N}$ is
the total infectious aerosol concentration for all $d_{0}$ 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\eta$ [41, 42]
where $\eta=(\nu_{a}^{3}/\epsilon)^{1/4}$ is the Kolmogorov length scale of
the turbulence. This means that as long as $St\ll 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)\approx 1$. For a low dissipation rate of $\epsilon=1$
$\mathrm{mW}\text{\,}{\mathrm{kg}}^{-1}$; $St=0.06$ for a $d_{M}$ sized
aerosol and the number density limit is $\mathcal{N}\ll 4\times 10^{5}$
${\mathrm{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 $\epsilon=1$
$\mathrm{W}\text{\,}{\mathrm{kg}}^{-1}$; $St=2.0$ for a $d_{M}$ sized aerosol
and the number density limit is $\mathcal{N}\ll 7\times 10^{7}$
${\mathrm{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 $\epsilon=10$
$\mathrm{W}\text{\,}{\mathrm{kg}}^{-1}$; $St=6.3$ for a $d_{M}$ sized aerosol
and the number density limit is $\mathcal{N}\ll 4\times 10^{8}$
${\mathrm{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.

## Model Solution and Simplification

### General

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

$\frac{d\vec{n}}{dt}=\mathbf{A}(d_{0},t)\vec{n}(d_{0},t)+\vec{\beta}(d_{0},t)\quad,$
(47)

where $\vec{n}(d_{0},t)$ and $\vec{\beta}(t)$ are the $n_{k}(d_{0}t)$ and
$\beta_{k}(d_{0},t)$ for $k>0$ in vector form and

$\mathbf{A}\equiv\begin{bmatrix}-\alpha(d_{0},t)-\gamma(t)&2\gamma(t)\\\
&-\alpha(d_{0},t)-2\gamma(t)&3\gamma(t)\\\ &&\ddots&\ddots\\\
&&&\ddots&M_{c}\gamma(t)\\\
&&&&-\alpha(d_{0},t)-M_{c}\gamma(t)\end{bmatrix}\quad.$ (48)

is an upper bidiagonal $M_{c}\times M_{c}$ square matrix. For any fixed
$d_{0}$ or bin of $d_{0}$, the resulting system of ODEs is particularly
amenable to efficient numerical solution even for very large $M_{c}$ because
$\mathrm{A}$ is sparse with only one or two elements per row.

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

$\vec{n}(d_{0},t)=\exp{\left[\int_{t_{0}}^{t}{A(d_{0},x)dx}\right]}\vec{n}_{0}(d_{0})+\int_{t_{0}}^{t}{\exp{\left[\int_{s}^{t}{A(d_{0},x)dx}\right]}\vec{\beta}(d_{0},s)ds}\quad.$
(49)

Working this out using the structure of the diagonalization of $\mathbf{A}$ in
S1 Appendix., the general solution for each $k$ is

$n_{k}(d_{0},t)=\exp{\left[-\int_{t_{0}}^{t}{\alpha(d_{0},x)dx}\right]}\exp{\left[-k\int_{t_{0}}^{t}{\gamma(x)dx}\right]}\\\
\bullet\sum_{p=k}^{M_{c}}{\binom{p}{k}n_{0,p}(d_{0})\left[1-\exp{\left[-\int_{t_{0}}^{t}{\gamma(x)dx}\right]}\right]^{p-k}}\\\
+\sum_{p=k}^{M_{c}}\binom{p}{k}\int_{t_{0}}^{t}\beta_{p}(d_{0},s)\\\
\bullet\exp{\left[-\int_{s}^{t}{\alpha(d_{0},x)dx}\right]}\exp{\left[-k\int_{s}^{t}{\gamma(x)dx}\right]}\left[1-\exp{\left[-\int_{s}^{t}{\gamma(x)dx}\right]}\right]^{p-k}ds\quad,$
(50)

where $\binom{k}{m}=k!/(m!(k-m)!)$ is the notation for the binomial
coefficient $k$ choose $m$.

### Coefficients Constant in Time

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

$\displaystyle\vec{n}_{\infty}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha}\vec{\beta}\quad,$ (51) $\displaystyle\vec{n}$
$\displaystyle=$
$\displaystyle\vec{n}_{\infty}+\left(\vec{n}_{0}-\vec{n}_{\infty}\right)e^{(t-t_{0})\alpha}\quad,$
(52) $\displaystyle\int_{t_{0}}^{t}{\vec{n}(v)dv}$ $\displaystyle=$
$\displaystyle(t-t_{0})\vec{n}_{\infty}+\frac{1}{\alpha}\left(\vec{n}_{0}-\vec{n}_{\infty}\right)\left[1-e^{(t-t_{0})\alpha}\right]\quad.$
(53)

For the trivial case that both $\gamma=0$ and $\alpha=0$, the solution is
instead

$\displaystyle n_{\infty,k}$ $\displaystyle=$
$\displaystyle\begin{cases}0&\text{if $\beta_{k}=0$}\quad,\\\
+\infty&\text{otherwise}\quad,\end{cases}$ (54) $\displaystyle\vec{n}$
$\displaystyle=$ $\displaystyle\vec{n}_{0}+(t-t_{0})\vec{\beta}\quad,$ (55)
$\displaystyle\int_{t_{0}}^{t}{\vec{n}(v)dv}$ $\displaystyle=$
$\displaystyle(t-t_{0})\vec{n}_{0}+\frac{1}{2}(t-t_{0})^{2}\vec{\beta}\quad.$
(56)

But for the general case of $\gamma\neq 0$, the solution is instead (see S1
Appendix.)

$\displaystyle n_{k}(d_{0},t)$ $\displaystyle=$ $\displaystyle
n_{\infty,k}+z^{s}\left[U_{k}(d_{0},\vec{\beta}(d_{0}),z)+V_{k}(\vec{n}_{0}(d_{0}),z)\right]\quad,$
(57) $\displaystyle\int_{t_{0}}^{t}{n_{k}(d_{0},v)dv}$ $\displaystyle=$
$\displaystyle(t-t_{0})n_{\infty,k}(d_{0})$ (58) $\displaystyle\quad-
U_{k}(d_{0},\vec{n}_{0}(d_{0}),1)+z^{s}U_{k}(d_{0},\vec{n}_{0}(d_{0}),z)$
$\displaystyle\quad-\frac{1}{\gamma}W_{k}\left(d_{0},\vec{\beta},z\right)\quad,$

where

$\displaystyle z(t)$ $\displaystyle=$ $\displaystyle
e^{-(t-t_{0})\gamma}\in(0,1]\quad,$ (59) $\displaystyle s(d_{0})$
$\displaystyle=$ $\displaystyle\frac{\alpha(d_{0})}{\gamma}+k\quad,$ (60)
$\displaystyle V_{k}(\vec{y},x)$ $\displaystyle=$
$\displaystyle\sum_{i=k}^{M_{c}}\binom{i}{k}y_{i}(1-x)^{i-k}\quad,$ (61)
$\displaystyle U_{k}(d_{0},\vec{y},x)$ $\displaystyle=$
$\displaystyle-\frac{1}{\gamma}\sum_{i=k}^{M_{c}}\binom{i}{k}y_{i}\sum_{p=0}^{i-k}\binom{i-k}{p}\frac{(-1)^{p}x^{p}}{s+p}\quad,$
(62) $\displaystyle W_{k}(d_{0},\vec{y},x)$ $\displaystyle=$
$\displaystyle\int_{1}^{x}{dv\,v^{s-1}U_{k}(d_{0},\vec{y},v)}\quad,$ (63)
$\displaystyle=$
$\displaystyle-\frac{1}{\gamma}\sum_{i=k}^{M_{c}}\binom{i}{k}\beta_{i}(d_{0})\sum_{p=0}^{i-k}\binom{i-k}{p}\frac{(-1)^{p}\left(z^{s+p}-1\right)}{(s+p)^{2}}\quad,$
(64)

and $n_{\infty,k}(d_{0})$ is the concentration density as $t\rightarrow\infty$
which is

$n_{\infty,k}(d_{0})=-U_{k}\left(d_{0},\vec{\beta},1\right)=\frac{1}{\gamma}\sum_{i=k}^{M_{c}}{\binom{i}{k}\beta_{i}(d_{0})\sum_{p=0}^{i-k}{\binom{i-k}{p}\frac{(-1)^{p}}{s+p}}}\quad,$
(65)

Note that $s$ is a function of $k$ and $e^{-(\alpha+k\gamma)(t-t_{0})}=z^{s}$.

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

$\mu_{j,k}(t)=\lambda_{S,j}\int_{d_{m,k}}^{d_{M}}d\phi
E_{S,r,j}(\phi,w,\lambda_{S,j})\left(1-E_{S,m,in,j}\left(w\phi\right)\right)\int_{t_{0}}^{t}{n_{k}(\phi,v)dv}\quad.$
(66)

Calculation of $\vec{n}_{k}(d_{0},t)$, $\vec{n}_{\infty,k}(d_{0})$,
$\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 $\vec{n}_{\infty,k}(d_{0})$ 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

$\displaystyle U_{k}(d_{0},\vec{y},x)$ $\displaystyle=$
$\displaystyle\begin{cases}-\frac{y_{M_{c}}}{\gamma s}&\text{if
$k=M_{c}$}\quad,\\\ \frac{(k+1)x}{s}U_{k+1}(d_{0},\vec{y},x)-\frac{1}{\gamma
s}V_{k}(\vec{y},x)&\text{otherwise}\quad,\end{cases}$ (67) $\displaystyle
W_{k}(d_{0},\vec{y},x)$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{y_{M_{c}}}{\gamma
s^{2}}\left(1-x^{s}\right)&\text{if $k=M_{c}$}\quad,\\\
\frac{1}{s}\big{[}(k+1)W_{k+1}\left(d_{0},\vec{y},x\right)&\\\
\quad+x^{s}U_{k}\left(d_{0},\vec{y},x\right)-U_{k}\left(d_{0},\vec{y},1\right)\big{]}&\text{otherwise}\quad,\end{cases}$
(68) $\displaystyle U_{k}\left(d_{0},\vec{y},1\right)$ $\displaystyle=$
$\displaystyle\begin{cases}-\frac{y_{M_{c}}}{\gamma s}&\text{if
$k=M_{c}$}\quad,\\\
\frac{(k+1)x}{s}U_{k+1}(d_{0},\vec{y},1)-\frac{y_{k}}{\gamma
s}&\text{otherwise}\quad,\end{cases}$ (69) $\displaystyle n_{\infty,k}$
$\displaystyle=$ $\displaystyle\begin{cases}\frac{\beta_{M_{c}}}{\gamma
s}&\text{if $k=M_{c}$}\quad,\\\ \frac{1}{\gamma
s}\left[\beta_{k}+(k+1)\gamma\,n_{\infty,k+1}\right]&\text{otherwise}\quad.\end{cases}$
(70)

This recursive analytical solution for $\vec{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^{-12}$. 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.

## Determining The Cutoff Mc

In order to reduce the number of equations that have to be solved, we need to
find a suitable cutoff $M_{c}<M$ if at all possible, whether for the whole
diameter range or for each diameter bin (advantage of doing a separate one for
each bin is that $M_{c}$ tends to be small for the small diameter bins), such
that the contribution of all higher multiplicities is less than a threshold
$T\in(0,1]$ fraction of the total contribution from all multiplicities. In
many cases, this depends only on the $\rho_{p,j}$ of the infectious
individuals and one can skip directly to Eq (80) for the value of $M_{c}$ to
use (shown in Fig 2 for a few $\rho_{p,j}$). However, some cases such as when
one starts the model after some number of infectious individuals have left the
environment, when there is significant transport from other rooms, etc.
require additional heuristics. These heuristics are developed below.

A cutoff is suitable if the total contribution for all $k>M_{c}$ to the
average pathogen dose and therefore infection risk is small compared to the
total contribution for $k\leq M_{c}$. It is almost always true that $M_{c}<M$,
and in many cases it can even be $M_{c}=1$. This depends on the distribution
of exhaled aerosol sizes and the pathogen concentration $\rho_{p}$ in the
respiratory tract fluid where the aerosols are produced. For very low pathogen
loading, one can use $M_{c}=1$. Let $d_{-}$ and $d_{+}$ be the bounds in
$d_{0}$ of the bin (or whole range in which case $d_{-}=d_{m,1}$ and
$d_{+}=d_{M}$) being considered.

The most reliable way to determine $M_{c}$ is to use the model with the cutoff
$M$ and determine $M_{c}$ afterwards using the result, but that defeats the
point of finding $M_{c}$ since the effort one wants to save has already been
expended. So we need heuristics to determine $M_{c}$ ahead of time. All of
them consider the dose contribution from high multiplicity aerosols and
consider a simplified $k\mu_{j,k}$ from Eq (46) with a particular
concentration density multiplied by the average absorption efficiency of
susceptible individuals. For each heuristic, we will define this parameter to
be $\mathcal{H}_{h,k}(t)$ where the $h$ denotes the particular heuristic.
Then, the heuristic for $M_{c}$ is that we must find the $M_{c}$ such that

$\sum_{k=1}^{M_{c}}{\mathcal{H}_{h,k}(t)}\gg\sum_{k=M_{c}+1}^{\infty}{\mathcal{H}_{h,k}(t)}\;\;\forall\;\;h,t\geq
t_{0}\quad.$ (72)

Note that we must take the largest $M_{c}$ out of the values suggested by the
individual heuristics.

An equivalent way to express this heuristic is to look at the ratio of the sum
of $\mathcal{H}_{h,k}$ after the cutoff ($k>M_{c})$ to the total, defined as

$J_{h,M_{c}}(t)\equiv\frac{\sum_{k=M_{c}+1}^{\infty}{\mathcal{H}_{h,k}(t)}}{\sum_{k=1}^{\infty}{\mathcal{H}_{h,k}(t)}}\quad.$
(73)

Now, $J_{h,M_{c}}(t)\in[0,1]$ and is approximately the ratio of the
contribution of the higher multiplicities $k>M_{c}$ aerosols to the total,
which we want to be small. An equivalent statement of the heuristics is that
one must find the $M_{c}$ such that $J_{h,M_{c}}\ll 1\;\forall\;h,t\geq
t_{0}$. One way to determine $M_{c}$ is to say pick some threshold
$T\in(0,1]$, and then find the smallest $M_{c}$ such that $J_{h,M_{c}}\leq T$
for all heuristics. Let $M_{c,h}(T)$ be the smallest value of $M_{c}$ that
satisfies $J_{h,M_{c,h}}(t)\leq T$, which makes it the single heuristic value
of $M_{c}$. Then, $M_{c}$ is just the maximum $M_{c,h}$.

First, we define the average absorption efficiency of the susceptible
individuals as

$A_{S}(d_{0},t)\equiv\left<E_{S,r}(d_{0},w(d_{0},t),\lambda_{S})\left[1-E_{S,m,in}\left(w(d_{0},t)d_{0}\right)\right]\right>_{S}\quad.$
(74)

If the $\alpha$, $\beta$, $\gamma$, and $w$ are constant in time; it is a lot
less effort to calculate $n_{\infty,k}(d_{0})$ using Eq (65) than
$n_{k}(d_{0})$. Then, each $\mu_{j,k}\sim A_{S}n_{\infty,k}$. If $q_{r}(t)$
and $n_{r,k}(d_{0},t)$ are non-zero, the doses from them have a similar
scaling. If the initial concentration density includes a lot of aerosols with
high multiplicities, we will need to set $M_{c}$ to be large enough to include
them even if they won’t matter after the initial time. We need to consider
this if $n_{0,k}\gg n_{\infty,k}$ for any $k>1$, and they will have a similar
scaling. These heuristics are

$\displaystyle\mathcal{H}_{\infty,k}$ $\displaystyle=$ $\displaystyle
k\int_{d_{-}}^{d_{+}}{A_{S}(\phi,t)n_{\infty,k}(\phi)d\phi}\quad,$ (75)
$\displaystyle\mathcal{H}_{r,k}(t)$ $\displaystyle=$ $\displaystyle
k\int_{d_{-}}^{d_{+}}{A_{S}(\phi,t)n_{r,k}(\phi,t)d\phi}\quad,$ (76)
$\displaystyle\mathcal{H}_{0,k}(t)$ $\displaystyle=$ $\displaystyle
k\int_{d_{-}}^{d_{+}}{A_{S}(\phi,t)n_{0,k}(\phi)d\phi}\quad.$ (77)

The last heuristic is similar but considers the infectious individuals inside
the environment instead of the concentration density. This has the advantage
of not needing to determine $n_{\infty,k}(d_{0})$. We essentially take the
average over the $d_{0}$ interval of $\beta_{I,k}(d_{0})$ from Eq (36) times
the absorption efficiency of the average susceptible individual. We thus
define the infectious individuals heuristic parameter

$\mathcal{H}_{I,k}(t)\equiv
k\int_{d_{-}}^{d_{+}}d\phi\,A_{S}(\phi,t)\sum_{j=1}^{N_{I}}{\lambda_{i,j}(t)n_{I,j,k}(\phi,t)\left[1-E_{I,m,out,j}(\phi)\right]}\quad.$
(78)

But there are practical difficulties in using it directly. So instead, we will
define the heuristic for each individual infectious individual using the
largest diameter in the range $d_{+}$, and one would use the maximum $M_{c}$
indicated by all of these. This has the advantage that there is a simple form
for the required $M_{c}$, which is derived in S4 Appendix.. It is

$M_{c,I,j}(d_{+},T)=1+C_{P}^{-1}\left(\left<k\right>\\!\left(d_{+},t\right)_{j},\left(1-T\right)C_{P}\left(\left<k\right>\\!\left(d_{+},t\right)_{j},K_{m}(d_{+})-1\right)\right)\quad,$
(79)

where $C_{P}$ is the CDF (Cumulative Distribution Function) of the Poisson
distribution and $C_{P}^{-1}(\mu,c)$ is the inverse CDF to find the smallest
$k$ for which $C_{P}(\mu,k)\geq c$. Note that when $K_{m}(d_{+})\gg 1$ and
$K_{m}(d_{+})\gg\left<k\right>\\!\left(d_{+},t\right)_{j}\gg 1$,
$C_{P}\left(\left<k\right>\\!\left(d_{+},t\right)_{j},K_{m}-1\right)\simeq 1$
and

$M_{c,I,j}(d_{+},T)\simeq
1+C_{P}^{-1}\left(\left<k\right>\\!\left(d_{+},t\right)_{j},\left(1-T\right)\right)\quad.$
(80)

When the assumptions don’t apply, this will give an overestimation, so it is
usable to get the value of $M_{c}$ to use. It will just give a bigger value
than necessary.

Fig 2 shows $M_{c,I,j}$ as a function of $d_{0}$ for several different
$\rho_{p,j}$. Increasing $\rho_{p,j}$ approximately just shifts the curves for
$M_{c,I,j}$ to the left on a log-scale. Notice the very strong effect of
$\rho_{p,j}$ on $M_{c}$, with values a little under 7000 being required for
the largest diameter bin for $\rho_{p,j}=10^{11}$ ${\mathrm{cm}}^{-3}$ and a
value of 2 being required for the same bin for $\rho_{p,j}=10^{6}$
${\mathrm{cm}}^{-3}$. Since $M_{c}$ increases with $d_{0}$, the vast majority
of the effort to determine the concentration density and the infection risk
will be spent on the largest bins except for small values of $\rho_{p,j}$.

Fig 2: Required Mc Based on Pathogen Concentration in Infectious Individuals.
$M_{c,I,j}$ required to capture 99% of pathogen production for each diameter
at aerosol production $d_{0}$ from an infectious individual, with each line
being a different pathogen concentration in their respiratory tract fluid
$\rho_{p,j}$ (see legend).

## Example for SARS-CoV-2 with High Viral Load

### 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_{0},t)$ need to be applied at any source or sink, nor
in the calculated absorbed doses. Let the room have volume $V=200$
${\mathrm{m}}^{3}$ with a height of $h=4$ $\mathrm{m}$, with ventilation
$q_{r}=0$, $q_{v}=0$, and $q_{o}=0.5$ ${\mathrm{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_{0}$. We ignore deposition ($\alpha_{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
$\lambda_{S,j}=0.3$ ${\mathrm{m}}^{3}\text{\,}{\mathrm{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}$–$10^{11}$ ${\mathrm{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})=0$. Stage 1; at $t=t_{0}=0$, one
infectious individual enters the room who is speaking, wearing no mask,
breathing at a rate $\lambda_{I,j}=0.5$
${\mathrm{m}}^{3}\text{\,}{\mathrm{hr}}^{-1}$ (just below an 0.54
${\mathrm{m}}^{3}\text{\,}{\mathrm{hr}}^{-1}$ average value for reading out
loud [15]), and has a high respiratory tract fluid pathogen concentration of
$\rho_{p,j}=10^{10}$ ${\mathrm{cm}}^{-3}$. Stage 2; then at $t=3$
$\mathrm{hr}$, one more infectious individual enters the room who is
continuously coughing while wearing a simple2 mask, breathing at a higher rate
of $\lambda_{I,j}=2.0$ ${\mathrm{m}}^{3}\text{\,}{\mathrm{hr}}^{-1}$, and has
a higher respiratory tract fluid pathogen concentration of
$\rho_{p,j}=10^{11}$ ${\mathrm{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 ${\mathrm{m}}^{3}\text{\,}{\mathrm{hr}}^{-1}$ from which we chose 2.0
${\mathrm{m}}^{3}\text{\,}{\mathrm{hr}}^{-1}$.

We use mask filter efficiencies of the functional form

$E_{C,m,in,j}(d)=E_{C,m,out,j}(d)=E_{\infty}-\left(E_{\infty}-E_{0}\right)e^{-d/D_{m,c}}\quad,$
(81)

where $E_{\infty}$ is the aerosol filtering efficiency as
$d\rightarrow\infty$, $E_{0}$ is the aerosol filtering efficiency as
$d\rightarrow 0$, and $D_{m,c}$ is the scale of the mask efficiency
transition. We will use $D_{m,c}=10$ $\mathrm{\SIUnitSymbolMicro m}$. We
consider individuals wearing no masks or one of two types of masks. Their
filtering efficiencies are

none (no mask)

$E_{0}=E_{\infty}=0$.

mask simple1

$E_{0}=0.2$ and $E_{\infty}=0.8$.

mask simple2

$E_{0}=0.95$ and $E_{\infty}=0.99$.

The filtering efficiencies of both the simple1 and simple2 masks are shown in
S6 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}(t))=\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\times 10^{-3}$ and the same
value of $r$ as found for HCoV-229E in humans which is $r=5.39\times 10^{-2}$
[34]. We use $\gamma=0.64$ ${\mathrm{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
$\mathrm{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 $\mathrm{nm}$ to 100
$\mathrm{nm}$). This is used with Eq (34) and (35) to get the $\beta_{I,k}$.
They are shown in the top-right panel of Fig 3. The aerosol size distributions
have two peaks at approximately 2 $\mathrm{\SIUnitSymbolMicro m}$ and 100
$\mathrm{\SIUnitSymbolMicro 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_{0}$ at $t=6$ $\mathrm{hr}$ compared to the aerosol concentration densities
being exhaled by speaking and coughing individuals from Johnson et al. [22]
scaled by $10^{-4}$ 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\times 10^{-3}$ (Bottom-Right) $r=5.39\times
10^{-2}$.

### 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$ $\mathrm{\SIUnitSymbolMicro m}$ and $d_{M}=50$
$\mathrm{\SIUnitSymbolMicro 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, $\left.\beta_{I,k}\right|_{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 $\alpha$ and $\beta_{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^{-3}$. 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 $\left.n_{k}\right|_{i}(t)$ and
$\left.\mu_{j,k}\right|_{i}(t)$ are solved analytically if $M_{c}\leq 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^{-4}$ $\mathrm{hr}$, which is required for stability and an accurate
solution with the large $\left.\alpha\right|_{i}+M_{c}\gamma$ values in the
largest bin. After determining $\left.\alpha\right|_{i}$ and
$\left.\vec{\beta}\right|_{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_{0}$ (see top-right panel of Fig 3); we can see how as $d_{0}$
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_{0}>15$ $\mathrm{\SIUnitSymbolMicro m}$ the increasing $\alpha$ 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
$\alpha$ and total sinks $\alpha+k\gamma$ are shown in Fig 4 and we can see
that settling causes $\alpha$ to increase by over a factor of 10 from 100
$\mathrm{nm}$ to 50 $\mathrm{\SIUnitSymbolMicro m}$. 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 $\alpha$ without inactivation, $\alpha+\gamma$ for $k=1$,
and $\alpha+M_{c}\gamma$ for $k=M_{c}$ (different values for Stage 1 and 2).

The pathogen concentrations as a function of $d_{0}$ and $k$ right after the
beginning and at the end of each Stage are shown in S7 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
S8 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 $\alpha(d_{0})$ and $\beta_{k}(d_{0})$ 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\times 10^{-3}$, but more
noticeable but still small for the larger $r=5.39\times 10^{-2}$.

## Discussion

### 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
$\beta_{k,new}=0\;\forall\;k\neq 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_{0}$ 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 $\rho_{p}$ of the first
speaking infectious individual was adjusted, and we took the steady state case
where $\vec{n}_{0}=\vec{n}_{\infty}$ and calculated the constant
$d\mu_{j,k}/dt$ for each susceptible individual. Then using the constant
$d\mu_{j,k}/dt$ and an initial dose of zero, we found the time, $\tau_{50}$,
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 $\tau_{50}$ to be underestimated (overestimation of risk).
The underestimate of $\tau_{50}$ is shown in Fig 5.

Fig 5: Effect of Ignoring Multiplicity, Full Version. Full version of Fig 1
with more $\rho_{p}$ and the effect of masks. Plot of the ratio of the time
required to reach a 50% infection risk when multiplicity is ignored
$\tau_{50,ignore}$ to when it is fully accounted for $\tau_{50,full}$ for
different respiratory tract fluid pathogen concentrations $\rho_{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 $\mathrm{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 $\rho_{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 $\rho_{p}$
increases, the expected multiplicity range for each $d_{0}$ increases which
makes ignoring multiplicity underestimate $\tau_{50}$ more. For the $r$ values
considered here, $\rho_{p}\leq 10^{9}$ ${\mathrm{cm}}^{-3}$ underestimates
$\tau_{50}$ by at most 20% and $\rho_{p}\leq 10^{8}$ ${\mathrm{cm}}^{-3}$
underestimates it by at most 12%. But for $\rho_{p}=10^{11}$
${\mathrm{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 $\left<\Delta\right>=r^{-1}$ and all
infectious aerosols have the exact same multiplicity $k$. Then, the $\mu$ for
all other multiplicities is zero and $\mu_{k}=\left<\Delta\right>/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\ll r^{-1}$, the effect of multiplicity on $\mathfrak{R}_{E}$
is small. It starts to rapidly decrease near $k\sim r^{-1}$ 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
$\left<\Delta\right>=r^{-1}$ for different single pathogen infection
probabilities $r$. (Right) The infection risk as a function of the dose scaled
by $r$ ($\left<\Delta\right>r$) for different $k$ and the same fixed
$r=10^{-2}$ ($r^{-1}=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\left<\Delta\right>$ for fixed $r$. This is shown on the right side of Fig 6
for $r=10^{-2}$. For low $k\ll r^{-1}$, the infection risk curves are nearly
identical. For $k\geq r^{-1}$, 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^{-1}$, 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
$\rho_{p,j}$ is large. Large aerosols where $\left<k\right>\gtrsim r^{-1}$
will contribute less to the infection risk than would otherwise be expected
from their resulting average pathogen dose $\left<\Delta_{k}\right>$. While we
must have $M_{c}>\left<k\right>$, $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\times 10^{-3}$ as was done in the
example, Fig 2 shows that this would be important for $d_{0}>15$
$\mathrm{\SIUnitSymbolMicro m}$ for a high viral concentration of
$\rho_{p,j}=10^{11}$ ${\mathrm{cm}}^{-3}$ and $d_{0}>30$ for the lower but
still high viral concentration of $\rho_{p,j}=10^{10}$ ${\mathrm{cm}}^{-3}$.
If we were to consider $r=5.39\times 10^{-2}$ as was also done in the example,
Fig 2 shows that this would be important for $d_{0}>5$
$\mathrm{\SIUnitSymbolMicro m}$ for a high viral concentration of
$\rho_{p,j}=10^{11}$ ${\mathrm{cm}}^{-3}$ and $d_{0}>10$ for the lower but
still high viral concentration of $\rho_{p,j}=10^{10}$ ${\mathrm{cm}}^{-3}$.

Going back to the risk overestimation from ignoring multiplicity in Fig 5,
decreasing $r$ decreases the underestimation in $\tau_{50}$ because the ratio
of the average multiplicity in the larger diameter bins to $r^{-1}$ 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 $\alpha_{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

$q_{p}(t)=\frac{1}{V}\left[\sum_{j=1}^{N_{I}}{\lambda_{I,j}(t)}+\sum_{j=1}^{N_{S}}{\lambda_{S,j}(t)}+\sum_{j=1}^{N_{O}}{\lambda_{O,j}(t)}\right]=\frac{\sigma_{A}}{\left<h\right>}\left<\lambda_{A}\right>_{A}\quad,$
(82)

where $\sigma_{A}$ is the horizontal area density of all individuals and
$\left<h\right>$ is the average height of the environment.

The mean adult breathing rates from sedentary/passive to high intensity
activity ranges between $0.25$ ${\mathrm{m}}^{3}\text{\,}{\mathrm{hr}}^{-1}$
and $3.2$ ${\mathrm{m}}^{3}\text{\,}{\mathrm{hr}}^{-1}$ [14]. For sitting, it
would be hard to get $\sigma_{A}$ to be more than 1 ${\mathrm{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
$\left<h\right>=4$ $\mathrm{m}$, this density limit would yield
$\max(q_{p})\in\left[0.063,0.8\right]$ ${\mathrm{hr}}^{-1}$. If the
environment is poorly ventilated (total ventilation rate $q_{v}+q_{o}+q_{r}$
less than 1 ${\mathrm{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
$\alpha_{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 $\mathrm{m}$ social distancing, the maximum $\sigma_{A}$ are $0.14$ and
$0.080$ ${\mathrm{m}}^{-2}$ respectively. For a typical room height of
$\left<h\right>=4$ $\mathrm{m}$, this density limit would yield
$\max(q_{p})\in\left[0.005,0.11\right]$ ${\mathrm{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}\sim\sigma_{A}$. Since $\beta_{k}\sim N_{I}$ and $\alpha_{C,f}\sim
N_{C}$ but the non $\alpha_{C,f}$ terms of $\alpha$ stay constant, the source
increases faster than the sinks meaning that $n_{k}$ increases and therefore
$\mathfrak{R}$ increases. So, increasing $\sigma_{A}$ with everything else
held constant increases the risk for the susceptible individuals. Thus,
deliberately making $\alpha_{C,f}$ non-negligible is not a viable strategy to
decrease risk. If the $\alpha_{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 $\beta_{I,k}$, the sinks
$\alpha_{C,f}$, and the total dose over time $\mu_{j,k}$. Masks can
substantially improve the total filtering efficiency of the people in
$\alpha_{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
$\alpha_{C,f}$ will have only a small effect on the total sink $\alpha$.
Instead, the main contribution is to reducing $\beta_{I,k}$ and $\mu_{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 ${\mathrm{m}}^{-3}$ after 3 $\mathrm{hr}$ (Fig 3). During Stage 2, an
addition infectious individual has entered the room. The second infectious
individual’s $\rho_{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
$\rho_{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 ${\mathrm{m}}^{-3}$.

The reduction in the average dose $\mu_{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 $\alpha$ is negligible
($\alpha_{C,f}$ is generally small compared to the other sinks) and
$\beta_{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}\sim S_{I,m,out}$ since
$\beta_{I,k}\sim S_{I,m,out}$. Then, $\mu_{j,k}\sim S_{S,m,in}n_{k}\sim
S_{S,m,in}S_{I,m,out}$, which is a quadratic term. Now $\alpha_{C,f}\sim
S_{C,m,in}S_{C,m,out}$ makes the effect stronger (usually only slightly
stronger) than quadratic since it only serves to increase $\alpha$ 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_{0}$, then if exposure starts at steady state,
$\mu_{j,k}\sim Sn_{k,\infty}\sim S\beta_{k}/\alpha\sim S^{2}/(1-cS^{2})$ where
$c\in[0,1)$ is a constant that depends on the relative importance of the
$\alpha_{C,f}$ in the total $\alpha$. In this form, it is easier to see how
$\mu_{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 $\beta_{I,k}$ or reducing $\mu_{j,k}$ plus the effect on
$\alpha_{C,f}$). If only non-susceptible non-infectious individuals wear
masks, there is still a reduction in the dose but it is small since
$\alpha_{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_{0},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_{0},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_{0},t)d_{0})$. 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_{0},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_{0},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_{0},t)\left[1-E_{I,m,out,j}(d_{0})\right]$. These partial
corrections could be done for specific cases (e.g. susceptible individual 2 is
1 $\mathrm{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_{0},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 $\zeta$ and solute composition, and
therefore the same $w(d_{0},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
$\mu_{k,j}$ summed over the individual solutions. This would also be the
solution if $\zeta$ varies in different locations in the respiratory tract
where infectious aerosols are produced for an infectious person. If $\zeta$
changes over time for the sources but the solute composition is constant, then
one could generalize the model to additionally track $\zeta$ (or equivalently
$d_{D}$) and initial diameter at production $d_{0}$ separately.

Another problem is the choice of diameter limits $d_{0}\in[d_{m,k},d_{M}]$ for
each multiplicity. We have neglected the fact that the solute concentration is
much greater for $d_{0}$ 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$ $\mathrm{\SIUnitSymbolMicro m}$. To look at it,
we took the example case and re-calculated it for 23 equal log-width bins
between 100 $\mathrm{nm}$ and 100 $\mathrm{\SIUnitSymbolMicro 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
$\mathrm{\SIUnitSymbolMicro m}$, $54.8$ $\mathrm{\SIUnitSymbolMicro m}$, and
$30.1$ $\mathrm{\SIUnitSymbolMicro m}$. The time step for the numerical
solution had to be reduced to $5\times 10^{-6}$ $\mathrm{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\times 10^{-2}$,
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\times 10^{-3}$,
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