A Combined Computational and Experimental Study on Nanoparticle Transport and Partitioning in the Human Trachea and Upper Bronchial Airways

Received: January 12, 2020
Revised: April 10, 2020
Accepted: June 16, 2020

 Copyright The Author(s). This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are cited.

Cite this article:

Ma, J., Su, W.C., Chen, Y., Shang, Y., Dong, J., Tu, J. and Tian, L. (2020). A Combined Computational and Experimental Study on Nanoparticle Transport and Partitioning in the Human Trachea and Upper Bronchial Airways. Aerosol Air Qual. Res. 20: 2404–2418. https://doi.org/10.4209/aaqr.2020.01.0015

HIGHLIGHTS

• Nanoparticle transport and partition to airway morphology was investigated.
• Small-scale wall shear fluctuation is only seen in the anatomical model.
• Nanoparticle mixing in the anatomical model was significantly more efficient.
• Flow and particle partition in lung lobes were not necessary in proportion.
• Simulated flow partition among 5 lobar airways was dependent on boundary condition.

Keywords: Airway morphology; Nanoparticle transport; Deposition and partitioning; Tracheobronchial tree; Computational modeling; Experimental measurement.

INTRODUCTION

Aerosol transport and deposition in human respiratory airways has been an important research area in past few decades. The knowledge of particle trajectory, particle-flow interaction and deposition site is of significant value to respiratory toxicity study, and also to the therapeutic design of efficient inhalable drug delivery devices (Kleinstreuer and Zhang, 2003; Asgharian and Price, 2007; Inthavong et al., 2010; Longest and Holbrook, 2012; Gu et al., 2019; Tian et al., 2019). The internal flow problem requires appropriate reconstruction of the bounding surface that most realistically represent the airway geometry. For human extra-thoracic airways (nasal cavity/pharynx/larynx), anatomical models were frequently rebuilt via image scans or cadaver casts (Cheng, 1999; Tu et al., 2013); however, such methodology is not practical to reconstruct human whole lung which contains millions of airway branches with terminal dimension drops to hundreds of micrometers. Instead, whole lung modeling can only be achieved through mathematical approximation guided by key indices obtained from airway samples. Weibel (1963) proposed such a whole lung model (23 generations) completely assembled by sequence of 2D symmetric bifurcations based on selected sampling of non-randomly distributed structures. Mathematical equations for human lung morphometry were prescribed to quantitatively correlate the anatomy to its physiological functions. Based on the principles of: 1) volumetric proportionality between airway branches and the delivered pulmonary space; 2) homogeneous distribution of airway terminal branches, Kitaoka et al. (1999) formulated deterministic algorithms for 3D human tracheobronchial tree modeling with morphometric characteristics very similar to those obtained from airway cast measurements. The algorithms were further packaged into a software which can be freely downloaded for the purpose of studying lung anatomy, physiology, and pathophysiology (Kitaoka et al., 2013).

Apart from whole lung modeling, significant research was conducted in regional airways, where sectional details of the transport processes, critical to the understanding of underlying physics, can be analyzed. Focused on tracheobronchial cartilaginous rings, Martonen et al. (1994) investigated the morphology-induced variation on fluid dynamics in a simplified 2D airway passage. Using a symmetric 2D bifurcation model, Li and Ahmadi (1995) simulated the micro- and nanoparticle dispersion and deposition in human trachea and the main bronchi, where main deposition mechanisms were identified. Comer et al. (2001) performed the computational investigation of flow structures and particle deposition patterns in a double-bifurcation airway model (G3–G5). Zhou and Cheng (2005) measured the regional particle deposition efficiency and deposition patterns in a human airway cadaver replica including the oral cavity, pharynx, larynx, trachea, and four generations of bronchi. Based on a single bifurcation element (Heistracher and Hofmann, 1995; Phillips and Kaye, 1997), Tian and Ahmadi (2012) proposed a multi-level asymmetric lung bifurcation model (G0–G3), flexible to be extended, for airflow and particle transport investigation. More recently, up to the 6^th, 7^th and 15^th generation of anatomical tracheobronchial tree models were used by Koullapis et al. (2016), Lintermann et al. (2017), Frederix et al. (2018), Xu et al. (2019) and Dong et al. (2019a, b) for regional aerosol deposition studies. Clearly seen from these research investigations, both anatomical reconstruction and the mathematical simplified modeling approaches were engaged. While significant insights were rendered from these pioneer studies, fundamental analysis seeking morphometry-induced variations in flow pattern, particle trajectory and particle deposition statistics in human tracheobronchial airways were rare. In particular, intelligence of particle and flow partition at the airway bifurcating sites, critical to address navigated targeting delivery was not fully investigated.

