Analyzing the dominant SARS-CoV-2 transmission routes toward an ab initio disease spread model.

Identifying the relative importance of the different transmission routes of the SARS-CoV-2 virus is an urgent research priority. To that end, the different transmission routes and their role in determining the evolution of the Covid-19 pandemic are analyzed in this work. The probability of infection caused by inhaling virus-laden droplets (initial ejection diameters between 0.5 µm and 750 µm, therefore including both airborne and ballistic droplets) and the corresponding desiccated nuclei that mostly encapsulate the virions post droplet evaporation are individually calculated. At typical, air-conditioned yet quiescent indoor space, for average viral loading, cough droplets of initial diameter between 10 µm and 50 µm are found to have the highest infection probability. However, by the time they are inhaled, the diameters reduce to about 1/6th of their initial diameters. While the initially near unity infection probability due to droplets rapidly decays within the first 25 s, the small yet persistent infection probability of desiccated nuclei decays appreciably only by O ( 1000 s ) , assuming that the virus sustains equally well within the dried droplet nuclei as in the droplets. Combined with molecular collision theory adapted to calculate the frequency of contact between the susceptible population and the droplet/nuclei cloud, infection rate constants are derived ab initio, leading to a susceptible-exposed-infectious-recovered-deceased model applicable for any respiratory event-vector combination. The viral load, minimum infectious dose, sensitivity of the virus half-life to the phase of its vector, and dilution of the respiratory jet/puff by the entraining air are shown to mechanistically determine specific physical modes of transmission and variation in the basic reproduction number R 0 from first-principles calculations.

INTRODUCTION

One of the longstanding questions of pandemics involving respiratory droplets is identifying their dominant mode of transmission. The most well recognized pathways for infectious respiratory diseases are (i) the direct contact/inhalation of the relatively larger infectious droplets (>5 µm) commonly known as the droplet mode of transmission, (ii) airborne or aerosol transmission, which is presumed to be caused by inhalation of very small infectious droplets (<5 µm) floating in air, and (iii) contact with infectious surfaces–fomites. For the present Covid-19 pandemic, while the droplet mode of transmission is well established, evidence for aerosol transmission renders identifying the dominant transmission route an intriguing scientific problem with extremely high implications for human health and public policy. On July 9, 2020, the World Health Organization issued a scientific brief stating “Urgent high-quality research is needed to elucidate the relative importance of different transmission routes; the role of airborne transmission in the absence of aerosol generating procedures.” In this paper, we establish a fundamental theoretical framework where the relative strength of the individual transmission routes is analyzed from first-principles with idealizing assumptions. Many biological aspects of the disease transmission including but not limited to effects of immune response are beyond the scope of this paper and will not be addressed here with an exclusive focus on the physical aspects of the disease transmission. Physics is involved in at least four levels in a Covid-19 type pandemic evolution: micro-scale droplet physics, spray/droplet-cloud physics, collision/interaction between the spray/cloud and the susceptible individuals, and deposition and absorption of the inhaled droplets/droplet nuclei. The first three are addressed in this paper at different levels of complexity. We adopt the convention that respiratory droplets (all liquid phase droplets of all sizes, typically 0.5 µm–750 µm, including both airborne and ballistic droplets) cause disease transmission by “droplet or d” route, whereas dried or desiccated droplet nuclei, which, in this paper, refer to the semi-solid/crystalline residue that remains after the droplet liquid evaporates, are responsible for the “dried droplet nuclei or n” route of transmission. Thus, the d route invariably includes droplets less than as well as greater than 5 µm, instead of resorting to the rather arbitrary threshold to distinguish between droplets and nuclei. The reason of our choice is that the distinct thermodynamic phase of the transmission vector (liquid vs solid/semi-solid) is expected to be a much better identifier to delineate the different pathways. Furthermore, the virus survivability within the dried droplet nuclei could be well different from that of the liquid droplet. Small and medium sized droplets do remain airborne after their ejection for substantially long periods of time due to the fact that the droplet size continuously changes, except in highly humid conditions (RH > 85%), due to evaporation. Respiratory droplets are ejected during different expiratory events: breath, cough, sing, sneeze, or talk(ing), when the droplets are ejected with different droplet size distributions. In violent expiratory events like coughing or sneezing, the droplets co-move with a turbulent jet of exhaled air. The trajectories of the jet and the droplets could diverge due to droplet inertia and gravity effects. Nevertheless, experiments by Bourouiba et al. have shown that these droplets can travel rather large distances initially within a turbulent jet, which later transition to a puff or a cloud due to the lack of a continuous momentum source. Depending on the ambient conditions and droplet size, these droplets evaporate at different times. However, while water, the volatile component of the mucosalivary liquid, evaporates, the non-volatile components, salt, protein, mucus, and virus particles present, separate out by crystallization processes. It is to be noted that a homogeneous solution droplet without any admixture or impurities undergoes homogeneous nucleation toward salt crystallization only below a certain threshold relative humidity called efflorescence relative humidity, which is 45% for NaCl at room temperatures. However, in the presence of different components, especially 100 nm sized virus particles, there is no dearth of nucleation sites in infectious mucosalivary liquid. Hence, in this paper, we will consider salt crystallization beyond a certain threshold supersaturation to form droplet residues at all relative humidities. These droplet residues, typically about 10%–20% of the initial droplet diameter are called dried droplet nuclei, remain floating as aerosols and are believed to be responsible for the airborne mode of droplet transmission. While only very few experiments have so far probed the structure of these dried droplet nuclei, the first of its kind work by Vejerano and Marr provided critical insights on the distribution of virus particles inside the dried droplet nuclei. Marr et al. offered mechanistic insights on the role of relative humidity in respiratory droplet evaporation but the question on the role of evaporation, the resulting chemistry inside the dried droplet nuclei on virus survivability persisted. This was explored by Lin and Marr. It was found that virus survivability inside sessile droplets is a non-monotonic function of ambient relative humidity and of course dependent on the specific virus type as well. Questions on the survivability of the SARS-CoV-2 virus inside the dried droplet nuclei from contact free droplets, as it would happen for respiratory sprays, are not yet settled. In this paper, first, we present a model to identify the probability of infection transmission for two different routes d and n by accounting for the corresponding droplet size distribution, viral load, and virus half-life. Next, we briefly present the droplet/nuclei cloud aerodynamics and respiratory droplet evaporation physics. The 1% (w/w) NaCl–water solution is used as a surrogate for the mucosalivary fluid. This is followed by modeling the generalized infection rate constant, which can be used in theory for any kind of expiratory event or for any mode of transmission. The rate constant is then incorporated into a Susceptible-Exposed-Infectious-Recovered-Deceased (SEIRD) model using the formalism of a chemical reaction mechanism. Finally, the Results and Discussion are presented followed by the Conclusions.

