In Silico Quantification of Intersubject Variability on Aerosol Deposition in the Oral Airway

The extrathoracic oral airway is not only a major mechanical barrier for pharmaceutical aerosols to reach the lung but also a major source of variability in lung deposition. Using computational fluid dynamics, deposition of 1–30 µm particles was predicted in 11 CT-based models of the oral airways of adults. Simulations were performed for mouth breathing during both inspiration and expiration at two steady-state flow rates representative of resting/nebulizer use (18 L/min) and of dry powder inhaler (DPI) use (45 L/min). Consistent with previous in vitro studies, there was a large intersubject variability in oral deposition. For an optimal size distribution of 1–5 µm for pharmaceutical aerosols, our data suggest that >75% of the inhaled aerosol is delivered to the intrathoracic lungs in most subjects when using a nebulizer but only in about half the subjects when using a DPI. There was no significant difference in oral deposition efficiency between inspiration and expiration, unlike subregional deposition, which shows significantly different patterns between the two breathing phases. These results highlight the need for incorporating a morphological variation of the upper airway in predictive models of aerosol deposition for accurate predictions of particle dosimetry in the intrathoracic region of the lung.


Introduction
The extrathoracic upper airway acts as the first line of defense to prevent inhaled toxicants from reaching the lungs but also as a barrier to the delivery of inhaled drugs. The proportion of aerosol deposited in the oral cavity and throat depends on the flow field and the size of inhaled particles [1]. The intricate anatomy of the upper airway produces complex flow patterns and particle trajectories and is thus also an important factor affecting aerosol deposition in the human respiratory system [2]. Several studies have shown that deposition in the throat is a major source of variability in lung deposition [3][4][5].
Aerosol deposition in the oral airway has been studied both computationally [2,6] and experimentally [7][8][9]. Previous in vitro work with steady flow rates presented empirical correlations for predicting oral deposition [10,11]. However, previous studies mostly evaluated total upper airway deposition and have not distinguished deposition in the oropharyngeal and laryngeal regions. This distinction could be important. For instance, laryngeal deposition of inhaled corticosteroids (ICS), a mainstay in the treatment of chronic reactive airway disease, elicits local side effects, including dysphonia [12,13], that could be minimized with optimized flow rate for specific upper airway anatomy.
Understanding the mechanics of particle deposition in the upper airway is useful to design new inhalation therapies for respiratory diseases [14]. It also helps estimate exposure risks to inhaled toxicants [15] or airborne transmission of SARS-CoV-2-laden droplets [16].

Reconstruction of Human Upper Airway Geometries
CT scans in DICOM format were used to create three-dimensional (3D) realistic models of the oral airway anatomy using Mimics 23.0 (Materialise Inc., Leuven, Belgium). Briefly, by setting a previously identified optimal threshold between −1024 and −300 Hounsfield Units (HU) [21,22], the 3D anatomy of the oral cavity, including the mouth, oropharyngeal region, larynx, vocal cord, and upper tracheal sections, was selected as the region of interest (ROI). The nasal cavity and paranasal sinuses were excluded. The oral airway geometries were exported in stereolithography (stl) format to ICEM-CFD 21.0 (ANSYS 2021 R2, Inc., Canonsburg, PA, USA), where the planar inlet and outlet surfaces were defined. All geometries were oriented with the oral cavity floor parallel to the Z-axis ( Figure 1). Planar divisions between the posterior mouth, oropharynx, area of larynx and trachea were Pharmaceutics 2023, 15, 160 3 of 18 created for recording region-specific particle deposition in airflow simulations. Consistent with experimental conditions, a 245 mm-long tube was added to the mouthpiece. oropharyngeal region, larynx, vocal cord, and upper tracheal sections, was selected as the region of interest (ROI). The nasal cavity and paranasal sinuses were excluded. The oral airway geometries were exported in stereolithography (stl) format to ICEM-CFD 21.0 (ANSYS 2021 R2, Inc, Canonsburg, Pennsylvania), where the planar inlet and outlet surfaces were defined. All geometries were oriented with the oral cavity floor parallel to the Z-axis (Figure 1). Planar divisions between the posterior mouth, oropharynx, area of larynx and trachea were created for recording region-specific particle deposition in airflow simulations. Consistent with experimental conditions, a 245 mm-long tube was added to the mouthpiece.