To fill the gap, Shang et al. (2018) performed a numerical comparison of nanoparticle transport and deposition in human upper tracheobronchial airways (trachea and the first bifurcation truncated at five lobar partitions) by using an anatomical realistic reconstruction (through CT scans) and a mathematical simplified airway model. The investigation identified major irregularity (at head airway connection)-induced variation in airflow and particle deposition patterns in the trachea, while particle fate in the bronchi were less affected. The findings correlated well with that of Zhang and Finlay (2005), though their experiments focused more on effect of cartilaginous rings and particles in the micrometer range. Accounting for this geometric irregularity at head airway connection, a modified idealized airway model (based on that of Shang et al., 2018), closer to morphometry of the anatomical construction, was developed in the current study. Using the anatomical airway model (Shang et al., 2018) as a benchmark, new boundary considerations mimicking lobar bronchi air partition were applied to the newly developed idealized model, and major trachea-irregularity-induced enhancement of nanoparticle deposition was quantified. In particular, flow and particle partitioning in the branching airways, an important postulation for aerosol distribution in whole lung dosimetry modeling (Phalen et al., 2016), was examined, and the detailed flow and particle profiles at five lobar exists were analyzed.

A paralleled experimental measurement of sodium chloride (NaCl) nanoparticle penetration and deposition rates in the corresponding pair of anatomical and idealized airway models was developed. The computational model predictions were compared with the experimental measurement, general agreement and valuable information were obtained. It should be noted that, due to instrumentation limitation, current experiment only included nanoparticles of sizes greater than 10 nm. However, to obtain the relationship of particle partition, deposition and penetration characteristics with respect to particle sizes, current computational simulation included a wider range of nanoparticles (1–100 nm) assuming negligible agglomeration.

Current research contributes to the understanding of morphometry sensitivity of flow and particle trajectories in human tracheobronchial airways. It is of significance in context of the increased use of tracheobronchial airway modeling in particle study in human lung. In addition, particle partitioning, an important feature for navigated delivery of nanoparticles in airway passages, was investigated. Lastly, the paralleled experiment provided valuable validation data for current and future computational studies.

 
METHOD

 
Human Upper Tracheobronchial Airway Geometric Models

An anatomical reconstruction of an adult human tracheobronchial tree up to the 15^th generation (Fig. 1(a)) was made via CT scans (Shang et al., 2018). A truncation, containing trachea and the main bronchi, was extracted and used for a morphometry comparative study with an idealized airway model developed by Shang et al. (2018). Extending from the study, a modification was made to the existing idealized model by introducing the major irregularity at head-airway connection mimicking that of the anatomical reconstruction (Fig. 1(b)). As shown in the figure, the sectional upper tracheobronchial airway terminated at entrance of the five lobes leading to deeper lung. Further details on model specification and development can be found in the work of Shang et al. (2018).

Fig. 1. Human airway geometric modeling: (a) 15^th generation anatomical tracheobronchial tree model from CT scans; (b) truncated anatomical model containing trachea and the first bifurcation; (c) mathematical idealized model.

 
Computational Mesh

Computational mesh was generated by using Ansys Fluent Meshing (18.0), where 789,310 and 750,134 polyhedron elements of dimension 0.5 mm were created for the anatomical and idealized models as shown in Fig. 2. With the first grid point located 0.05 mm away, 10 prism layers were attached to the airway surface to resolve near-wall features. The prism layers contained hexahedron cells, growing at a rate of 1.2 from the boundary to the core region. Mesh independent study was performed on the newly developed idealized model, where velocity profiles along the line and across a cross section near the airway bifurcation (Fig. 3(b)) were obtained and compared by using mesh resolutions of 0.7, 0.6, 0.5 and 0.4 mm, mesh sizes of 370,588, 512,786, 750,134 and 1,189,719 respectively. It was shown in Fig. 3 that velocity profiles along the selected line and cross section were insensitive to the four mesh resolutions, while Fig. 3(b) implied the optimal mesh size of 750,134 with mesh dimension of 0.5 mm be used for the computational study.

Fig. 2. Computational mesh with polyhedron meshing schemes.

Fig. 3. Mesh independent study: (a) velocity profiles along the line of a cross section near the bifurcation; (b) area-averaged velocity magnitude on the cross section near the bifurcation.

 
Fluid Flow Simulation

Steady inhalation (30 L min^–1), assuming the particle deposition mainly occurred during inhalation phase (Inthavong et al., 2010), was employed. Breathing air was predominantly turbulent due to the presence of strong irregularities at the laryngeal (Tian and Ahmadi, 2012). Reynolds stress transport model together with the “two-layer zonal” boundary condition was used to reproduce the breathing air, as it has been shown to be the most accurate in resolving near-wall flow features, critical to micro- and nanoparticle deposition studies (Tian and Ahmadi, 2007). Details on the governing equations of turbulence mean flow, turbulence fluctuations, and turbulence near-wall modeling can be found in the work of Tian and Ahmadi (2007), and Shang et al. (2018).