At present, widely used epidemiological models for infectious respiratory diseases do not account for the underlying flow physics of disease transmission. In most epidemiological models, the rate constants (parameters of the SEIRD differential equations) that lead to are obtained by fitting available data on the number of new infections. Indeed, these types of epidemiological models have provided immense insights on Covid-19 and the necessity of non-pharmaceutical public health interventions. However, it is to be recognized that the data for a Covid-19 type pandemic are almost always under-reported due to a large number of asymptomatic cases. Furthermore, the actual rate constants and could depend on several physical factors such as temperature, relative humidity, UV-index, and so on. Thus, with changing conditions, parameters obtained from fitting recent past data may not be adequate, standalone, to predict the nature of future outbreaks. Therefore, there is a pressing need to develop a framework to understand and calculate from ab initio calculations while being cognizant of the idealizing assumptions and limitations involved in the process. To the best of the authors’ knowledge, this is the first time that the aerodynamics and thermodynamics of the droplets/nuclei, viral load, half-life, and minimum required viral dose for infection have been systematically accounted for to obtain the rate constants and of a SEIRD model, ab initio.

THE MODEL

Probability of infection by different transmission routes

In this subsection, we estimate the probability of infection by different transmission modes for any expiratory event. Consider an infected person I exhaling a droplet laden jet that quickly transforms into a droplet cloud in the vicinity of susceptible individuals S, as shown in Fig. 1.

FIG. 1.

Schematic of the interaction between S and the droplet/dried droplet nuclei cloud ejected by I, resulting in E. The red ellipse marks the control volume containing the turbulent cloud (colored by red) on average, analyzed for computing the probability of infection .

The instantaneous diameter of the respiratory jet/puff/cloud is given by σ, its velocity with respect to S is given by , while the effective diameter of the hemispherical volume of air inhaled by S, for every breath, is given by σ. σ = (σ + σ)/2. The primary objective of this subsection is to estimate , which denotes the time dependent probability of infection of S for the given expiratory event α and the type of transmission vector β. Thus, α denotes one among breathing, coughing, singing, sneezing, or talking, while β denotes one among droplets, dried droplet nuclei, or fomites. It is to be recognized that are not only functions of time but also dependent on α and β, though their subscripts have been dropped for brevity. As the droplet/nuclei cloud entrains surrounding air, it grows in size with concomitant dilution of the particles inside, thereby reducing . Of course, must also be determined by the viral load and droplet size distribution. could be obtained by solving the transport equations. Instead, here we take a Lagrangian approach of tracking and analyzing the control volume of the air-droplet cloud ejected by I and droplets/droplet nuclei within. At the moment of the onset of the expiratory event denoted by t = 0, a log-normal distribution could be used to describe the probability density function (pdf) of the initial droplet size distribution f of the ejected respiratory spray,where D is the sample space variable of the initial droplet diameter D, and μ and σ are the mean and standard deviation of ln(D), respectively. If N is the total number of droplets ejected for the expiratory event α, the number of droplets within the interval dD is given by Nf(D)dD. Therefore, for a given f and ρ, the viral load in the number of copies per unit volume of the ejected liquid, the cumulative number of virions in droplets between sizes D1 and D2 is given by

For the SARS-CoV-2 virus, Wölfel reported the average viral load in sputum to be ρ = 7 × 106 copies/ml, while the maximum is given by ρ = 2.35 × 109 copies/ml. Utilizing Eq. (2), we can define (time dependent number of virions inhaled from droplets). For an infection to occur, some non-zero number of active virions must be found in the droplets present in the total volume of air inhaled. The maximum time for S to cross the volume with diameter σ is given by t = (σ + σ)/V, while the number of breaths per unit time is N ≈ 16/60 s−1 and the volume inhaled per breath is . Therefore, the total volume of air inhaled while crossing the respiratory cloud is . The fraction of virion population surviving within the droplets or dried-droplet nuclei at time t is given by ψ(t) and can be assumed to decay as . is the half-life of the SARS-CoV-2 virus in d or n. While the half-life of the SARS-Cov-2 within aerosols, in general, could be estimated from Refs. 27 and 28, the distinction of the half-life of the virus within airborne droplets or dried droplet nuclei is not yet available to our knowledge. While this information is indeed most critical, in view of its absence, for the present work, we will mostly assume ψ(t) = ψ(t) = ψ(t) and , unless specifically mentioned. At typical indoor conditions, T = 21.1 °C and RH = 50% and with UV index =1 on a scale of 10, min. Accounting for these, is given by the following equation:Here, D1(t) and D2(t) are the minimum and maximum of initial droplet diameters, respectively, available in the droplet cloud after time t, as shown in Fig. 2. The non-linearity of t vs D (in the log–log plot of Fig. 2) occurs due to the phase transition of the droplet population. Beyond τ, all droplets have been converted into dried droplet nuclei. At time t, droplets with D < D1(t) have evaporated and have been converted to dried droplet nuclei, while D > D2(t) have escaped by gravitational settling and have been converted to potential fomites. Clearly, as σ increases with time, decreases due to dilution effects and also because droplet numbers are depleted by evaporation and settling. Since the maximum evaporation time of the airborne droplets , the effect of virus half-life on is negligible. The details of the methodology to derive the parameters concerning the respiratory jet and droplet dynamics from the conservation principles of mass, momentum, energy, and species could be found briefly in Subsections II B and II C, respectively. Further details could be found in the work of Chaudhuri et al.