Flow Simulations
CFPD simulations were performed by solving Navier-Stokes equations with Lagrangian particle tracking. The volume fraction of the particles in the air was in the range of 10 −7 to 10 −9 . Particle number density was thus sufficiently low so that particle motion did not affect airflow and a one-way coupling between airflow and particle transport was used [23]. The Navier-Stokes equation for incompressible flow is where = ( , , , ) is the fluid velocity vector, = 1.204 kg/m 3 is the fluid density, = 1.825 × 10 −5 kg/ (m.s) is the dynamic viscosity of the Newtonian fluid (i.e., air), is pressure and is time. To solve the Navier-Stokes equation, mesh independence analysis has been performed, and geometries were meshed in ICEM-CFD 21.0, a preprocessor tool from Ansys Fluent, with approximately twelve million tetrahedral elements containing eight prism layers with a total prism zone thickness of 0.2 mm. Prism mesh was generated   simulations were performed by solving Navier-Stokes equations with Lagrangian particle tracking. The volume fraction of the particles in the air was in the range of 10 −7 to 10 −9 . Particle number density was thus sufficiently low so that particle motion did not affect airflow and a one-way coupling between airflow and particle transport was used [23]. The Navier-Stokes equation for incompressible flow is where u = u(x, y, z, t) is the fluid velocity vector, ρ = 1.204 kg/m 3 is the fluid density, µ = 1.825 × 10 −5 kg/ (m.s) is the dynamic viscosity of the Newtonian fluid (i.e., air), p is pressure and t is time. To solve the Navier-Stokes equation, mesh independence analysis has been performed, and geometries were meshed in ICEM-CFD 21.0, a preprocessor tool from Ansys Fluent, with approximately twelve million tetrahedral elements containing eight prism layers with a total prism zone thickness of 0.2 mm. Prism mesh was generated near the wall to ensure capturing of the steep air velocity gradient and of accurate particle deposition patterns in the simulations. Inhalation rates representative of resting and fast breathing were considered. Mesh quality was checked by considering the quality criteria of greater than 0.1 to ensure that the mesh had no distorted elements [24]. Experimental evidence suggests that flow in the extrathoracic airway is laminar at resting breathing rates [25,26]. Thus, for slow-breathing condition (18 L/min), Equation (1) was solved in the steady-state laminar flow regime by employing Ansys Fluent 2021 R2 (ANSYS, Inc., Canonsburg, PA, USA) and choosing SIMPLEC algorithm for pressurevelocity coupling and second-order upwind discretization method. Steady-state flow simulations were conducted by imposing boundary conditions on the inlet mass flow rate of the air and pressure outlet boundary conditions applied at the outlet boundary.
To ensure numerically accurate flow simulation results, the Shear-Stress Transport k-ω turbulent model (SST k-ω) was employed to simulate the inspiratory and expiratory airflow for fast-breathing conditions (45 L/min) [27]. The turbulence length scale was considered 0.001 m, and turbulent intensity was assumed to be 5% at the inlet and outlet [28,29].
The following boundary conditions for the flow were applied: (1) zero velocity at the walls (no-slip boundary condition), (2) inlet mass flow at the mouthpiece's inlet to set airflow of 18 and 45 L/min, and (3) zero pressure at the outlet.