Ansys Fluent 18.0 was used to perform the steady flow simulation. In physically realistic breathing cycles, airflow was driven by the pressure gradient between ambient environment and the thoracic cavity. Therefore, prescribed pressure difference between airway entrances and the exits was frequently enforced as the boundary condition in breathing airflow simulations (Koullapis et al., 2016; Shang et al., 2018). In addition to this approach, current study took an additional boundary consideration where assumed flow partition among the five lobar branches was prescribed. As the regional airway was distantly upstream the chest diaphragm where the exact pressure at each branching exists was unknown, the new boundary consideration might be more physically appropriate. Comparisons of particle simulation results were made by using the two different boundary considerations.

Particle Simulation

Assuming particles are sufficiently dilute that the airflow field is not affected by the presence of particles, and particle-particle collision is not considered, the governing equation of motion for each particle is given as:

where u_p is the particle velocity vector, t is the time, g is the gravitational acceleration, ρ_p and ρ are the particle and fluid densities respectively. As the current study is focused on nanoparticles, gravitational and the buoyancy forces can be neglected. F_D is the particle-flow slip force given by 18µ(u_p – u)/(d²ρ_p), with d being the particle diameter. In Eq. (1), C_c is the Cunningham correction for nanoparticles which is given as:

where λ is the gas mean free path. In Eq. (1), F_L is the Saffman lift force, given as:

where K = 2.594 is the constant coefficient (Saffman, 1965), ν is the fluid kinematic viscosity, S is the particle-to-fluid density ratio, u and d_ij are the fluid velocity and strain-rate tensor respectively. The Brownian diffusion force F_B in Eq. (1) is modeled as a Gaussian white noise random process with spectral density given by Li and Ahmadi (1992):

where k_b = 1.38 × 10^–23 J K^–1 is the Boltzmann constant and T is the absolute temperature.

Particle simulation was performed by using the discrete phase model (DPM) of Ansys Fluent 18.0. 10,000 monodispersed ultrafine particles from 1 to 100 nm were uniformly released from the trachea entrance. Particle trajectories, deposition rate and penetration statistics were calculated. It was assumed that particle deposited onto the airway surface upon touching.

 
Deposition Intensity Factor

To quantify particle deposition patterns in the tracheobronchial airway surface, the deposition intensity factor (DIF) is defined as the ratio of the number of deposited particles onto 1 mm² surface area of the local region to the total number of released particles:

The parameter is similar to a prior defined deposition enhancement factor (DE), formulated as the regional particle deposition within a defined area relative to the total deposition normalized by the surface area of the entire region (Balásházye et al., 1999; Zhang et al., 2005; Longest and Holbrook, 2012). The DIF is continuously defined across all surface area of the tracheobronchial airway model, where particle deposition density can be easily visualized. The choice of 1 mm² as the local region is determined by the optimal visual representation of the particle deposition results in the current study.

Experimental Setup

In this study, lung deposition experiments were conducted by delivering sodium chloride nanoparticles (10, 20, 50, and 100 nm) into the physical models of the airway (anatomical model and the mathematical idealized Model) (Fig. 4). The airway physical replicas were printed using a high-resolution 3D printer (LulzBot TAZ 6; Aleph Objects, Inc., Loveland, Colorado, USA) with conductive PLA filaments (Proto-pasta; ProtoPlant, Vancouver, WA, USA). By obtaining the penetration efficiency of the NaCl particles for all lung tubes on the airway replica, NaCl particle deposition fraction (%) in a certain airway bifurcation can be calculated. The experimental approach has been used previously to estimate lung deposition of airborne carbon nanotubes and graphene (Su and Cheng, 2015; Su et al., 2016). To facilitate the measurement of NaCl particle penetration efficiencies at a certain lung tube on the airway replica, a set of trimmed airway replicas were made similar to those shown in Su and Cheng (2015). Prior to the experiment, the surface of the airway replica was coated with a thin layer of silicon oil to prevent the resuspension of the deposited particles. To generate the nanoscaled NaCl particles (10, 20, 50, and 100 nm), an atomizer aerosol generator (3079A; TSI Inc., Shoreview, MN, USA) was employed using sodium chloride water solution (0.5% v/v). NaCl aerosols generated by the atomizer were immediately delivered to a series of drying columns to allow water mist in the NaCl aerosols to completely evaporate. Then, nanoscaled NaCl particles were charge neutralized and delivered into the airway replicas (untrimmed and trimmed airway replicas) for deposition measurement. It should be noted that, due to current technical limitation, generation and classification of NaCl nanoparticles below 10 nm is not possible.