FIG. 2.

Wells curves modified by accounting for droplet cloud aerodynamics, thermodynamics, and droplet desiccation for 1% NaCl–water droplets at T = 21.44 °C and RH = 50%. After time t, the droplets with initial size D < D1(t) have been converted to dried droplet nuclei; those within D1(t) ≤ D ≤ D2(t) are in the liquid droplet state; D > D2(t) have settled and could be potential fomites. τ = 22.87 s. The red, blue, and yellow shaded regions denote the regimes of droplets, dried droplet nuclei, and fomites, respectively.

It is essential to note that in many diseases, uncertainty exists over transmission routes. For example, in the case of Covid-19, the transmission by dried droplet nuclei is not certain. As such, we do not know for sure if the SARS-CoV-2 virus survives within dried droplet nuclei and remain infectious, though there is evidence that SARS-CoV-2 and some other virus do survive quite well inside dried droplet nuclei. Even if they do, their half-life and infection potential could be different with respect to those inside droplets. Thus, as would be shown later, it is essential to define different rate constants for different transmission modes. Since the dried droplet nuclei are a product of the droplets, itself, the two routes are highly coupled and are not independent.

For the droplet nuclei, is given byD(t) = D1(t) if t < = τ, and D(t) = D2(t) if t > τ. decreases with time due to the increase in σ with time, i.e., the dilution effect as well due to virus half-life. decreases with time due to the increase in σ with time, i.e., the dilution effect as well due to virus half-life.

The generalized probability of infection as a function of infectious dose inhaled while crossing the droplet/nuclei cloud , can now be expressed as

The total probability of infection for the expiratory event α could be defined as

The form of Eq. (6) is based on the dose response model by Haas, which has been used by Nicas, Sze To et al., and many other authors to calculate the infection probability. Mathematically, it is also similar to the Wells–Riley equation used by Buonanno et al. to assess the aerosol risk of SARS-CoV-2 during talking and breathing. However, in contrast to these works, here, droplet cloud aerodynamics (Subsection II B) coupled with the detailed droplet evaporation-nuclei production mechanism (Subsection II C) and droplet settling dynamics are utilized in a semi-analytical framework to calculate the time varying inhaled virion number and corresponding probability of infection, probably for the first time. Eventually, as shown later, this framework will be used to calculate the basic reproduction number . As such, the form of this equation is also validated by the results from the work of Zwart et al., where r is a constant for a particular virus. For this paper, we will use r = 0.5 such that inhaling at least ten virions by d and/or n route would result in an infection probability , unless specifically mentioned. In the absence of this exact r for the SARS-CoV-2 virus at the time of writing this paper, this is an educated guess. Hence, r = 0.05 and r = 0.005 corresponding to minimum infectious doses of 100 and 1000 virions, respectively, will also be eventually explored near the end of the paper.

Aerodynamics of droplets and nuclei

The droplets when ejected during respiratory events mostly follow the volume of exhaled air. Due to continuous entrainment, the exhaled air volume grows in diameter, and as a result, its kinetic energy decays with time. Bourouiba et al. identified that for a short duration, the exhaled droplets evolve inside a turbulent jet, which transitions to a puff at later stage. Since the respiratory droplets or the dried nuclei experience aerodynamic drag, it is essential to identify the evolution of the surrounding jet or puff. Based on the literature of transient turbulent jets and puffs, the following evolution equations for the axial location, velocity, and radial spread could be used:andwhere the subscripts j and pf denote the jet and puff, respectively, R0 and U0 are the average radius and axial velocity at distance x0, and K is a characteristic constant for turbulent jet and is reported to be 0.457. At the inception of the respiratory event (t = 0), the jet is assumed to have a velocity of U = 10 m/s and a radius of R = 14 mm—the average radius of human mouth. For analytical tractability, we assume that all droplets of all sizes are ejected at time t = 0 and would not consider time variation in the ejection of the droplets. This is a safe assumption since the expiratory event like cough lasts less than a second and the turbulence of the jet and the air entrained will rapidly disperse the ejected droplets into the jet/puff in any case. The characteristic constants for puff are a ≈ 2.25 and m = (xa)/(3R). Since the continuous ejection of air from mouth lasts only for the duration of the corresponding respiratory event, the jet behavior persists only for this period. Beyond this time (≈1 s), the puff behavior is observed. Hence, the velocity and the radial spread of the air surrounding the exhaled droplets will beThe horizontal displacement (X) of the exhaled droplet and its instantaneous velocity (U) due to the drag can be solved withHere, ρ and ρ are the vapor and liquid phase densities, respectively, R is the instantaneous radius of the droplet, and C is the drag coefficient. We can assume C = 24/Re for the gas phase Reynolds number, Re = (2ρ|U − U|R)/μ < 30. Re for the respiratory droplets are typically less than 0.1. At the time of ejection (t = 0) from respiratory cavities, the droplets are assumed to have a velocity (U) close to that of the surrounding air (U), and hence, Re ≈ 0. Thus, for t = 0, we use U − U = 0.01U. The time for a droplet with initial diameter D [given by t(D)] to fall a height of 1.8 m is calculated using Stokes’s settling velocity. Mathematically, t is obtained by the following equation:By solving Eqs. (7)–(10) over the droplet and nuclei lifetime, the axial distance traveled by them, X, which is the distance of the center of the cloud, can be evaluated. As the velocity of the individual droplets approaches the surrounding gas velocity within a very short time, we assume the absolute instantaneous velocity of the droplet/nuclei cloud is given by V = U from Eq. (9). Since the droplets and nuclei are dispersed within the jet/puff, the diameter of the droplet cloud ejected by I can be approximated as twice the radial spread of the exhaled air, σ(t) = 2R(t). The puff will continue to grow in a large volume without obstruction. For a small, nearly cubic room with strong mixing and poor ventilation, a constant σ for t ≥ t could be estimated from the volume of the room V using . This implies a volume filling state of the ejected cloud and a homogeneously mixed state of the aerosols. is the time at which grows to .