Particle Transport Simulations
Forces affecting particles for the size range used in this study included drag force and gravity (consistent with experimental orientation), allowing modeling of deposition mechanisms of inertial impaction and sedimentation. The equation of motion governing the trajectory of a particle is where u p is the particle velocity, g is the gravity field, ρ p is the density of the particle, F D u − u p is the drag force per unit particle mass and where d p is the particle diameter, Re p , the particle Reynolds number and C D is the drag coefficient [24]. The Lagrangian-based Discrete Phase Model (DPM) was used to predict particle deposition in the eleven anatomically realistic models. Deposition of particles with mass median aerodynamic diameters (MMAD) of 1-30 µm (1 µm increments) was investigated during resting and fast-breathing conditions. Particles were considered inert with a spherical shape and density of 1000 kg/m 3 so that the particle diameter corresponds to the aerodynamic diameter. Particles were released from the inlet of the tube. As the tube and mouthpiece were not in the region of interest, the boundary condition "reflect" was applied in the DPM simulation. Boundary condition "trap" was considered for the mouth cavity, oropharyngeal, laryngeal, and tracheal region to predict aerosol deposition. The "escape" boundary condition was used at the outlet. The trajectories of 10,000 particles were simulated for each particle size and inhalation flow rate, allowing for particle number-independent predictions. Particles were injected at the inlet with a blunt profile and random spatial distribution. Increasing to 100,000 particles changed predicted deposition by less than 0.38%.
The percentages of inhaled particles deposited in each anatomical region ( Figure 1) were quantified. For instance, particle deposition in the laryngeal region was defined as 100 × (N L /N I ), where N I is the total number of particles inhaled and N L is the number of particles deposited on the surface region mapped as the laryngeal region on the airway model. Similarly, the total deposition efficiency was defined as 100 × (N T /N I ), where N T is the total number of particles deposited in all airway regions combined.

Comparison of Subject-Specific In Silico Predictions with Experimental Results
Whole-lung deposition predictions were obtained for breathing parameters and functional residual capacity (FRC) that matched on a subject-by-subject basis those measured in 7 healthy subjects during controlled breathing of aerosols [20]. In these experiments, subjects were asked to target a tidal volume of 1000 mL of particle-laden air (1 and 2.9 µm) at constant inhalation flow rates of 18 and 45 L/min. Actual tidal volumes and flow rates measured over five consecutive breaths of controlled exposure for each experimental condition are listed in Table 2 along with each subject FRC. These data were used in the in silico predictions used in the comparison. The overall retained fraction was calculated as the sum of the deposition fraction in the oral cavity and retained fraction in the intrathoracic lung. The intrathoracic deposition was estimated with the latest version of the Multiple-Path Particle Dosimetry (MPPD) model that includes a mechanistically based model component for alveolar mixing of particles and that accounts for multiple breaths of aerosol intake [30]. Deposition in the oral cavity was obtained from CFD simulations as described above. Deposition in the oral cavity was also obtained from the semi-empirical equation of Stahlhofen et al. [31]: where d a is the aerodynamic diameter expressed in µm and Q is the overall volumetric flow rate expressed in mL/s. The mass of particles injected in the MPPD model was set as (1 − η oral ) × C inh where C inh is the inhaled particle concentration. Deposition in the oral cavity during exhalation was based on the mass of particles exiting MPPD (η oral C exh.MPPD ).