Fig. 4. The experimental setup of nanoparticle lung deposition measurements.

Two units of nanoparticle sizer (SMPS+C; GRIMM Aerosol Technology, Germany) were used simultaneously, with one unit connected to the trachea inlet, and the second unit connected to the outlet of the specified lung tube along the airway replica. The two particle sizers were set at the same particle size channel to simultaneously measure the number concentration of NaCl particles at the specified size. Consequently, penetration efficiency at a specific nanoparticle size for a certain lung tube can be obtained. Measurements of NaCl particle penetration efficiency were carried out in turn for all lung tubes along the airway replica. Fig. 4 shows the experimental setup of nanoparticle lung deposition experiments. When penetration efficiencies of all lung tubes on the airway replica were acquired, the deposition fraction of the nanoscaled NaCl particles in each airway bifurcation can be reasonably calculated given air flow rate in each lung tube. Details of the approach (methods and equations) can be found in the work of Su and Cheng (2015) and Su et al. (2016).

RESULTS AND DISCUSSION

 
Breathing Airflow Pattern

Streamwise airflow in the anatomical and newly developed idealized airway models were displayed in Fig. 5. Main features of the pressure and velocity distribution conformed to each other in the two representations. Both captured the velocity acceleration and pressure drop in vicinity of the head airway connection due to bending and contraction. Airway contraction was more radical in the anatomical model that led to a higher pressure drop, flow acceleration and larger flow recirculation zone behind the bend. At carina ridge, elevated pressure was observed in both the anatomical and idealized models, corresponding to the flow stagnation region where partitioning of the breathing air took place before entering into the daughter branches. It was observed in Fig. 5 that bulk air penetrated into the inner walls and was more pronounced in the daughter branch with sharper branching angle (left bronchi). Higher level of mixing of the breathing air was observed in daughter branches of the anatomical model than that in the idealized one.

Fig. 5. Streamwise airflow patterns (pressure and velocity magnitude) at breathing rate of 30 L min^–1 in the central plane of the tracheobronchial airways: (a) truncated anatomical model; (b) mathematical idealized model.

Velocity vectors of the secondary flow at selected cross sections along the airway was shown in Fig. 6. Due to space limitation, selected sample locations were only demonstrated on the anatomical model in the figure. It was noted that the color map includes the streamwise velocity component, but length of the vectors only account in plane secondary-flow components. Corresponding to the observation in Fig. 5, level of the secondary flow in the anatomical model was relatively higher than that in the idealized representation throughout the airways. It was seen from Fig. 6 that secondary flow was mainly induced by the abrupt changes of airway morphology (e.g., bend, contraction, and the bifurcation), while effect of the boundary layer formation was insignificant. Higher level of mixing in the anatomical model, reflected by the more uniform secondary-flow distribution, might be attributed to the wavy wall roughness in the trachea due to cartilage rings and the deformation of airway passage from circular to highly oblate cross-sectional areas at transition from trachea to the daughter branches.

Fig. 6. Secondary flow pattern at cross sections along the trachea, left and right bronchi: (a) anatomical model; (b) idealized model.
 

Turbulence

Fig. 7 compared the turbulent kinetic energy at the breathing rate of 30 L min^–1 in the anatomical and newly developed idealized airway models. Significant turbulence production was observed at the trachea contraction in the anatomical model (Fig. 7(a)), where turbulence filled the entire region up to the carina ridge. A secondary production, at much smaller scale, was induced by flow partition immediately following the bifurcation. In the main channels (trachea and daughter branches), turbulence was continuously dissipated as the airflow traveled downstream (Fig. 7(a)). When looking at the idealized model (Fig. 7(b)), main feature of the turbulence activity was captured, however to a lesser degree where trachea contraction was introduced near the entrance. The generated turbulence was seen to propagate to a smaller region, and dissipated more quickly than that in the anatomical model. In both models, left bronchi hosted slightly higher-level turbulence activities than that in the right bronchi.

Fig. 7. Turbulence kinetic energy at breathing rate of 30 L min^–1 in central plane of the tracheobronchial airways: (a) truncated anatomical model; (b) mathematical idealized model.

Wall Shear Stress

Fig. 8 displayed the wall shear stress of the airway surfaces mapping to a 2D domain for superior visualization. Details on surface mapping technique and the applications can be found in the work of Inthavong et al. (2014), Tian et al. (2017a, b), and Shang et al. (2018). As shown in the figure, global high-shear zones agreed in both models, which were at the trachea entrance and in the vicinity of the airway bifurcations. Due to bending and contraction near trachea entrance, both models experienced high level of shear stress in upper sections of the airways (Figs. 8(a) and 8(b)). Significant difference between the anatomical and idealized models was the small-scale local wall shear fluctuation, which was observed throughout the anatomical airway surface (Fig. 8(a)), however, not seen in the idealized model (Fig. 8(b)). Very low and uniform distribution of the wall shear stress in the trachea and bulk region of the daughter branches in the idealized model (Fig. 8(b)) implied a smooth airway surface with low friction. On the contrary, sporadic and streak-patterned wall shear stress in the anatomical model (Fig. 8(a)) implied a rough surface corresponding to the cartilage rings and predominantly wavy interior wall. Wall shear distribution was closely related to the formation of near-wall eddies that contribute to particle deposition processes.