The exhaled volume of air is initially at a temperature (T = 33.25 °C) and vapor mass fraction (Y1, corresponding to RH0 = 71.6%) different from the ambient. The values mentioned are averaged quantities measured over several subjects according to Mansour et al. The instantaneous temperature and vapor mole fraction that the droplet would encounter as its own ambient are the temperature [T(t)] and vapor mass fraction [Y1,(t)] of this volume of air during its evolution. This can be expressed with the following scaling relation:where ΔT = T − T and ΔY1, = Y1, − Y1,.

Droplet evaporation

In this paper, we use 1% NaCl–water droplets as the model respiratory droplet and adopt a slightly revised evaporation model (with respect to that presented in Ref. 29) for predicting the droplet evaporation time t. It is to be recognized that the droplets are surrounded by exhaled air volume, as described in Subsection II B, and hence, it serves as the “ambient condition” for the droplet. The evaporation mass flux for quasi-steady state conditions is given byHere, is the droplet mass loss rate due to evaporation, R is the instantaneous droplet radius, ρ is the density of water vapor, and D is the binary diffusivity of water vapor in air, and α is the thermal diffusivity of surrounding air. B = (Y1, − Y1,)/(1 − Y1,) and B = C(T − T)/h are the Spalding mass transfer and heat transfer numbers, respectively. Y is the mass fraction with the numerical subscripts 1, 2, and 3 denoting water, air, and salt, respectively. Additionally, the subscripts s, g, and ∞ denote the location at droplet surface, surrounding gas, and at very far field ambient, respectively. h and C are the specific latent heat of vaporization and specific heat of the droplet liquid, respectively. Unlike in a pure water droplet, the vapor pressure at the surface of droplets with non-volatile dissolved substances as in respiratory droplet/salt solution droplets could be significantly suppressed. Raoult’s law provides the modified vapor pressure at the droplet surface for the binary solution, P(T, χ1,) = χ1,P(T), where χ1, is the mole fraction of evaporating solvent (here water) at the droplet surface in the liquid phase and χ1, = 1 − χ3,. The far field vapor concentration, on the other hand, is related to the relative humidity of the ambient. Considering the effects of Raoult’s law and relative humidity, the vapor concentrations at the droplet surface and far field can be expressed aswhere M1 and M2 are the molecular weights of water and air, respectively. Instantaneous Y1, is evaluated from Eq. (12). The latent heat required for evaporation is provided by the droplet’s internal energy and/or surrounding ambient. It has been verified that the thermal gradient in the liquid phase is rather small. Therefore, neglecting the internal thermal gradients, T is obtained from the energy balance,where T is the instantaneous droplet temperature, and are the instantaneous mass and surface area of the droplet, ρ is the density of the binary mixture of salt (if present) and water, and k is the conductivity of air surrounding the droplet. is the thermal gradient at the droplet surface and can be approximated as (T − T)/R. Due to the continuous loss of water, the solution would become supersaturated in most occasions leading to the onset of crystallization. The crystallization kinetics is modeled with a one-step reaction. The validation of the model (1% NaCl–water solution) with saliva droplet experiments (average of three runs) from a healthy subject is shown in Fig. 3. The experiments were performed in a contact free condition in an acoustic levitator at T = 28 °C and RH = 41%. The reader is referred to Refs. 29 and 46 for details of the experimental configuration. Here, a difference of up to 15% could be found for the different stages of droplet drying, between the model prediction and the saliva droplet drying curve from experiments. It should be noted that human saliva contains mucus and varieties of salts and electrolytes along with compositional variations, which are difficult to model very accurately.

FIG. 3.

Comparison of the model output for 1% NaCl–water droplets with human saliva droplet experimental data averaged over three runs. D (normalized with D shown in the left ordinate) is the effective diameter of the droplet accounting for both the solute and solvent, while D1 (shown in the right ordinate) is the effective diameter only accounting for the solvent (water) mass of the droplet.

Ab initio rate constants for SEIRD model

With the methodology to calculate the probability of different transmission routes identified, we proceed to evaluate the respective “rate constants.” Widely used epidemiological models do not account for the flow physics of disease transmission. Within the framework of the well known SIR-model, the model constants proposed by Stilianakis and Drossinos included overall droplet cloud features like the number of droplets per unit volume of the cloud but did not include crucial physics like cloud aerodynamics, evaporation, or crystallization dynamics that lead to droplet-nuclei formation. As such, these control the time evolution of the droplet cloud, and as shown later, the spatio-temporal evolution of the cloud and the constituent droplets plays a major role in determining the critical rate constants of the problem. A model connecting the macro-scale pandemic dynamics with the micro-scale droplet physics accounting for droplet-cloud aerodynamics, evaporation, and crystallization physics has been recently presented by Chaudhuri et al. Drawing inspiration from the well known molecular collision theory of reactions due to collisions, a chemical reaction mechanism was obtained where three elementary reactions described the pandemic evolution. Adopting the notations of the SEIRD model, one of the reaction rate constants that determined the conversion of a susceptible individual S to an exposed individual E, upon contact with the droplet cloud ejected by the infectious person I, was denoted by k1 (or k1, as opposed to the new rate constant to be defined here) and was called infection rate constant. This k1, was modeled using the molecular collision theory. From Chaudhuri et al., the expected number of collisions per unit time between S and of I is given by resulting in the infection reaction I + S → I + E. n and n are the number of infected I and susceptible S people in unit volume. The infection rate constant of this reaction is then given byσ is the jet/puff diameter, which is also assumed to be the diameter of the droplet cloud. is the relative velocity of the droplet cloud with respect to S, while τ is the droplet lifetime. V can be approximated as the preferred walking speed, which according to Refs. 49 and 50 equals to 1.3 ± 0.3 m/s, and hence, for the current study, we will assume V = 1.3 m/s. t is the average time period between two expiratory events. t is calculated as t = 3600 × 24/Nexp where Nexp is the average number of infecting expiratory events per person per day. We assumed N = 3 from the measured coughing frequency of 0–16 in normal subjects.