Statistical Analysis
The curve of best fit of in silico predictions of aerosol deposition in the oral cavity was calculated by fitting β 1 and β 2 to a sigmoidal function where d 2 a Q is the same impaction parameter used in the Stahlhofen equation (Equation (4)). β 1 , β 2 , and the percentage of variance explained (R 2 ) by the curve of best fit were compared to the parameters adapted in the Stahlhofen equation.
To compare in silico predictions with experimental data, a one-way ANOVA test for correlated samples was used with the Tukey Multiple Comparison post hoc test. The paired t test was used to compare in silico predictions of oral deposition between inspiration and expiration. Significant differences were accepted at the p < 0.05 level.
MATLAB R2022a was used for curve fitting and confidence interval estimation, and R v4.0.2 was used for statistical analysis and visualization. Figure 2A illustrates oral deposition efficiency as a function of the commonly used impaction parameter d a 2 Q, where d a (µm) is the aerodynamic diameter and Q (L/min) is the flow rate. In agreement with previous in vitro and in silico studies [8,9,18,32], the high correlation ( Figure 2A, best fit, R 2 = 86.02%) between oral deposition and the impaction parameter d a 2 Q shows that inertial impaction is the dominant deposition mechanism in the upper airway. Total oral deposition of inhaled particles increased with increasing particle size and inhalation flow rate. For instance, in subject H1 and for a flow rate of 18 L/min, oral deposition increases from 0.26% for 1 µm particles (gray × symbol at d a 2 Q = 18 µm 2 L/min, Figure 2A) to 1.49% for 5µm particles (gray × symbol at d a 2 Q = 450 µm 2 L/min) and to 20.71% for 10 µm particles (gray × symbol at d a 2 Q = 1800 µm 2 L/min). In other words, 99.74% of 1 µm particles, 98.51% of 5µm particles and 79.29% of 10 µm particles traversed the airway model and were delivered to the intrathoracic region of the lung. Increasing the inhalation rate caused a substantial increase in the percentage of particles deposited in the upper airway. Simulations for a 45 L/min inhalation rate predicted 1.21% of 1 µm particles (black × symbol at d a 2 Q = 45 µm 2 L/min, Figure 2A) and 76.78% of 10 µm particles (black × symbol at d a 2 Q = 4500 µm 2 L/min) depositing in the upper airway model, respectively. While most of 1 µm particles were delivered to the trachea at both flow rates, <25% of inhaled 10 µm particles were delivered to the intrathoracic lungs at the higher flow rate compared to~80% at the low flow rate. Despite larger intrasubject variability in oral deposition, similar trends were observed for all subjects ( Figure 2). For instance, deposition of 3 µm aerosols ranged from 0.51% to 11.59% (median = 0.97%) between subjects at 18 L/min ( Figure 2B) and from 0.84% to 52.32% (median 6.24%) at 45 L/min ( Figure 2C). There was no significant difference in oral airway deposition between healthy and mild-to-moderate COPD subjects.

Total Oral Deposition Efficiency
Predictions of oral deposition in both healthy and COPD cohorts show a good agreement with the empirical curve previously obtained by Stahlhofen et al. [31] from controlled in vivo experiments (red curve, Figure 2A). The equation of best fit to patient-specific CFD data was [1 best fit = 0.8602).

Effect of Particle Impaction on Regional Deposition
Regional deposition fractions were computed as a function of d a 2 Q for all the subregions of the upper airway as defined in Figure 1. There was significant intersubject variability in the distribution of deposited particles among all three subregions of the upper airway ( Figure 3A-C). Deposition in the mouth cavity increased with increasing particle size and increasing flow rate ( Figure 3A). Deposition in both the oropharyngeal and the laryngeal region followed a bell shape, with the maximum deposition varying largely between subjects: maximum deposition occurred for an impaction parameter ranging between 10 3 and 1.2 × 10 4 in the oropharyngeal region ( Figure 3B) and between 400 and 5000 in the laryngeal region ( Figure 3C).

Comparison of Inspiratory and Expiratory Particle Deposition in the Upper Airway
The effect of flow direction (inspiratory versus expiratory flow) on oral deposition was also investigated. For expiratory flow simulations, particles were injected at the tracheal outlet with a blunt profile and random spatial distribution. Figure 4 shows the spatial distribution of deposited particles in two subjects with highly different upper airway shapes following inhalation or exhalation of 3 µm particles. In subject H5, hotspots of deposited particles were mainly found on the posterior laryngeal wall and at the level of the vocal cords following inhalation ( Figure 4A) and on the anterior laryngeal wall following exhalation ( Figure 4C). In contrast, in subject H6, most deposited particles were located on the posterior oropharyngeal wall following inhalation ( Figure 4B) and on the anterior laryngeal wall and on the hard palate (mouth cavity) following exhalation ( Figure 4D). The hotspots were associated with regions of high airflow velocities ( Figure A2), again suggesting that inertial impaction is the dominant deposition mechanism in the upper airway. These data not only highlight the large variability in deposition patterns between subjects but also between breathing phases, i.e., between inspiration and expiration. airway ( Figure 3A-C). Deposition in the mouth cavity increased with increasing particle size and increasing flow rate ( Figure 3A). Deposition in both the oropharyngeal and the laryngeal region followed a bell shape, with the maximum deposition varying largely between subjects: maximum deposition occurred for an impaction parameter ranging between 10 3 and 1.2 × 10 4 in the oropharyngeal region ( Figure 3B) and between 400 and 5000 in the laryngeal region ( Figure 3C).

