ABSTRACT
How human respiratory physiology and inhaled airflow therein proceed to impact transmission of SARS-CoV-2, leading to the initial infection, is an open question. An answer can help determine the susceptibility of an individual on exposure to a COVID-2019 carrier and can also quantify the still-unknown infectious dose for the disease. Combining computational fluid mechanics-based tracking of respiratory transport in anatomic domains with sputum assessment data from hospitalized COVID-19 patients and earlier measurements of ejecta size distribution during regular speech – this study shows that the regional deposition of virus-laden inhaled droplets at the initial nasopharyngeal infection sites, located in the upper airway, peaks over the droplet size range of 2.5 – 19 µ; and reveals that the number of virions that can potentially establish the infection is, at most, of 𝒪(102).
SEVERE acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has been identified as the causative agent for coronavirus disease 2019 (COVID-19), that has inflicted a global pandemic with over 17 million confirmed infections and nearly 700,000 deaths worldwide (as of late-July 2020).1
As is well-known by now, transmission of respiratory infections such as COVID-19 occurs through carriage of pathogens via droplets of different sizes produced during sneezing, coughing, singing, normal speech, and even, breathing.2 Accordingly, the means of person-to-person infection are projected to be three-way3: (a) inhalation of virus-laden droplets emitted by an infected individual at close-range; (b) inhalation of vaporized droplet nuclei that can float in air for hours; and (c) contaminating the respiratory mucosa through physical contact to external surfaces (fomites) with droplet deposits sitting on them. While (a) is valid for short-distance exposures to the COVID-19 carrier, transmission through modes (b) and (c) can happen over larger distances and longer time scales. However, clustering trends of infection spread (e.g. in industrial units4, in closed groups5, and inside households6) suggest that close-range exposures can be a critical determinant in worsening the pandemic. A follow-up question might be – what entails an exposure? A key component therein are the respiratory droplet sizes one is exposed to. Coughing and sneezing typically generate droplets with length-scales of 𝒪 (102) to 𝒪 (103) µ, while oral droplets ejected during normal speaking can range over ∼ 0.1 − 500 µ.3,7 The main competing effects determining the fate of these droplets are the ambient temperature and humidity (e.g. low relative humidity induces fast evaporation and shrinkage of the droplets), and the size of the droplet that controls its inertia and the gravitational force acting on it. While smaller droplets would stay airborne for longer, the larger droplets tend to fall fast ballistically; with the critical size for this transition being in the vicinity of 100 µ.8,9 Of note here, this study does not insist on any nomenclatural distinction between “aerosols” and “droplets” owing to ambiguities10 in common perception, and simply refers to all expiratory liquid particulates as droplets.
For mechanistic tracking of what range of virus-bearing droplet sizes might be more potent for transmission and to eventually induce infection, it is key that we identify the initial infection sites. A recent study11 used SARS-CoV-2 isolates to examine and compare the infection process of cultured cells from different parts of the human airway. The data reveals a striking pattern of continuous variation from a relatively high infectivity in ciliated epithelial cells along the nasal passage lining, to less infectivity in cells lining the throat and bronchia, and finally to relatively low infectivity at the lung cells. The pattern is governed through angiotensin-converting enzyme 2 (ACE2), which is a single-pass type I membrane protein and is the cell surface receptor that the virus utilizes to intrude into cells. ACE2 is abundant on ciliated epithelial cells, but is relatively scarce on the surface of the lower airway cells. While the findings are for in vitro samples, deposition of virus-laden droplets along the anterior nasal airway might not be so effective as to launch an infection despite the presence of ciliated cells, since the mucus layer provides some protection against virus invasion and infection.3 This leaves out the nasopharynx (the region along the respiratory tract posterior to the septum and comprising the upper part of the pharynx, see Panel A in Figure 1) as the main initial infection site; it acts as the seeding zone for subsequent infection of the lower airway via aspiration of virus-laden boluses of nasopharyngeal fluid. The ansatz is supported by data12 on the efficacy of nasopharyngeal swab testing for COVID-19 diagnosis, when compared to oropharyngeal swabs. So at this point, a valid question to ask would be: what are the typical droplet sizes that are making their way to the nasopharynx during inhalation?
A. The sagittal, axial, and coronal views of the CT-based nasal domains, shown respectively in a–c for anatomic reconstruction 1 (AR1) and d–f for anatomic reconstruction 2 (AR2). Nasopharynx is marked in red. B. Visuals of heat-maps for inspiratory transmission trends, showing the percentage of droplets of each size undergoing nasopharyngeal deposition (NPD). Data for different inhaled airflow rates are arranged along separate rows. Tracked droplet sizes are along the horizontal axis (positioned between the two heat-maps for AR1 and AR2). NPD peaks for droplets sized between 2.5 – 19 µ in AR1 and 2.5 – 15 µ in AR2. C. Correlation between RANS-based SST k-ω and LES results for the higher airflow rates i.e. 30, 55, and 85 L/min; therein the first three frames (bottom-left) are for AR1, the other three frames (bottom-right) correspond to data for AR2. The frames are on an aspect ratio of 0.5; ρ represents the Pearson’s correlation coefficient.
Respiratory droplets, on being expelled, undergo dehydration and consequent shrinkage, the extent of which partially depends on the fraction of non-volatile material in the droplets, such as dehydrated epithelial cell remnants, white blood cells, enzymes, DNA, sugars, electrolytes etc. So, although sputum is 99.5% water; ejected droplets, on dehydration, would have a higher density of 1.3 g/ml,13 which is what has been used for droplet tracking simulations here. This considers that the non-volatile weight fraction is in the 1 – 5% range. Such dehydration shrinks the expelled droplet diameter to 27 – 34% of the initial size. Therefore, assuming a mean 30% shrinkage and considering 100 µ as the critical size that prompts ballistic sedimentation, this study tracks inhaled droplet sizes in the range of 0.1 µ to 30% of 100 µ, i.e. 30 µ, for close-range exposure. The choice of the minimum in this size range, i.e. 0.1 µ, is dictated by SARS-CoV-2 size, which is in between 0.08 – 0.2 µ, with an average physical diameter of 0.1 µ.14
Next piece in this puzzle relates to the breathing parameters. Allometric relations15 approximate the minute inhalation at 18.20 L/min for a 75-kg male and 15.05 L/min for a 75-kg female, for gentle steady breathing while sitting awake. In general, the inspiratory flow rates can stretch over a wide range of ∼ 15 – 85 L/min, based on whether the individual is inhaling gently or breathing in forcefully. This study tests droplet transmission at four different inhalation rates, viz. 15, 30, 55, and 85 L/min; notably these discrete flow rates are the ones traditionally used in literature16 to gauge filtration capacities of protective mask respirators. Flow physics un-dergo a transition over this range; e.g. the 15 L/min flow through anatomically realistic nasal conduits lies in the laminar regime,17–28 the flow structures however devolve into turbulence29–32 at higher inhalation rates. Numerical schemes, that are used in this study to track respiratory transport, have been selected accordingly (see the supplementary methods section). Note that with nasopharynx being the initial infection zone, inhaled transport via mouth poses less likely for virion transmission and the modeling framework here focuses on the nasal domain in the upper airway.
Computed tomography (CT) scans of the nasal cavity from two healthy adult subjects were used in this study, to develop: anatomic reconstruction 1 (AR1) and 2 (AR2). Refer to the methods, for details of the in silico geometry building. Note that while it might be possible to access many research scans with some form of sinonasal irregularity, as they make part of disease-specific studies; it is relatively rare to have nasal CT scans from otherwise healthy subjects and hence limits the size of this study’s test cohort. For unbiased conclusions, it is, however, critical that we use disease-free CT-normal nasal airway domains to track the virion transport characteristics that would lead to the infection.
Results and discussion
Droplet size range that targets the nasopharynx
The overall droplet size range of 2.5 – 19 µ (in AR1: 2.5 – 19 µ, in AR2: 2.5 – 15 µ) registers the peak, in terms of the percentage of droplets of each size that are deposited at the nasopharynx. The range is determined by a cut-off of at least 5% deposition for around 3000 tracked droplets (viz. 3015 in AR1, 3000 in AR2) of each size. Panel B in Figure 1 displays the heat-maps for nasopharyngeal deposition (NPD) for different droplet sizes, during inhalation at the four tested airflow rates. The discrete droplet sizes, that were tracked, have been marked along the horizontal axis of the heat-maps. The patch bounded by the grey lines can, in fact, be a definitive graphical technique to delineate the hazardous droplet size range for various airborne transmissions. Note that these findings assume that the post-dehydration density of the respiratory droplets (expelled by the carrier and now being inhaled by the exposed individual) is at 1.3 g/ml. If there is little or no dehydration and as such the ejected droplet density remains at ∼ 1 g/ml, the inhaled droplet size range for peak NPD upscales to 3 – 20 µ (results in the online data repository33); since the slightly lighter droplets can now penetrate further into the intra-nasal airspace, the transport process being aided by the ambient inspiratory streamlines.
Also at this point, to think of a realistic exposure to a COVID-19 carrier, the vulnerable individual can be considered to inhale at different airflow rates over the duration of exposure. In such context, the Panel A in Figure 2 extracts the averaged nasopharyngeal deposition for the different tested inhalation rates in the two test subjects. Such inhalation-averaged transmission presents an approximate dehydrated droplet size range of 2.5 – 15.0 µ, for a minimum 2% NPD for each droplet size.
A. Graphical representation of the percentage of droplets of each size undergoing deposition at the nasopharynx; averaged from AR1 and AR2, for the tested inhalation rates. B. Distribution of droplet sizes ejected each minute during normal speaking, the numbers are calculated from earlier studies on expelled droplet tracking with food coloring.9 C. Pie diagram showing which droplet sizes are dominant contributors for virion transmission at the nasopharynx, for ejecta size distribution as in panel B. Symbol 𝔻 is the inhaled droplet diameter. The numbers assume that the droplets have undergone dehydration before being inhaled into the nasal airspace. D. Probabilistic interpretation of a droplet to contain at least 1 virion, based on whether the droplet size at inhalation is pre-dehydration or post-dehydration. E. Estimated number of virions that are deposited at the susceptible individual’s nasopharynx via dehydrated inhaled droplets, during close-range exposure to a COVID-19 carrier.
Droplets that are better at carrying the virions
The next pertinent question is: how effective are these droplets at carrying virions? SARS-CoV-2 belongs to a family of single-stranded RNA viruses, and virological assessments34 done on the sputum of hospitalized COVID-19 patients show an averaged viral load of 7 × 106 RNA copies/ml of oral fluid, with the peak load being 2.35 × 109 copies/ml. For the average load, simple calculations (see methods) show that the probability that a dehydrated 10-µ droplet (contracted from its original size of ∼ 33 µ) will carry at least 1 virion is 13.6%. The same number is 45.8% for a post-shrinkage 15-µ droplet. The probability drops exponentially to 0.2% for a 2.5-µ dehydrated droplet. Now, with existing data on the size distribution of expelled droplets during normal speaking (see Panel B, Figure 2), the proportion of virion deposits at the nasopharynx by different droplet sizes can be computed (see Panel C, Figure 2) using the transmission data presented in Figure 1. The deposition trends are again for droplets that are being inhaled post-dehydration.
As a note in this context, but albeit conspicuously enough, in the absence of dehydration, the probability of 1 virion being embedded in, for instance, a 10-µ droplet plummets to 0.37% (see Panel D, Figure 2). This rationalizes why in geographic regions with high humidity (and hence relatively less dehydration and shrinkage of respiratory ejecta), the pandemic’s spread has been somewhat measured.35,36
What could be COVID-19’s infectious dose?
The infectious dose is a fundamental virological measure quantifying the number of virions that can go on to start an infection; the value of which is still not conclusively known for SARS-CoV-2.37 Theoretically, according to the independent action hypothesis,38 even a single virion can potentially establish an infection in highly susceptible systems. Whether the hypothesis is true for humans and specifically for SARS-CoV-2 transmission is as yet undetermined. The rapid spread of COVID-19 though a priori suggests a small infective dose for the disease, that is ameliorating inter-human transmission.
Since it is unethical to expose subjects to SARS-CoV-2 (especially in the current absence of a remediating therapeutic – as of July 2020), this study introduces a strategy synergizing computational tracking and virological data, to quantify the infectious dose. Based on the nasopharyngeal transmission trends (Figure 1) and the virion transmission data (Panels B-C of Figure 2), for a 5-minute exposure: the number of virions depositing at the susceptible individual’s nasopharynx is 11 (considering average RNA load in the carrier’s sputum). If the infecting individual is in more serious stage of morbidity and hence considering the peak RNA load, as many as 3835 virions will be deposited on the nasopharynx of the exposed individual over 5 minutes (see Panel E in Figure 2).
Now to derive an estimate of infectious dose, consider the March 2020 Skagit Valley Chorale superspreading incident5 in the United States, where an infected person transmitted the disease to 52 other individuals in a 61-member choir group. Exposure time there was reported to be 2.5 hours; consequently, for an average RNA load (assuming that the carrier had mild-to-moderate symptoms), the number of virions depositing at a closely-positioned susceptible individual’s nasopharynx over that duration approximates to (11/5) × 2.5 × 60 ≈ 330. So, ∼ 300 can be considered a conservative upper estimate for the COVID-19 infective dose, the order agreeing with preliminary estimates from replication rates of the virus.39
That the number of virions needed to establish the infection is, at most, of 𝒪 (102) is indeed remarkable! The scale clearly manifests how communicable this disease is, especially if discerned in the perspective of infectious doses for other airborne transmissions, e.g. the infective dose for influenza A virus, when administered through aerosols to human subjects lacking serum neutralizing antibodies, ranges between 1950 – 3000 virions.40 The discovery of infectious dose, together with this study’s detection of droplet sizes that target the infection-prone nasopharynx, could provide a useful resource in mitigating the COVID-19 pandemic and in the mechanistic design of topical anti-viral therapeutics41–43 and targeted intra-nasal vaccines44,45 that can induce prophylaxis through mucosal and systemic immune responses.
Methods
Anatomic geometry reconstruction
Medical imaging-based in silico model generation was accomplished according to relevant guidelines and regulations, with the anatomic geometries being reconstructed from existing de-identified imaging data from two CT-normal subjects. The use was approved with exempt status by the Institutional Review Board at UNC Chapel Hill. The test subjects include a 61 year-old female (AR1) and a 37 year-old female (AR2). The nasal airspaces were extracted from the medical grade scans (with the CT slices collected at coronal depth increments of ∼ 0.4 mm) over a delineation range of -1024 to -300 Hounsfield units, and was complemented by careful hand-editing of the selected pixels to ensure anatomic accuracy. For this step, the DICOM (Digital Imaging and Communications in Medicine) scans for each subject were imported to the image processing software Mimics 18.0 (Materialise, Plymouth, Michigan).
The reconstructed geometries were imported as stere-olithography files to ICEM-CFD 15.0 (ANSYS, Inc., Canons-burg, Pennsylvania), and then meshed spatially into minute volume elements. As per mesh refinement-based protocols,46,47 each computational grid contained more than 4 million un-structured, graded tetrahedral elements (e.g. 4.54 million in AR1, 4.89 million in AR2); along with three prism layers of 0.1-mm thickness at the airway walls, with a height ratio of 1. The nostril inlet planes comprised 3015 elements in AR1 (1395 elements on left nostril plane, 1620 elements on right nostril plane) and 3000 elements in AR2 (1605 on left nostril plane, 1395 on right nostril plane).
Numerical simulations
The study considers droplet transport for four different inhaled airflow rates, viz. 15, 30, 55, and 85 L/min. The lower flow rate (i.e. 15 L/min) corresponds to comfortable resting breathing, with the viscous-laminar steady-state flow physics model standing in as a close approximation.17–28 At higher flow rates (extreme values of which may sometimes lead to nasal valve collapse), the shear layer separation from the tortuous walls of the anatomic geometries results in turbulence.29–32 While accounting for the turbulent characteristics of the ambient airflow, the study averages the droplet deposition percentages from implementation of two distinct categories of numerical schemes, viz. (a) shear stress transport (SST) based k-ω model, which is a sub-class under Reynolds-averaged Navier Stokes (RANS) schemes that parameterize the action of all turbulent fluctuations on to the mean flow; and (b) Large Eddy Simulation (LES). The two numerical techniques depict high correlation in terms of droplet deposition at the nasopharynx (see Figure 1, Panel C). However, it should be noted that while the SST k-ω scheme, a 2-equation eddy-viscosity model, is computationally less expensive; it averages the short time-scale flow artifacts, such as the transient vortices, that might emerge in these complicated domains; and hence the prediction of droplet transport affected by the simulated ambient airflow may at times contain errors (e.g. with the low-pressure vortex centers acting as droplet attractors). LES is computationally more expensive, it separates the turbulent flow into large-scale and small-scale motions, and accounts for the small fluctuations through a sub-grid scale model (in this study, Kinetic Energy Transport Model was used as the sub-grid scale model48). We take the averaged estimates of NPD from the two schemes, to minimize statistical and algorithmic biases.
The computational schemes employed a segregated solver on ANSYS Fluent, with SIMPLEC pressure-velocity coupling and second-order upwind spatial discretization. Solution convergence was monitored by minimizing the mass continuity and velocity component residuals, and through stabilizing the mass flow rate and static pressure at the airflow outlets. For the pressure-driven flow solutions: typical convergence run-time for a laminar simulation with 5000 iterations was approximately 5–6 hours for 4-processor based parallel computations executed at 4.0 GHz speed. The corresponding run-time for a RANS simulation was ∼ 12 hours; for an LES computation, it was 4–5 days. Note that for the LES work, the simulated flow time was 0.5 sec for the 30 L/min case, with 0.0002 sec as the time-step49 and it was 0.25 sec for the 55 and 85 L/min flow rates with the time-step at 0.0001 sec. ρ = 1.204 kg/m3 was used as air-density, and µ = 1.825 × 10−5 kg/m.s was assumed as air’s dynamic viscosity.
Following set of boundary conditions were enforced during the simulations: (i) zero velocity at the airway-tissue interface i.e. at the walls enclosing the nasal airspace (other-wise commonly referred to as the no slip condition), along with “trap” boundary condition for droplets whereby a droplet would come to rest after depositing on the walls; (ii) zero pressure at nostril planes, which were the pressure-inlet zones in the simulations, with “reflect” boundary condition for droplets to mimic the effect of inhalation on the droplet trajectories if they are about to fall out of the anterior nasal domain; and (iii) a negative pressure at the airflow outlet plane, which was the pressure-outlet zone, with “escape” boundary condition for droplets, i.e. allowing for the outgoing droplet trajectories to leave the upper respiratory airspace. Mean inlet-to-outlet pressure gradients were -9.01 Pa at 15 L/min, -26.65 Pa at 30 L/min, -73.73 Pa at 55 L/min, and -155.93 Pa at 85 L/min. For a reference on the general layout of the anatomic regions, see Panel A in Figure 1.
On convergence of the airflow simulations, inhaled droplet dynamics were tracked by Lagrangian-based discrete phase inert particle transport simulations in the ambient airflow; with the localized deposition along the airway walls obtained through numerically integrating transport equations that consider contribution of the airflow field on the evolution of droplet trajectories, along with the effects for gravity and other body forces such as the Saffman lift force that is exerted by a flow-shear field on small particulates moving transverse to the streamwise direction. Also, the droplet size range is considered large enough to discount Brownian motion effects on their spatial dynamics. Note that the study simulated the transport for 3015 droplets of each size in AR1 and 3000 droplets of each size in AR2, the numbers being same as the number of elements on the nostril inlet planes which were seeded with the to-be-tracked droplets for the droplet transport simulations. For the numerical tracking, the initial mass flow rate of the inert droplets moving normal to the inlet planes into the nasal airspace was required to be non-zero, and was set at 10−20 kg/s. After the transport simulations, the post-processing of the droplet deposition data along the airway walls provides the NPD. This numerical protocol has been rigorously validated50 by comparing the deposition trends along the inner walls of similar in silico digital models to the in vitro spray tests performed in 3D-printed solid replicas of the same reconstructions.
Statistical methods and data interpretation
Panel C in Figure 1 plots the NPD values from RANS (along horizontal axis) and LES (along vertical axis) schemes, implemented for the higher inhalation rates (i.e. 30, 55, and 85 L/min). The simulation outputs are linearly correlated with an average Pearson’s correlation coefficient of 0.98 for 30 L/min, 0.93 for 55 L/min, and 0.95 for 85 L/min. Subsequent check of the slope m for the linear best-fit trendline, through the scatter plots of RANS and LES-based NPD data, indicates how similar the estimates are quantitatively; the mean measures therein being m = 1.113 for 30 L/min, m = 1.052 for 55 L/min, and m = 1.177 for 85 L/min; with the value 1 signifying exact equivalence. The statistical operations were carried out on Wolfram Mathematica.
For the ejecta size distribution (Panel B, Figure 2), the study divides up the reported7 percentages for each size bin (i.e. 0–5, 5–10, 10–15 µ etc.) uniformly and apportions them to the discrete droplet sizes (belonging to the same size bin) that are tracked (see horizontal axis between the NPD heat-maps in Panel B, Figure 1), to estimate how many droplets of each size would be ejected by the carrier during unit time. The referenced article7 described the size bin limits as A − B; for consistency, this study interpreted that as droplet sizes (in µ) that are ≥ A and < B. Also, note the use of parentheses and square brackets to define the size bins on Panel B in Figure 2; e.g. [α, β) implies, as per set theory notations, the range of droplet sizes (in µ) that are ≥ α and < β.
Probabilistic estimation of virion contamination in expiratory droplets
Suppose the viral load in a COVID-19 carrier has been assessed to be V copies of RNA in each ml of sputum fluid. Let a typical expelled droplet diameter from the carrier be 𝔻 µ. With SARS-CoV-2 being a single-stranded RNA virus, the average number of virions embedded in each droplet can then be computed as (π/6) V 𝔻3 ×10−12. Therefore, every 100 droplets of the same size would have (π/6) V 𝔻3 ×10−10 virions; which, in other words, also represents the probability (in %) for a droplet of diameter 𝔻 µ, of containing at least 1 virion.
The above mathematical approach, however, by the nature of it, provides a simplistic estimate of viral load in the ejected droplets, by assuming the spatial distribution of virions to be uniform in the sputum. In reality, these considerations also lead to a vigorous area of research – on how the rheology of oral fluids might affect the ejecta generation and break-down, and the resultant volumetric concentration of virions in the expiratory remnants.
Data availability
This project has generated simulated quantitative, de-identified data on regional deposition over nasal tissues. The data-sets (including Fluent .cas and .dat files) and the numeric protocols; along with MATLAB codes, Wolfram Mathematica notebooks, and Microsoft Excel spreadsheets used for data post-processing – are available from the corresponding author on request (by contacting via the correspondence email addresses), through a shared-domain Google Drive folder.33
Author contributions
S.B. conceived the study design, performed research, analyzed the data, and drafted the manuscript.
Competing interests
The author declares no competing financial interests.
Additional information
A version of this manuscript has been screened for content and posted on medRxiv. doi: https://doi.org/10.1101/2020.07.27.20162362
Acknowledgements
This material is based upon work supported by a National Science Foundation (NSF) RAPID grant (Award Number 2028069, with S.B. as the Principal Investigator) for COVID-19 research. Any opinions, findings, and conclusions or recommendations expressed here are, however, those of the author and do not necessarily reflect NSF’s views. Supplemental assistance for the project came from S.B.’s faculty start-up package at South Dakota State University.
The author acknowledges A.J. Kimple, MD, PhD (attending rhinologist at the Department of Otolaryngology / Head and Neck Surgery, UNC Chapel Hill); L.P. Chamorro, PhD (faculty at the Department of Mechanical Science and Engineering, University of Illinois at Urbana-Champaign); and S. Jung, PhD (faculty at the Department of Biological and Environmental Engineering, Cornell University) for several fruitful discussions. Acknowledgements are also due to J.S. Kimbell, PhD (research faculty at the Department of Otolaryngology / Head and Neck Surgery, UNC Chapel Hill) for mentoring and support.
Computing facilities at both SDSU and UNC Chapel Hill were used for the work.
Footnotes
Additional details have been added to the Methods section, to clarify the statistical techniques and strategies on data interpretation. The reference list has been revised and updated as well. Other major edits are in the sub-section for probabilistic estimation of virion contamination of expiratory droplets.
References