Equation (16) is limited by several simplifying assumptions. In this paper, we develop a comprehensive model beyond these limitations, which is also rendered capable of delineating the relative dominance of the different disease transmission pathways. In particular, in the following, we derive a new rate constant accounting for (i) the finite viral loads, finite viral lifetime, and the corresponding probability of infection, (ii) the evolution of the collision volume with time, (iii) transmission by droplets of any sizes and the corresponding dried droplet nuclei, and (iv) the inhomogeneity of infection spreading. Furthermore, in this paper, we generalize the infection rate constant equation to account for transmission by any expiratory event. To that end, a generalized reaction mechanism that accounts for different modes of infection transmission as well as different forms of expiratory events is presented. This is followed by a comprehensive modeling of the individual infection rate constants, following which we arrive at an overall infection rate constant. In view of the above discussion, the basic reaction mechanism of Ref. 29 could be generalized to a comprehensive one where almost all possible expiratory events and modes of transmission could be included to yield

In [R1], α varies over different expiratory events, namely, breath, cough, sing, sneeze, and talk, while β varies over different modes of transmission, namely, droplet, droplet nucleus, and fomite. Thus, [R1] essentially represents 15 reactions. The rate constants of these individual reactions are defined in Table I.

TABLE I.

Infection rate constants for different expiratory events and modes of transmission.

k1,αβ	Droplet	Nucleus	Fomite
Breath	k1,bd	k1,bn	k1,bf
Cough	k1,cd	k1,cn	k1,cf
Sing	k1,gd	k1,gn	k1,gf
Sneeze	k1,sd	k1,sn	k1,sf
Talk	k1,td	k1,tn	k1,tf

Here, each of the parameters should be obtained for the respective combination of α, β. Including [R2] and [R3], in total, there are 15 + 2 = 17 reactions to be included in a complete model. As such, further granularity could be added by adding a location parameter γ. In that case, we can have k1,, where γ could represent home, school, office, transport, restaurants, and so on. In this paper, we will only consider two selected transmission modes: cough-droplets and cough-dried droplet nuclei with the rate constants k1, and k1,, as shown in Table I. These are expected to play the more dominant roles in disease transmission. However, the approach here could be used for any other transmission routes too, with the corresponding droplet size distribution. For Covid-19, fomites are being considered as a secondary source of infection and need to be dealt separately.

Individual rate constants will allow us to delineate the different modes of transmission on average. Furthermore, the definition of individual rate constants enables quantitative investigation of the relative dominance of each mode of transmission. This constitutes one of the major goals of this paper. The infection rate constants are generalized by inclusion of the probability for infection , averaging the collision volume over a characteristic time alongside including the dried droplet nuclei mode of transmission, in addition to the droplet mode of transmission. The revised rate constant for any expiratory event α, vector of transmission β, and location γ is given by the following equation:

Specifically, by utilizing Eqs. (3) and (5), the rate constant for the droplet mode of transmission d ejected during any expiratory event α could be calculated as

For the droplet nuclei, Eqs. (4) and (5) yield the rate constant k1,,

Note that we introduced a new parameter ψ(t) to calculate , which denotes the fraction of the infectious virion population active within the dried droplet nuclei population at time t. As t ∞, ψ(t) 0. Thus, the integration is performed by up to about —the largest evaporation time of the droplet set considered. The details on the survivability of specific SARS-CoV-2 inside dried droplet nuclei are not known. Hence, for now, we will assume that ψ(t) is independent of d and n, except when we will estimate its sensitivity in specific cases. According to the reaction mechanism given by [R1], [R2], [R3], E is formed by several parallel pathways. Therefore, the corresponding rate constants become additive. Hence, the location (γ) dependent infection rate constant can be defined as

While the rate constant k1, is derived from first-principles, it still results in the same infection rate constant for a given set of ambient temperature T, RH, and population density. As shown by Lloyd-Smith et al., the individual infectiousness distribution around the basic reproduction number is highly skewed. This suggests that a small fraction of infected individuals (superspreaders) are responsible for a large number of infections. Hence, the final challenge of this modeling effort is to include this effect. Such “superspreading” events could be results of (i) high local population density n, (ii) highly mobile infected individuals, and (iii) most importantly, high viral loading of the ejected respiratory droplets ρ. We will see that large viral loading ρ = ρ leads to very high infection probability, which would lead to large k1,. Since the rate constant is directly proportional to n, its effect is understandable. Thus, the effect resulting from mobility needs to be accounted.

Understanding and modeling human mobility at both the individual level and the population level have garnered recent interest. See a recent review by Barbosa et al. for a detailed exposition on this topic. Kölbl and Helbing used the statistical data of the UK National Travel Surveys collected for 26 years by the Social Survey Division of the Office of Population Census and Surveys to arrive at a generalized distribution of human daily travel behavior. They showed that for different modes of transport i ranging from walking, cycling, car driving, and so on, the travel time t normalized by the average travel time for the corresponding mode of travel and defined as , a common distribution for τ, irrespective of the mode of transport could be obtained. This was also argued from an energy point of view, where ; kJ per person per day—the average travel energy budget of the human body according to Ref. 54. In any case, the pdf, , after dropping the i given its universality isThe following constants were provided: for the universal curve.

Given two infected people I, it is reasonable to expect that the one with the higher mobility has more chance to infect others since they have greater exposure to the population and can infect people at different locations, all other conditions remaining fixed. Therefore, we can assume that the final infection rate constant should be proportional to τ.

The corresponding infection rate constant summed over all possible types of expiratory event α, transmission mode β, and location γ is, thus, given byClearly, k1 is now a function of two random variables τ, n, and the extreme individual realization of each could correspond to the superspreading events.

Finally, the average, overall infection rate constant k1 is given byThis is because and that ∑k1, is independent of τ. In the rest of the paper, we will mostly focus on this ensemble averaged rate constant k1. With the framework established, the individual effects of mobility and population inhomogeneity could be taken up in future works.

From the reactions [R1], [R2], [R3], we can obtain the set of ordinary differential equations of the SEIRD model, which would govern the evolution of [S], [E], [I], [R], and [], where the rate constants appear as respective coefficients. Here, S denotes susceptible, I denotes infected, E denotes exposed, R denotes recovered, and denotes deceased. The square brackets denote the number of the particular population type normalized by the total population, e.g., [I] = n/n,