Comparison of Inspiratory and Expiratory Particle Deposition in the Upper Airway
The effect of flow direction (inspiratory versus expiratory flow) on oral deposition was also investigated. For expiratory flow simulations, particles were injected at the tracheal outlet with a blunt profile and random spatial distribution. Figure 4 shows the Distribution of 1, 3 and 5 µm deposited particles at 18 L/min and 45 L/min breathing conditions among subregions of the subject-specific oral airway is shown for all subjects in Figure 5 both for inhalation and exhalation. These CFD results do not show any consistent trend when deposition occurring during inhalation is compared to that during expiration, with some subjects showing higher deposition during inhalation, others showing higher deposition during exhalation and a third group showing similar values between inspiration and expiration. As a result, there was no significant difference in oral deposition between inspiration and expiration in this group of subjects. This was also true for subregion deposition except for particles ≥ 5 µm ( Figure 5F) where most particles deposited in the laryngeal region during expiration, leaving few particles to travel and potentially deposit in the oropharynx and mouth cavity. Regional deposition for all particle sizes (1-30 µm) and expiratory flow rate (18 and 45 L/min) is shown in Figure 6. Deposition in the laryngeal region increased with increasing particle size and increasing expiratory flow rate ( Figure 6C), while deposition in both the oropharyngeal region and the mouth cavity followed a bell shape, with the maximum deposition varying largely between subjects ( Figure 6A,B).
Pharmaceutics 2023, 15, 160 9 of 18 4D). The hotspots were associated with regions of high airflow velocities ( Figure B1), again suggesting that inertial impaction is the dominant deposition mechanism in the upper airway. These data not only highlight the large variability in deposition patterns between subjects but also between breathing phases, i.e., between inspiration and expiration.  showing higher deposition during exhalation and a third group showing similar values between inspiration and expiration. As a result, there was no significant difference in oral deposition between inspiration and expiration in this group of subjects. This was also true for subregion deposition except for particles ≥5 μm ( Figure 5F) where most particles deposited in the laryngeal region during expiration, leaving few particles to travel and potentially deposit in the oropharynx and mouth cavity. Regional deposition for all particle sizes (1-30 µ m) and expiratory flow rate (18 and 45 L/min) is shown in Figure 6. Deposition in the laryngeal region increased with increasing particle size and increasing expiratory flow rate ( Figure 6C), while deposition in both the oropharyngeal region and the mouth cavity followed a bell shape, with the maximum deposition varying largely between subjects ( Figure 6A,B).

Comparison of In Silico Predictions with In Vivo Experimental Data
Whole-lung deposition predictions were obtained by coupling subject-specific CFD results with MPPD predictions as described in Section 2.4. Figure 7A displays these predictions against experimental data obtained by Darquenne et al. [20]. Experimental data are displayed as deposition measured over five consecutive breaths (mean ± standard deviation (SD)). Whole-lung deposition was also obtained by coupling predictions from Equation (4) with MPPD ( Figure 7B). The regression line between in silico predictions (y) and experimental data (x) were y = 0.931x + 0.030; R 2 = 61%, for MPPD/CFD results ( Figure  7A, dashed line) and y = 0.995x − 0.028; R 2 = 73%, for MPPD/empirical results ( Figure 7B). Figure 8 compares whole-lung deposition predictions for the seven subjects grouped