Fig. 8. Wall shear stress on the unwrapped 2D surface models of the tracheobronchial airways: (a) truncated anatomical model; (b) mathematical idealized model.

Particle Deposition Pattern

To visualize particle deposition pattern, DIF (Eq. ((5)) of 10,000 uniformly released monodispersed unit density solid particles (1 and 100 nm) were calculated and presented in Fig. 9. DIF quantified the local particle deposition density, and facilitated the identification of deposition hot spot and variations in the domain. For the 1 nm particles, deposition pattern in the newly developed idealized model (Fig. 9(b)) was similar to that in the anatomical model (Fig. 9(a)), where elevated deposition was observed sporadically in vicinity of the bending near trachea entrance. However, the deposition “hot region” in the anatomical model exhibited higher intensity and covered a larger area. For the 100 nm particles, neither the anatomical nor idealized model displayed deposition enhancement induced by the trachea bending, however, the slightly higher deposition throughout the airway surface as seen in Fig. 9(a) could be attributed to the local wall roughness (e.g., cartilage rings) in the anatomical model. It was clearly shown in Fig. 9 that regional deposition of nanoparticles in the upper tracheobronchial airways was size dependent. 1 nm particles were more likely to be deposited in the trachea, while significant fewer deposited in the main bronchi. For the 100 nm particles, no distinctive deposition preference was observed in both models, while slight higher deposition intensity was seen near the bifurcation ridge, especially in the idealized model (Fig. 9(b)). Still, slight higher level of deposition of 100 nm particles in the trachea was observed in the anatomical model than that in the idealized model.

Fig. 9. Particle deposition intensity factor (DIF) for 1 and 100 nm particles in the 2D unwrapped surface models: (a) truncated anatomical model; (b) mathematical idealized model.

Finally, a comparison between Figs. 8 and 9 showed a remarkable resemblance between the airway wall shear stress and particle DIF for the 100 nm particles. This was not seen with the 1 nm particles, indicating the presence of a molecular-diffusivity-dominated transport process other than a transport mechanism dictated by the particle-flow inertia interaction as seen in the earlier case.

 
Particle Penetration Profile

Figs. 10 and 11 displayed the penetrating profiles of 1 and 10 nm particles at the five bronchi exits before entering into deeper lung. For both sized particles, penetrating profiles varied in the anatomical and idealized models. In general, improved mixing of the inhaled nanoparticles was observed in the anatomical model, where particles were seen to disperse to all regions at the cross sections though preferential concentration might be exhibited. This was not observed in the idealized model (Figs. 10(b) and 11(b)), where all cross sections had space void of inhaled nanoparticles. In particular, large void space, free from particles, was clearly seen in the right middle (RM) and right upper (RU) exits.

Fig. 10. Particle distribution (1 nm) at 5 bronchi exits (LL: left lower lobe; LU: left upper lobe; RL: right lower lobe; RM: right middle lobe; RU: right upper lobe): (a) truncated anatomical model; (b) mathematical idealized model.

Fig. 11. Particle distribution (10 nm) at 5 bronchi exits (LL: left lower lobe; LU: left upper lobe; RL: right lower lobe; RM: right middle lobe; RU: right upper lobe): (a) truncated anatomical model; (b) mathematical idealized model.

A unique vortex-like penetrating profile was observed for the 10 nm particles as they exited from the left upper (LU) branch in both models (Fig. 11). This indicated the evolution of a vortex-like secondary flow which lead to the left upper lobe. An interesting observation was that the 1 nm particles did not respond to this unique secondary-flow pattern and particle distribution remained random across the LU exits (Fig. 10). It was also noted that level of inhomogeneity in particle distribution was increased when particle size grew from 1 to 10 nm (Figs. 10 and 11). Apparently, particles responded more to the flow with increased particle inertia, and preferential particle trajectories were starting to be formed. It was also noted that the largest discrepancies (in particle penetrating profile) between the idealized and anatomical models were seen in the RM and RU branches.

 
Particle Deposition Efficiency

Fig. 12 displayed the total and regional deposition efficiency of the nanoparticles entering from trachea of the tracheobronchial airway models. Particle deposition efficiency is defined as the ratio of the deposited particles in a region to the total number entering the domain. To gain further insight of particle-flow relationship, one more set of simulation were performed in the idealized model, with new boundary condition assuming prescribed flow partition in the five lobes following that of the anatomical model (Table 1). Simulation results from Shang et al. (2018), by using a prior developed idealized model without accounting for the geometric irregularity at head-airway connection, was also included for comparison. Experimental measurement on nanoparticle deposition in human tracheobronchial airway replicas in current and prior studies were plotted along with the computational simulations. It was shown in Fig. 12(a) that the simulated nanoparticle total deposition efficiency was the highest in the anatomical model (red -o-). By accounting for the geometric irregularity at head-airway connection, particle depositions in the newly developed idealized model (blue -□-, -◊-) were very close to that in the anatomical model. This was a significant improvement from simulation predictions without accounting for this geometric irregularity (green -*-, -☆-; Shang et al., 2018). However, slightly lower deposition was still seen in the newly idealized model (Fig. 12(a)), which could be attributed to the added deposition due to small-scale surface roughness in the anatomical model (Fig. 1).

Fig. 12. Total and regional particle deposition efficiency in the airway models: (a) total deposition; (b) deposition in the trachea; (c) deposition in the left bronchi; (d) deposition in the right bronchi.

A subtle yet important transition occurred at about 20 nm (Fig. 12(a)), where the deposition curve of the new idealized model deviated from that of the anatomical model, and approached to that of the prior idealized model (Shang et al., 2018), when particle size continued grew. This implied that, for larger nanoparticles (d > 20 nm), deposition was less sensitive to macroscale irregularities, however more affected by the small-scale surface roughness. The enhancement of particle deposition due to bending and contracting near the tracheal entrance gradually phased out for larger particles (d > 20 nm). Another interesting observation in Fig. 12(a) that total deposition in the truncated airways was insensitive to the prescribed boundary conditions (pressure driven versus prescribed flow partitioning among five lobes). Deposition was mainly governed by the global structure of the geometry and airflow.

Finally, the simulation agreed well with both current and prior experimental data in trend and magnitude. In particular, the measured total deposition efficiency in the anatomical reconstruction (Fig. 12(a), red +) was generally (except at 50 nm) slightly higher than that in the measurement in the idealized model (blue x), which supported the simulation results. It should be noted that the slight higher prediction of the deposition in current simulation and the experiment than that in the prior experiments (Cohen et al., 1990; Smith et al., 2001), was most likely to be caused by inter-subject variations in the two studies. Reported experimental data on nanoparticle deposition were extremely limited for size smaller than 10 nm. Current experimental capability also prohibited the characterization and measurement of nanoparticles smaller than 10 nm. Further experimental research in extreme submicrometer particle sizes (< 10 nm) are greatly in need.

Deposition efficiencies in the trachea (Fig. 12(b)) were similar to that of the total deposition in tracheobronchial airway models. Anatomical model captured the most inhaled nanoparticles, idealized model without accounting for the irregularity at head-airway connection (Shang et al., 2018) captured the least and the newly developed idealized model was in between. Key observations, such as deposition insensitivity to boundary conditions (pressure driven versus prescribed flow partitioning among five lobes) and the transition of deposition mechanisms at about 20 nm, were reproduced in Fig. 12(b). Nanoparticle deposition in the left and right bronchi was significantly lower than that in the trachea (5–8% versus 25–45%), and deposition efficiency was about the same in both the anatomical and idealized models (Figs. 12(c) and 12(d)). It should be noted that the lower level of deposition observed in the left bronchi of the anatomical model for d < 3 nm was most likely caused by fewer particles entering the branch due to high deposition in prior pathway (trachea). Fig. 12(c) and 12(d) indicated that the geometric variation at the trachea entrance had negligible effect towards nanoparticle deposition downstream in the first-level bronchi. Once again, current experimental measurement agreed well with the simulation in trend and magnitude in regional depositions in the trachea, left and right bronchi (Figs. 12(b), 12(c) and 12(d)). A consistently slight higher deposition in the bronchi (Figs. 12(c) and 12(d)) was measured in the anatomical replica than that in the replica of the newly developed idealized model, however, this was not perceived in the trachea (Figs. 12(b)).

Particle and Flow Partition in the Main Bronchi and the Five Lobar Branches

Fig. 13 examined the particle and flow partition (% distribution) in the main bronchi. To gain further insight of particle-flow partitioning relationship, one more set of simulation were performed in the idealized model of Shang et al. (2018) and newly developed idealized model (Figs. 13(d) and 13(e)), with new boundary condition assuming prescribed flow partition in the five lobes following that of the anatomical model (Table 1). Fig. 13 were plotted on the same scale for optimal comparison. Clearly seen in Figs. 13(b) and 13(c) that flow driven by the pressure gradient in the idealized models incurred different flow partitioning when compared to that in the anatomical model (Fig. 13(a)). It was also noted from Fig. 13 that, even with prescribed volume flow rate following the flow partitioning in the anatomical model, particle partitioning in all configurations was distinctive in magnitude as well as in patterns. Striking match between flow and particle partitioning was observed in Fig. 13(e) (newly developed idealized model with new flow boundary condition assuming prescribed flow partition in five lobes), however, it was not reproduced in any other configurations. In general, flow and particle partitioning in the main bronchi was not necessary in proportion, and percentage distribution of nanoparticles in the two bronchi was size dependent up to 10 nm (except for the newly developed idealized model with new boundary condition; Fig. 13(e)) and then it approached a steady state (10–100 nm). Nanoparticle partitioning in the two main bronchi mostly varied between 4 to 10%, and in the extreme case of 20% (Fig. 13(e)). Least variation or equal particle partitioning was observed in the anatomical model and in the pressure-driven newly developed idealized model (Figs. 13(a) and 13(c)).

Fig. 13. Flow and particle partition in the main bronchi: (a) anatomical model; (b) idealized model of Shang et al. (2018); (c) newly developed idealized model; (d) idealized model of Shang et al. (2018) with new flow boundary condition assuming prescribed flow partition among five lobes; (e) newly developed idealized model with new flow boundary condition assuming prescribed flow partition among five lobes. 

Experimental measurement of particle partition in the two main bronchi, at the specified branching volume flow rate (Table 1), was presented in Figs. 13(a) and 13(e). For the size of particles considered (10, 20, 50 and 100 nm), particle partition was seen to closely resemble that of the flow in both anatomical and the newly developed idealized airway replicas. The measured particle partition highly agreed with simulation results in the newly developed idealized model (assuming prescribed flow rate in five lobes; Table 1, Fig. 13(e)); however, measured partition in the anatomical replica did not agree with the simulation (Fig. 13(a)), where particles were almost equally parted into the two bronchi. The discrepancies might be related to the inhomogeneous particle penetration profiles (Figs. 10 and 11), which were not accounted for in the experiment with only one single reading of particle number concentration at each lung tube outlet. The discrepancy in Fig. 13(a) could also be related to boundary conditions, where the simulation was initiated by pressure difference while enforced flow partitioning was applied in the experiment. It was clearly shown in Fig. 13 that particle partition was most likely affected by both flow partitioning and particle penetration profiles. Fig. 13 also implied the sensitivity of particle partition to boundary conditions in the simulation. Further research is needed, in both the computational and experimental domains, to identify the most physical realistic boundary condition to be applied in particle partition study to recover the realistic transport processes.

Fig. 14 displayed the particle and flow partition (% distribution) in the five lobar branches. Similar to Figs. 12 and 13, one more set of simulation were performed in the idealized model of Shang et al. (2018) and newly developed idealized model (Figs. 14(d) and 14(e)), with new boundary condition assuming prescribed flow partition in the five lobes (Table 1). All subplots in Fig. 14 were on the same scale for optimal comparison. Clearly seen in Fig. 14, both the flow and particle partitioning in the RM branch was the lowest (5–10%) among all configurations. While flow partition in LU was the next lowest (~15–20%), particle partition in LU was relatively high in all simulations. In all configurations, particle partition in LU was over 25%. Agreeable trend of the particle and flow partition was also observed in the RU branch, where the highest percentage of the particle penetration correlated to the highest volume flow rate (Figs. 14(a), 14(d) and 14(e)). However, it should be noted that with pressure gradient boundary condition (Figs. 14(b) and 14(c)), volume flow rate in the RU branch was about 20%, roughly 10% lower than that in the anatomical model, yet particle partition was still among the highest in all branches. While particle partition in the left lower (LL) and right lower (RL) branches were among the lowest (< 20%), the volume flow rate was in the medium range in the anatomical model or the highest in the idealized models with pressure gradient boundary condition.

Fig. 14. Flow and particle partition in the five lobar branches: (a) anatomical model; (b) idealized model of Shang et al. (2018); (c) newly developed idealized model; (d) idealized model of Shang et al. (2018) with new flow boundary condition assuming prescribed flow partition among five lobes; (e) newly developed idealized model with new flow boundary condition assuming prescribed flow partition among five lobes.

Particle partition in five lobar branches were experimentally measured in the anatomical model and the newly developed idealized model (Figs. 14(a) and 14(e)). The experimental data were straightforward that all partition followed exactly the flow distribution with only minor deviations. The measurement supported simulation observation that the highest and lowest partition occurred in the RU and RM lobes respectively. However, simulated particle partition in other lobar branches were not in proportion to flow distribution as what was observed in the experiment. As pointed out earlier, with only one reading of the particle number concentration at each lung tube outlet, rather than accounting for the inhomogeneous particle penetrating profile as shown in Figs. 10 and 11, might attribute to the discrepancy between experiment and the simulation.

Important derivations from simulation results (Fig. 14) were:

• Flow partition among five lobar airways was dependent on airway boundary conditions. Flow driven by pressure gradient or flow driven by prescribed volume flow rate could yield different flow partition in subsequent airway branches.
• Particle partition was relative consistent in the five lobar airways. Independent of geometric modeling and boundary consideration in current study, the highest to lowest particle partition in the five lobes was RU, LU, LL, RL, and RM, respectively.
• Upper lobar branches (RU and LU) attracted significantly higher portion (> 70%) of the inhaled nanoparticles into deeper airways.
• Minimum nanoparticles (≤ 5%) passed through the RM (right middle) branch.
• Flow and particle partition in lung lobes were not necessary in proportion. Flow partitioning information cannot be used to predict nanoparticle partition in lung lobes.
• Partitioning of nanoparticles in the five lung lobes were size dependent up to 10 nm. Beyond that, partitioning was steady and size independent.

Important derivations from experiment (Fig. 14) were:

• Particle partition among five lobar airways closely follow flow distribution and independent of geometric models.

In summary, particle partition among lung lobes and the downstream branches contains valuable information for the design of active controlled medical devices for particle delivery in human respiratory systems. Further research is needed to identify a clear relationship between particle and flow partitioning in respiratory airways. The effect of inhomogeneity of flow and particle distribution along lung airways towards particle partition transport needs further investigation.

CONCLUSIONS

This study employed an anatomically realistic model of a human tracheobronchial airway as well as a newly developed idealized extension including the trachea and two main bronchi to investigate the transport and partitioning of nanoparticles and the influence of airway morphology. Experimental measurements were also conducted in parallel. Based on our results, we drew the following conclusions:

• To accurately simulate the transport of nanoparticles (d < 20 nm) in the human trachea, the major geometric features of the airway must be incorporated, although those of the secondary bronchi exert less influence.
• Local small-scale wall shear fluctuations were observed across the surface of the airway in the anatomical model but not the idealized model; this phenomenon may be related to the higher deposition efficiency of particles in the former, which was especially prominent for the fraction larger than 20 nm. The physical mechanisms underlying the relationship between small-scale superficial roughness and particle deposition in the human airway require further investigation.
• The anatomical model exhibited lower shear stress and lower deposition at the bifurcation site (for larger nanoparticles), suggesting that transport is more efficient in the anatomical design than the idealized. Further investigation is necessary to determine whether this difference is related to the extreme irregularity of the airway geometry.
• While traveling downstream, the nanoparticles were far more thoroughly mixed in the anatomical model than the idealized model. The effect of physical properties that are associated with this phenomenon, such as the airway geometry, small-scale superficial roughness of the channel, and bifurcation structure, requires further investigation.
• A stronger response to the local flow patterns of larger nanoparticles in certain regions, such as the upper left bronchial exit, which displayed vortex-like nanoparticle penetration, may shed light on enhancing the efficiency of targeted delivery systems.
• Contrary to the experimental results, the particle and flow partitions in the lung lobes were not necessarily proportional in the numerical simulation. A clear relationship between these partitions must be established through additional research before using flow partitioning data to predict nanoparticle partitioning in this region.
• In the simulation, the flow partitioning in the five lobar airways was dependent on the boundary conditions, whereas the particle partitioning was relatively consistent, regardless of the boundary conditions or modeled geometry. The largest particle partition was found in the upper right lobe, followed by the upper left, lower left, lower right, and center right lobes.
• The nanoparticle partitioning in the five lung lobes was dependent on the particle size when the particles were 10 nm or less in diameter. With larger particles, the partitioning remained steady, regardless of the particle size.
• The upper (right and left) lobar branches attracted a significantly higher percentage (> 70%) of the inhaled nanoparticles.
• The deposition of larger (d > 20 nm) nanoparticles was less sensitive to large-scale geometric variations, such as bends and contractions near the tracheal entrance, than to the small-scale superficial roughness of the passage.
• The deposition of nanoparticles was significantly lower in the left and right bronchi (5–8%) than the trachea (25–45%). The deposition efficiency was approximately equal between the models (Fig. 12(c)). The geometric variation near the tracheal entrance produced only negligible effects on the downstream nanoparticle deposition in the first-level bronchi.

ACKNOWLEDGEMENTS

The financial supports provided by the National Natural Science Foundation of China (Grant No. 91643102) and the Australian Research Council (Grant Nos. DE 180101138 and DP160101953) are gratefully acknowledged.

Aerosol Air Qual. Res. 20 :2404 -2418 . https://doi.org/10.4209/aaqr.2020.01.0015