RESULTS AND DISCUSSION

From the measurements by Duguid, the droplet size distribution from cough could be described using a log-normal distribution. The initial distribution f and the number of virions present in each droplet size for the average viral load ρ = 7 × 106 copies/ml are shown in Fig. 4(a). Clearly, assuming a uniform ρ, larger droplets can contain multiple virions. Thus, it is imperative to model the droplet evaporation and the jet spreading dynamics to know their eventual distributions. The total number of droplets ejected is 5000. Figure 4(b) shows the time evolution of the droplet number distribution as a function of instantaneous diameter D. The shift of the distribution to the left, i.e., toward smaller D, is an effect of evaporation. In addition, the right branch of the number distribution gets eroded due to settling of the larger droplets. At the conditions of interest, T = 21.44 °C and RH = 50%, the modal diameter of the droplet nuclei is 2.3 µm at t = 1050 s starting from an initial modal diameter of 13.9 µm at t = 0 s. Note that here, by the modal diameter, D corresponding to the peak of the histogram shown in Fig. 4(b) is referred. Interestingly, since small droplets evaporate fast, the left branch (small sizes) of the distribution shifts fast to further smaller sizes. A droplet with initial diameter D = 10 µm is reduced to D = 1.96 µm within t = 0.42 s. As such, for the entire droplet set, a modal diameter of 2.7 µm, which is within 22% of the final modal diameter, is achieved within t = 1 s or within a distance of X = 1.8 m from the origin of the respiratory jet. X denotes the distance of the center of the respiratory jet/puff (with a diameter of σ) from its origin. Within t = 10.6 s, X = 2.88 m, the droplet size distribution is very close to the final distribution. Due to this sharp reduction in droplet size due to evaporation (for RH < 85%) combined with settling of large droplets, practically, for most of the time, the disease appears to be transmitted by droplets/nuclei of instantaneous diameter less than 10 µm, the most probable instantaneous diameter being between 2.14 µm and 2.7 µm. However, it is to be noted that this diameter could be 5–6 times smaller than the initial ejected diameter of the droplet D. In a viewpoint article, Fennelly reported that for most respiratory infections, the smaller droplets (<5 µm and collected at a finite distance from the origin of the respiratory spray) were found to be pathogenic. Furthermore, Chia et al. reported PCR (Polymerase Chain Reaction) positive SARS-CoV-2 particles with sizes >4 µm and also between 1 µm and 4 µm from air samples collected. Thus, our results appear to be consistent with these clinical research findings.

FIG. 4.

(a) Probability Density Function (PDF) f of droplet diameter for cough and the number of virions N as a function of the initial droplet size D at t = 0. Average and maximum viral load ρ = 7 × 106 copies/ml and ρ = 2.35 × 109 copies/ml of SARS-CoV-2, respectively, are assumed from Ref. 25. (b) Droplet and/or dried droplet nuclei size distributions at different time t and distance X from the origin of the respiratory jet for T = 21.44 °C and RH = 50%. The curves at the last three time instants end abruptly due to the loss of droplets due to settling beyond that particular diameter. For t > τ = 22.87 s or X > 3.52 m, all airborne droplets have been desiccated to the corresponding droplet nuclei. The size distribution remains invariant for t > τ.

Next, we analyze the time varying infection probability. Interestingly, in Fig. 5(a), at T = 21.44 °C and RH = 50% for t > 1 s, the total probability of infection scales as for droplets and dried droplet nuclei. This is a combined effect of droplet evaporation, virus decay, and dilution due to the entrainment of fresh air within the jet/puff, the diameter of which increases initially as t1/2 and then as t1/4, respectively. After the droplets evaporate, the decay of the infection probability for the dried droplet nuclei slows down with respect to their droplet predecessors. This is because, while the infection probability decay for droplets is due to evaporation, settling, and dilution, the probability decay due to nuclei is due to only dilution and finite virus lifetime. Here, we are considering a very large, poorly ventilated indoor space, such as a shopping mall or a conference center, with a large number of occupants. It is to be noted that the dilution effect is arrived with the assumption that all people are in motion, but their motion do not affect the cloud aerodynamics. In reality, such motion could lead to increased turbulence and mixing, resulting in further dilution. Therefore, the increase in σ and the decay of the probability of infection could be faster in reality. An interesting feature in the original situation of interest (very large space) is that does not continuously follow after all the droplets evaporate. There is an accumulation of the droplet nuclei due to the evaporation of smaller droplets beforehand leading to a small jump in at τ. Indeed, the overall probability by Eq. (6) decreases without any discontinuity. From Fig. 5(b), we find that for X > 2 m, the overall probability decreases as justifying the necessity of social distancing. However, only after about 5 m.

FIG. 5.

Probability of infection for droplet route d and dried droplet nuclei route n, as well as total probability for (a) as a function of time measured from the instant of the beginning of the expiratory event and (b) as a function of distance measured from the location of the origin of the expiratory event along the center of the jet/puff trajectory. T = 21.44 °C and RH = 50%. The bold lines represent ρ = 7 × 106 copies/ml with , where min. The gray shaded region denotes the lower limit and upper limit , respectively, with .

The effect of preventing ejection of droplets beyond particular initial sizes on is examined in Figs. 6(a) and 6(b) as a function of time and distance. For D = 50 µm, an infection probability of 0.6 is obtained for t → 0, suggesting droplets with D < 50 µm are slightly more responsible for infection, at all times, for the conditions under consideration than their D > 50 µm counterparts. This trend continues until τ when all airborne droplets evaporate. This is qualitatively consistent with the exposure analysis and results of Chen et al., who considered the dispersion and evaporation of water droplets with size distribution from Duguid. However, when D = 10 µm, the corresponding probability of infection is very small, suggesting that for the average viral loading, at early times, the droplets of initial diameter 10 µm < D < 50 µm are the most lethal in terms of their probability to infect. However, while they infect, their diameters are substantially smaller.