Comparison of In Silico Predictions with In Vivo Experimental Data
Whole-lung deposition predictions were obtained by coupling subject-specific CFD results with MPPD predictions as described in Section 2.4. Figure 7A displays these predictions against experimental data obtained by Darquenne et al. [20]. Experimental data are displayed as deposition measured over five consecutive breaths (mean ± standard deviation (SD)). Whole-lung deposition was also obtained by coupling predictions from Equation (4) with MPPD ( Figure 7B). The regression line between in silico predictions (y) and experimental data (x) were y = 0.931x + 0.030; R 2 = 61%, for MPPD/CFD results ( Figure 7A, dashed line) and y = 0.995x − 0.028; R 2 = 73%, for MPPD/empirical results ( Figure 7B). Figure 8 compares whole-lung deposition predictions for the seven subjects grouped together for both the MPPD/CFD and MPPD/empirical cases and the experimental data. Data are presented as median (minimum, maximum).

Intersubject Variability in Deposition in the Oral Airway
Oral deposition of particles in the aerodynamic size range of 1-30 µm was numerically predicted in distinct geometries of oral airways of seven healthy adults and four subjects with mild-to-moderate COPD. Lung diseases such as COPD significantly alter the deposition of inhaled particles in the intrathoracic lungs [33]. This is mainly because of alterations in the geometry of the airways and alveolar spaces that result in regional changes in the ventilation distribution of inhaled air and flow patterns in the airspaces [34]. To our knowledge, there is no report of significant alterations in the geometry of the extrathoracic airway between healthy subjects and those with mild-to-moderate COPD. For similar exposure conditions, we found no significant difference in oral airway deposition between healthy and mild-to-moderate COPD subjects.
For any given combination of particle size and inhaled flow rate, a large scatter was observed between subjects. For example, deposition ranged from 0-18%, 1-56% and 2-85% at an impact parameter of~200, 400 and 1000 µm 2 L/min, respectively (Figure 2A). These data compare well with in vitro measurements obtained in nine replicas of oral airways where deposition ranged from 0-30%, 0-60% and 5-95% at an impact parameter of 200, 400 and 1000 µ 2 L/min, respectively [9]. A similar scatter in oral deposition between subjects was also observed by Grgic et al. in a separate in vitro study [8].
Data from in vivo studies have also reported large intersubject variability [35][36][37][38][39][40][41]. Using these in vivo data, Stahlhofen et al. derived a semi-empirical equation based on particle size and breathing pattern characterized by the impaction parameter (Equation (4)). We developed a similar correlation based on our in silico predictions (Figure 2A). Despite pronounced intersubject variability in deposition predictions, the sigmoidal curve of in silico oral deposition versus impaction parameter remains statistically indifferent from the experimental curve predicted by the Stahlhofen equation, which alone explains the experimental data well with a best fit R 2 value of 85.98%. This compares to an R 2 of 86.02% for the sigmoidal curve based on in silico predictions.
The flow rates used in this study are representative of those achieved during aerosol delivered by nebulizers (18 L/min, tidal breathing) and dry powder inhalers (DPIs, 45 L/min, target 30-60 L/min). With an optimal size distribution of 1-5 µm for pharmaceutical aerosols, our data suggest that, with a nebulizer, more than 75% of the inhaled aerosol would be delivered to the intrathoracic lungs in most subjects (gray area, Figure 2B), with only a few of them having >25% of the larger particles (i.e., 5 µm) depositing in the upper airway. In contrast, there is a much larger variability with the use of a DPI with only about half of the subjects delivering >75% of the 1-5 µm aerosol to the intrathoracic lungs (gray area, Figure 2C). These data show that a more consistent aerosol dose among subjects can be delivered to the intrathoracic lungs when using nebulizers (that require a lower flow rate) than DPIs (that require a higher flow rate). In terms of in silico population studies of aerosol deposition in the upper airway and tracheobronchial trees, these data also suggest that, for aerosol exposure conditions that resulted in low inter-subject variability, a standardized oral upper airway model could be used, resulting in significant computational time savings [42].
Fewer studies have looked at regional deposition within the oral extrathoracic airway [8,43] but it is often assumed that, during inhalation, most deposition occurs at the level of sudden constrictions present in the oropharyngeal and laryngeal regions [31,38]. Our data suggest this to be the case for impaction parameters up to~1000-2000 µm 2 L/min. However, for larger values of the impaction parameter, deposition in the mouth also becomes significant ( Figure 3A). As a result, deposition in the oropharynx resembles a bell-shaped curve ( Figure 3B) rather than the sigmoid shape seen in the mouth. A similar bell-shaped curve was obtained in the larynx ( Figure 3C). Compared to deposition in the oropharynx, the curve was however further shifted towards the smaller particles due to additional upstream deposition. Finally, a comparison of our in silico predictions with in vitro data obtained in seven realistic mouth-throat geometries for impaction parameters ranging from 270 to 3800 µm 2 L/min shows relatively good agreement [8]. Indeed, these in vitro data not only demonstrated high intersubject variability in the regional deposition but also significant deposition in the mouth for experiments performed with impaction parameters around 1000 µm 2 L/min or larger.