FIG. 6.

Probability of infection for droplet route d at different cutoff droplet sizes (implying no droplets beyond that size) for (a) as a function of time measured from the instant of the beginning of the expiratory event and (b) as a function of distance measured from the location of the origin of the expiratory event along the center of the jet/puff trajectory. T = 21.44 °C and RH = 50%. ρ = 7 × 106 copies/ml with , where min.

at T = 10 °C and RH = 20% is shown in Figs. 7(a) and 7(b) as a function of time and distance. In comparison to the previous case, here, the droplet survives longer due to lower temperature, while droplet nuclei induce higher due to longer virus half-life. Thus, at lower temperatures, the higher infection probability could be expected. However, in both cases, at short time and distance from the expiratory event, droplets (both small and large) dominate transmission, and only after most droplets evaporate, the dried droplet nuclei route is significantly activated. While the transmission probability by dried droplet nuclei is always lower than that by droplets, their lifetime is theoretically infinite as opposed to the finite lifetime τ of the droplets. Hence, despite their low instantaneous probability of infection, cumulatively, they contribute significantly, assuming that the virus remains infectious for significant times within the dried droplet nuclei. If so, as will be shown below, the persistent dried droplet nuclei appear to be a major transmission mode of the virus. It is to be recognized that these results were obtained with average viral loading ρ = 7 × 106 copies/ml. If we consider ρ = 2.35 × 109 copies/ml, does not decay from the maximum fixed value of 1 until from about 100 s or from 5 m from the origin of the respiratory jet, along the center of the jet. Thus, it is expected that such a kind of viral loading could infect a large number of S potentially leading to a superspreading event.

FIG. 7.

Probability of infection for droplet route d and dried droplet nuclei route n (a) as a function time measured from the instant of the beginning of the expiratory event and (b) as a function of distance measured from the location of the origin of the expiratory event along the center of the jet/puff trajectory. T = 10 °C, RH = 20%, and ρ = 7 × 106 copies/ml.

Using thus obtained, we can evaluate the corresponding rate constants for d and n using Eqs. (18) and (19), respectively. These rate constants for ρ = 7 × 106 copies/ml are presented in Table II for four cases IA, IB, IC, and II. Cases I(A–C) correspond to T = 21.44 °C and RH = 50%, while case II represents T = 10 °C and RH = 20%. The population density in both cases is assumed to be 10 000 people/km2. In all cases, homogeneous mixing is assumed without any social distancing or lockdown. In all cases, k1, > k1,. Case IA represents no restriction and clearly high rate constant values are attained in this case. Case IB represents a hypothetical situation where the ejection of all droplets with D > 10 µm is restricted. This is hypothetically possible by stringent enforcement of population wide usage of ordinary face-masks without any exceptions. Furthermore, using k1, and k3, we can define the basic reproduction number . The calculated k1, and could be found in Table II. A very interesting trend emerges between cases IA and IB. We find that if the ejection of droplets even beyond 10 µm could be completely prevented, drops from 4.22 (0.33, 5.63) for case IA to 0.048 for case IB. For case 1A, the numbers in the brackets denote for the lower limit and upper limit , respectively. between cases 1A and 1B represent a two order of magnitude difference, and for at case IB, no outbreak is possible. The bifurcation point is attained for the critical droplet size D = 27 µm. This is shown in Table II as case IC. The implication is that preventing ejection of droplets with the initial size beyond 27 µm would just prevent the outbreak. Of course, it is to be recognized that we are only considering cough as the mode of droplet ejection alongside many idealizing assumptions. Furthermore, these results were arrived at with the average viral load ρ = 7 × 106 copies/ml. If we consider the maximum reported viral load ρ = 2.35 × 109 copies/ml with free mixing among I and S, , indicating a superspreading event. As such, it could be a combination of high mobility and large viral loading of I, which could lead to a superspreader.

TABLE II.

Calculated infection rate constant values for different modes of transmission for coughing with and without mask at typical indoor conditions. Cases IA, IB, and IC: T = 21.1 °C and RH = 50%. Case II: T = 10 °C and RH = 20%. For all cases, ρ = 7 × 106 copies/ml.

Case	Condition	k1,cdγ	k1,cnγ	R0,c
IA	Cough, no mask	0.0182	0.2743	4.2219
IB	Cough, Ds,0 = 10 µm	1.58 × 10−5	0.0033	0.0476
	cutoff mask for all
IC	Cough, Ds.0 = 27 µm	6.23 × 10−4	0.0686	0.9981
	cutoff mask for all
II	Cough, no mask	0.0229	0.3519	5.4087

calculated at the average and maximum viral loading ρ = 7 × 106 copies/ml and ρ = 2.35 × 109 copies/ml, respectively, for different droplet size cutoffs at T = 21.44 °C and RH = 50% is shown in Fig. 8. The cutoff D means that all droplets with sizes D > D are prevented from ejecting. Figure 8 also shows the sensitivity of the assumption on the results. In Fig. 8, the lower and upper limits represent the conditions and , respectively. If all droplets are allowed to be ejected at average viral loading, for , , while for , with the base for the typical indoor conditions assumed above. Clearly, the change in the lower limit of is much more sensitive than its upper limit. This is because even if the viral lifetime is much longer, dilution reduces infection probability. However, with ρ = 2.35 × 109 copies/ml for , , while for , around the base case of , all other conditions remaining the same. Interestingly, with D = 10 µm, reduces by a factor of 40 with respect to no cutoff condition. However, D = 5 µm reduces by another factor of 18 with respect to D = 10 µm condition for the base cases. At this condition, . For both viral loadings, the maximum and the averaged blocking droplets D ≥ 5 µm can theoretically yield . All the results, so far, have been obtained with r = 0.5, which implies a minimum infectious dose of 10 virions. Figures 9(a) and 9(b) show the corresponding for r = 0.05 and r = 0.005, respectively. These imply minimum infectious doses of 100 and 1000 virions, respectively. While the qualitative trend is similar, indeed, for these two cases are much lower in comparison to r = 0.5. As such, it seems likely that the minimum infectious dose of SARS-Cov-2 is . While the base for the average viral loading, r = 0.5 and no cutoff is consistent with that of the reported values for Covid-19, the order of magnitude larger values of obtained at maximum viral loading should be viewed in the context of superspreading events.

FIG. 8.

Comparison of for different conditions based on two different viral loadings ρ = 2.35 × 109 copies/ml and ρ = 7 × 106 copies/ml and if the ejection of droplets beyond the specified cutoff sizes is prevented. In each box, the red line denotes the equal half-life condition , irrespective of the phase. The lower and upper limit corresponds to and , respectively. All data are obtained at T = 21.44 °C and RH = 50% and with r = 0.5 implying a minimum infectious dose of 10 virions.

FIG. 9.

Comparison of for different conditions based on two different viral loadings ρ = 2.35 × 109 copies/ml and ρ = 7 × 106 copies/ml and if the ejection of droplets beyond the specified cutoff sizes is prevented. In each box, the red-line denotes the equal half-life condition , irrespective of the phase. The lower and upper limit corresponds to and , respectively. All data are obtained at T = 21.44 °C and RH = 50% and with (a) r = 0.05 implying a minimum infectious dose of 100 virions and (b) r = 0.005 implying a minimum infectious dose of 1000 virions.

Using the governing equation (24), the evolution of the pandemic for average viral loading and r = 0.5 for case IA is presented in Fig. 10(a). The growth rate of the infected population for case IA and case IB is shown in Fig. 10(b) with the assumption that the usage of face masks for the entire I population (which would practically be required for the entire population) is implemented after a fixed time from the onset of the outbreak. The effectiveness of facemasks has been mechanistically proven. As expected, in this SEIRD model with ab initio infection rate constants, the effect of the usage of masks is almost immediate since the assumed latency period is only one day. Leffler et al. analyzed Covid-19 data from 198 countries to conclude that government policies on mask wearing significantly reduced mortality. A significant feature of the results of this paper is that, while they are computed mechanistically from first-principles with assumptions and limitations, they produce physically meaningful outcomes.

FIG. 10.

SEIRD-droplet/nuclei model output with computed rate constants mentioned in Table II for case IA. No lockdown and free mixing are assumed. (a) Case IA, no mask usage and (b) case IB, all I wear masks from the (i) 100th day and (ii) 30th day of the onset of the pandemic that prevents the ejection of D > 10 µm droplets with ρ = 7 × 106 copies/ml. Calculated rate constants are shown in Table II for cases IA and IB. The masks for this case are assumed to prevent the ejection of all particles greater than 10 µm.

CONCLUSIONS

Current epidemiological models for infectious respiratory diseases do not account for the underlying flow physics of transmission. This work presents a SEIRD model developed from the mechanics of respiratory disease transmission, including but not limited to known SARS-CoV-2 virus characteristics, thermodynamics, transport, and aerodynamics of respiratory droplets/nuclei cloud leading to its interaction with a susceptible population. Starting with a well known cough droplet size distribution, we derived the time dependent probability of infection and infection rate constants for different routes accounting for viral load, virus stability, respiratory droplet cloud aerodynamics, evaporation, and crystallization for poorly ventilated conditions. Assumptions are inherent in all models. While the assumptions pertaining to the present model have been described throughout the paper at their respective points of discussion, for the sake of clarity, we revisit the three major model assumptions here. (1) Infectious respiratory droplets contain water, electrolytes, mucus, enzymes, and virus. Here, we model it using 1% NaCl–water droplets, thereby accounting for the first and second major components of the mucosalivary fluid. The evaporation characteristics are well modeled, as shown in Fig. 4. Viral load is assumed to be uniform across all droplet sizes. (2) The droplets are ejected in a turbulent cough jet, which transitions to a puff in a poorly ventilated, quiescent, large indoor space. Typical droplet size distribution, mean thermodynamic, and aerodynamic state of such jets and puffs are considered in which the droplets evaporate until their desiccation. While the results only reflect the cough jet and its infection dynamics, the model could be easily extended to account for similar respiratory jets. We do not consider the effect of turbulence or buoyancy on droplet evaporation or in mixing of the aerosols. The internal air-circulation effect has not been considered since it would be location specific. (3) Infection occurs when a certain number of virions are inhaled upon collision of susceptible individuals with the respiratory puff. The collision frequency combined with probability of infection yields the rate constants for the diseases spread ordinary differential equations. The indoor space volume considered is much larger than the puff volume over the time considered. However, the model could be easily extended to a situation where the indoor space is filled with the aerosols.

Considering a well known cough droplet size distribution, most number of droplets have an initial diameter of D = 13.9 µm, but within 1 s of their ejection, most number of droplets of the same set gets reduced to a diameter of 2.7 µm due to evaporation, accounting for the thermodynamic state of the exhaled air, at typical, air-conditioned yet quiescent indoor space. For the average viral loading, at early times, the droplets of initial diameter 10 µm < D < 50 µm are the most lethal in terms of their probability to infect. However, while they infect, their diameters could be 5–6 times smaller. Indeed, for most of the time, infection is spread by inhalation of small, airborne droplets or their desiccated nuclei. While the instantaneous probability of infection by droplets is significantly larger than its dried nuclei in the short time and range, the much longer persistence of the dried nuclei results in its stronger relative contribution to the infection rate constant, under the assumption that the virus half-life is independent of the phase of its vector. The infection rate constant is derived ab initio by calculating collision frequency between the droplets/nuclei cloud and the susceptible population for different ambient conditions including the probability of infection. The SEIRD model output obtained with the calculated rate constants for average viral loading, for the specific conditions of interest, shows that preventing ejection of droplets with the initial diameter greater than 10 µm can potentially prevent further outbreaks even for a minimum infectious dose of 10 virions. The critical droplet diameter, preventing ejection of droplets above which would result in , is found to be 27 µm. For maximum viral loading, the critical droplet diameter is 5 µm to just prevent the outbreaks. Furthermore, the strong sensitivity of on the variation of virus half-life at different phases of the droplet/aerosol as well as the minimum infectious dose is demonstrated.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.