Differences in Regional Deposition between Inhalation and Exhalation
Most of the studies on upper airway deposition during mouth breathing have focused on the inspiratory phase, e.g [2,5,8,9,32], with very few looking at a deposition during expiration [7,43]. This is mainly because inhalation drug therapies are designed to maximize deposition in the lungs prior to exhalation. Indeed, when using a dry powder inhaler (DPI) or a pressurized metered-dose inhaler (pMDI), proper drug inhalation techniques call for a deep inhalation followed by a breath hold to allow for particles to settle in the airspaces, minimizing the number of particles being exhaled [44]. However, poor inhalation technique can result in a significant fraction of exhaled aerosol. Additionally, drugs administered with nebulizers do not typically include an end-inspiratory breath-hold, which can result in significant exhaled aerosol fractions.
Due to the paucity of available data for the exhalation mode, some authors have suggested that deposition in the oral cavity during exhalation could be neglected [31] while others have proposed that inspiratory and expiratory deposition efficiency in the oral cavity should be considered to be equal [7,45]. Our data suggest the latter to be a reasonable hypothesis, at least for micron-sized particles. Indeed, there was no significant difference between inspiratory and expiratory deposition in the oral cavity for this size range ( Figure 5). This is also supported by data from Verbanck et al. [43], albeit performed in a single oral airway geometry, that showed almost identical oral deposition curves for inspiration and expiration for experiments performed with 3 and 6 µm particles.
In terms of regional deposition, intrasubject differences were observed when the airflow direction was reversed. While during inspiration, hot spots were mainly found in the oropharyngeal region and also in the mouth for the largest particles (10 µm and up), high deposition was preferentially found in the laryngeal region during expiration with minimal deposition in the mouth cavity for large values of d 2 a Q (Figures 4 and 6). When averaged over all subjects, these differences were statistically significant for particles ≥ 5 µm. This may be of minimal clinical relevance for the largest particles (>20 µm) that tend to have high intrathoracic deposition rates, leaving only a negligible particle fraction, if any, to be exhaled. In contrast, different deposition patterns may be of importance when assessing the side effects of an inhaled drug with particle size distribution in the range of 1-10 µm, typical of most pharmaceutical aerosols. For example, the delivery of inhaled corticosteroids (ICS) is highly effective at controlling the inflammatory component of chronic airway disease such as in asthma and COPD, while limiting systemic toxicities [46]. However, the efficacy of ICS is highly dependent upon its ability to bypass the upper airway. High extrathoracic deposition not only limits the amount of drug that can reach the lungs but can also result in unwanted local side effects. Indeed, the repeated deposition of corticosteroids in the larynx can cause a wide variety of clinical side effects including hoarseness, sore and/or dry throat, dysphonia, and candidiasis [12,47].

Comparison of Whole-Lung Deposition with Experiments
In an attempt to reproduce aerosol exposure studies performed on the same subjects from which the upper airway geometries were obtained, we predicted deposition for breathing maneuvers similar to that used in the experiments [20]. To do so, we coupled our predictions of oral deposition to an intrathoracic deposition obtained with an improved MPPD model [30] that also accounted for subject-specific lung volumes and subject-specific inhalation and exhalation flow rates ( Table 2). As the MPPD model assumes uniform ventilation among the different regions of the lung, we limited our comparison to whole-lung deposition data obtained in healthy subjects for which, unlike in COPD subjects, a uniform ventilation distribution is a reasonable assumption. Comparison of our MPPD/CFD predictions with experimental values shows relatively good agreement ( Figure 7A) as did the comparison between experimental data and MPPD predictions coupled with the Stahlhofen equation for oral deposition (MPPD/empirical predictions, Figure 7B). As the Stahlhofen equation only incorporates particle size and flow rate characteristics, less scatter was observed in the MPPD/empirical predictions than in the MPPD/CFD data, the latter also reflecting the effect of upper airway geometry on oral deposition (Figure 8).
There are a few limitations worth noting that could have affected our predictions. First, CFD simulations were performed in upper airway geometries with rigid walls that were based on CT images obtained at the end of a one-liter inspiration. Thus, the effect of any variation in the upper airway geometry that occurs during tidal breathing even in healthy subjects with no upper airway pathology [48,49] was neglected. Second, as CT imaging and aerosol studies were performed on two different occasions, it is highly probable that the position of the tongue within the oral cavity changed between the sessions. Furthermore, the position of the tongue was not controlled for during the aerosol studies and its position may well have moved during the experiment. This would affect deposition predictions in the oral cavity. Indeed, previous studies have shown that the position of the tongue significantly affects the delivery of aerosol in the trachea, highlighting the high intrasubject variability in aerosol delivery to the lungs [50][51][52].

Conclusions
Drug inhalation is a mainstay in the management of respiratory diseases. Its success does not only depend on the pharmacology of the drugs being inhaled but also on the site and extent of deposition in the respiratory tract. Both the physical properties of the inhaled aerosols and the subject characteristics (i.e., lung volume and geometry, breathing pattern, disease) strongly affect deposition. The upper airway is characterized by variable cross-sectional areas and sharp changes in a flow direction that is conducive to deposition by inertial impaction. Most pharmaceutical aerosols target the intrathoracic lung and as such, deposition in the oral passages should be kept to a minimum to circumvent local adverse effects and to maximize intrathoracic lung dose.
Deposition of 1-30 µm particles was predicted in eleven models of oral airways of adults at two flow rates, one flow typical of tidal breathing/nebulizer use and one flow in the range recommended when using a DPI. For an optimal size distribution of 1-5 µm for pharmaceutical aerosols, our data suggest that >75% of the inhaled aerosol is delivered to the intrathoracic lungs in most subjects when using a nebulizer but only in about half the subjects when using a DPI. This is also the first report of in silico predictions in such many geometries for both inhalation and exhalation. Averaged over all geometries, our data showed no significant difference in overall deposition efficiency in the oral cavity between the inspiratory and expiratory phases. In contrast, subregional patterns largely differed between the two phases, with areas of hot spots preferentially in the oropharynx during inspiration, and in the laryngeal region during expiration. Funding: This study was funded by the NIEHS at the NIH grant U01 ES028669.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki. Acquisition of human data was approved by the Institutional Review Board of the University of Washington, Seattle (Protocol 40712, initial date of approval: 23 October 2012). Use of de-identified data was approved by the University of California, San Diego (protocol 120661, latest date of approval: 15 June 2022).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on reasonable request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on reasonable request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
Appendix A-Individual geometries Figure A1. Upper airway geometry of all individuals. The same scale is used in all models. Table A1 lists the volume and surface area of each upper airway geometry shown in Figure A1.  Figure A1. Upper airway geometry of all individuals. The same scale is used in all models. Table A1 lists the volume and surface area of each upper airway geometry shown in Figure A1.