An In Silico Subject-Variability Study of Upper Airway Morphological Influence on the Airflow Regime in a Tracheobronchial Tree

Determining the impact of inter-subject variability on airflow pattern and nanoparticle deposition in the human respiratory system is necessary to generate population-representative models, useful for several biomedical engineering applications. Thus, the overall research objective is to quantitatively correlate geometric parameters and coupled transport characteristics of air, vapor, and nanoparticles. Focusing on identifying morphological parameters that significantly influence airflow field and nanoparticle transport, an experimentally validated computational fluid-particle dynamics (CFPD) model was employed to simulate airflow pattern in three human lung-airway configurations. The numerical results will be used to generate guidelines to construct a representative geometry of the human respiratory system.


Introduction
Configurations and dimensions of human respiratory systems may vary significantly among subjects, thereby influencing airflow as well as aerosol transport and deposition. Thus, it is of interest to assess the degree of inter-subject variability to develop population-representative lung-airway models with beneficial applications in toxicology and drug delivery. Such a subject-variability study can establish relationships between morphological dimensions and airflow characteristics, to guide drug inhalation for individuals with different lung physiological characteristics, as well as drug inhaler design with accurate dose delivery. For example, dry-powder inhalers (DPIs) and pressurized metered-dose inhalers (pMDIs) are the most widely used devices for the treatment of asthma and chronic obstructive pulmonary diseases (COPD). In particular, a pMDI is convenient to use, inexpensive, and allows for multi-dose applications [1]. Moreover, the analysis will generate insight into the physiology of the lung and assist in the development of non-invasive diagnostic tools to treat lung diseases by comparing geometric differences between healthy and diseased configurations. The methodology could also be extended to variability studies between humans and animals. Such inter-species variability investigations are essential to provide more realistic scale-up correlations based on different levels of particle deposition in human vs. animal upper lung airways.
Choi et al. [2] pointed out that there are two factors which significantly affect the airflow between subjects, i.e., the constriction ratio of the glottis concerning the trachea and the curvature and shape of the airways. They employed large eddy simulation (LES) for their turbulence modeling. Farkhadnia et al. [3] investigated the geometric influence on the laminar airflow field and particle transport in G3-G6 triple bifurcations with and without partial blockage due to COPD. Johari et al. [4] compared airflow fields in computed tomography (CT)-scan generated and simplified human airway models and found that over-simplified geometries will induce noticeable differences in numerical simulation results. They also stressed that the roughness of the realistic airway walls might influence the airflow field. Xi et al. [5] investigated particle depositions in different mouth-throat models which were reconstructed via modified morphological parameters of four prototypes. They discovered that the degree of realism of the airway models significantly affected particle deposition from the oral cavity to the glottis, while the effect of the oral airway curvature was minor.
In this study, airflow patterns in three human respiratory systems were investigated using a validated computational fluid-particle dynamics (CFPD) model. Numerical simulations were performed for different steady-state inhalation flow rates (Q in = 37 and 75 L/min). The ultimate goal is to identify morphological parameters that significantly influence the airflow and coupled nanoparticle transport characteristics, and pave the way to construct a representative human upper respiratory system geometry. It is the subject-specific upper lung-airway geometry which generates the "inlet conditions" for the rest of the respiratory tract, and hence determines subsequent airflow fields and nanoparticle deposition. In this paper, we test the hypothesis that despite structural differences in human airway geometries, the basic airflow patterns (i.e., local velocity distribution, glottis jet, and turbulence intensity), as well as nanoparticle depositions, are quite similar.

Human Respiratory System Geometries
Three human upper airway geometries (Figure 1a-c), were selected and connected to a subject-specific tracheobronchial tree from the middle of the trachea to bifurcations up to Generation 8 (G8) [6]. Each human upper airway geometry contains a 20-mm inhaler inlet, oral cavity, pharynx, larynx, and a portion of the trachea. Inlet centers are all located at (x,y,z) = (0,0,0), and the upper airway geometries end in the plane at x = 0.12 m. Farkhadnia et al. [3] investigated the geometric influence on the laminar airflow field and particle transport in G3-G6 triple bifurcations with and without partial blockage due to COPD. Johari et al. [4] compared airflow fields in computed tomography (CT)-scan generated and simplified human airway models and found that over-simplified geometries will induce noticeable differences in numerical simulation results. They also stressed that the roughness of the realistic airway walls might influence the airflow field. Xi et al. [5] investigated particle depositions in different mouth-throat models which were reconstructed via modified morphological parameters of four prototypes. They discovered that the degree of realism of the airway models significantly affected particle deposition from the oral cavity to the glottis, while the effect of the oral airway curvature was minor.
In this study, airflow patterns in three human respiratory systems were investigated using a validated computational fluid-particle dynamics (CFPD) model. Numerical simulations were performed for different steady-state inhalation flow rates (Qin = 37 and 75 L/min). The ultimate goal is to identify morphological parameters that significantly influence the airflow and coupled nanoparticle transport characteristics, and pave the way to construct a representative human upper respiratory system geometry. It is the subject-specific upper lung-airway geometry which generates the "inlet conditions" for the rest of the respiratory tract, and hence determines subsequent airflow fields and nanoparticle deposition. In this paper, we test the hypothesis that despite structural differences in human airway geometries, the basic airflow patterns (i.e., local velocity distribution, glottis jet, and turbulence intensity), as well as nanoparticle depositions, are quite similar.

Human Respiratory System Geometries
Three human upper airway geometries (Figure 1a-c), were selected and connected to a subject-specific tracheobronchial tree from the middle of the trachea to bifurcations up to Generation 8 (G8) [6]. Each human upper airway geometry contains a 20-mm inhaler inlet, oral cavity, pharynx, larynx, and a portion of the trachea. Inlet centers are all located at (x,y,z) = (0,0,0), and the upper airway geometries end in the plane at x = 0.12 m.
Specifically, Geometry A and Geometry B are two in-house subject-specific human upper airway geometries, which were reconstructed from computed tomography (CT) scans (Figure 1a,b). The tracheobronchial tree used for all three configurations came from the same set of computed tomography (CT) scans of Geometry A. Geometry C, i.e., the idealized human upper airway, is shown in Figure 1c, and was created based on the geometric dimensions measured by Cheng et al. [7].  Specifically, Geometry A and Geometry B are two in-house subject-specific human upper airway geometries, which were reconstructed from computed tomography (CT) scans (Figure 1a Table 1 lists geometric parameters at different cross-sections in the three upper airways (see Figure 2a-c). Locations of different cross-sections are also shown in Section 3.2.

Mesh Generation and Independence Test
Computational meshes were generated with a multi-layer core region consisting of dense hybrid tetrahedral elements and prism layers. Prism layers were generated near the wall surface to contain the viscous sub-layers fully and to resolve any geometric features present there (Figure 3a-c).
Such high local mesh resolution is also necessary to accurately calculate near-wall derivative values. The thickness of the first prism layer guarantees y + < 1 to capture the laminar and transitional boundary layers correctly, where y + is the dimensionless wall distance (or local near-wall Reynolds number). Specifically, y + < 1 is employed to resolve the viscous sub-layer for optimum accuracy using the shear stress transport (SST) transition model [9,10]. If the y + -value is too large, then the transition onset location will not be accurately predicted and will move upstream with increasing y + . The mesh topology was determined by refining the mesh until grid independence of the flow field solutions was achieved. As an example, the final mesh of Geometry A contains 11,601,350 elements and 3,403,103 nodes. Mesh details are shown in Figure 3a,b. Considering only the different upper parts of Geometries A to C, the necessary mesh elements are regenerated and tested for the upper airways. The final meshes of B and C contain 10,465,168 and 7,043,696 elements, respectively. The worst skewness is guaranteed at values higher than 0.4. These are critical parameters which may significantly influence airflow pattern and particle deposition as well. Specifically, circularity Cr is defined as Equation (1): Perimeter of an area equivalent circle Perimeter of a luminal area (1) The hydraulic diameter D h is defined as Equation (2): Table 1 lists geometric parameters at different cross-sections in the three upper airways (see Figure 2a-c). Locations of different cross-sections are also shown in Section 3.2.

Mesh Generation and Independence Test
Computational meshes were generated with a multi-layer core region consisting of dense hybrid tetrahedral elements and prism layers. Prism layers were generated near the wall surface to contain the viscous sub-layers fully and to resolve any geometric features present there (Figure 3a-c).
Such high local mesh resolution is also necessary to accurately calculate near-wall derivative values. The thickness of the first prism layer guarantees y + < 1 to capture the laminar and transitional boundary layers correctly, where y + is the dimensionless wall distance (or local near-wall Reynolds number). Specifically, y + < 1 is employed to resolve the viscous sub-layer for optimum accuracy using the shear stress transport (SST) transition model [9,10]. If the y + -value is too large, then the transition onset location will not be accurately predicted and will move upstream with increasing y + . The mesh topology was determined by refining the mesh until grid independence of the flow field solutions was achieved. As an example, the final mesh of Geometry A contains 11,601,350 elements and 3,403,103 nodes. Mesh details are shown in Figure 3a,b. Considering only the different upper parts of Geometries A to C, the necessary mesh elements are regenerated and tested for the upper airways. The final meshes of B and C contain 10,465,168 and 7,043,696 elements, respectively. The worst skewness is guaranteed at values higher than 0.4.

Governing Equations
The Euler-Lagrange method, enhanced with in-house C programs, is employed to simulate airflow and nanoparticle transport and deposition in the human respiratory systems in future publications. The steady-state airflow in the airways was assumed to be incompressible.

Governing Equations
The Euler-Lagrange method, enhanced with in-house C programs, is employed to simulate airflow and nanoparticle transport and deposition in the human respiratory systems in future publications. The steady-state airflow in the airways was assumed to be incompressible.

Airflow Field Equations
The airflow dynamics in the respiratory tract is always unsteady and driven by the pressure differences under the action of the cyclic breathing process. In this study, the inhalation conditions are simplified to steady-state. The conservation laws of mass and momentum can be written in tensor form as Equations (3) and (4): where u j represents the fluid velocity, p is the pressure, and g j are the body forces, including gravity, electromagnetic force, etc. The viscous stress tensor τ ij in Equation (4) is given by Equation (5): With typical inhalation rates, airflow through the oral airway region and the first few generations is incipient turbulent, becoming laminar again at the fourth to sixth generation and remaining so after that. Therefore, a SST transition model was adapted for this study, predicting "laminar-to-turbulent" transition onset, and providing computational efficiency and good accuracy when compared to large eddy simulation. The use of the SST transition model was validated with three-dimensional (3-D) in vitro velocity data in a physical model of a subject-specific human respiratory system [11].

Nanoparticle Dynamics Equations
A Lagrangian frame of reference for the trajectory computations will be employed for nanoparticle dynamics. Assuming nanoparticles to be spheres, in light of the large particle-to-air density ratio and negligible thermophoretic forces, the reduced particle trajectory Equation (6) reads: in which u P i and m p are the velocity and mass of the particle, respectively; F D i represents the drag force; F G i is the gravity [12]; F BM i is the Brownian motion induced force [13][14][15]; and F L i is the Saffman lift force [16].

Numerical Setup
The numerical solution of the governing equations with appropriate boundary conditions was achieved with a user-enhanced, commercial finite-volume based program, i.e., ANSYS Fluent and CFX 18.0 (ANSYS Inc., Canonsburg, PA, USA). The SST transition model was employed to solve the laminar-to-turbulence airflow fields. Numerical simulations were performed on a local Dell Precision T7500 workstation with 40 GB RAM and 12 3.33-GHZ CPUs. The steady-state solution of the flow field was assumed to be converged when the dimensionless mass residual <10 −4 . The user-enhanced C programs, i.e., user-defined functions (UDFs) can perform the following tasks for the full Euler-Lagrange model for future nanoparticle simulations: (1) Apply anisotropic correction on turbulence fluctuation velocities. (2) Select two steady-state inlet conditions representing different scenarios of drug inhalation, i.e., 37 L/min and 75 L/min. (3) Apply uniform gauge pressure at the terminal outlets. It is worth mentioning that the outlet boundary condition is not 100% realistic due to the limitation on measuring pressures at small airways. However, this is the widely used and commonly recognized boundary condition at airway outlets, which also will not influence the subject variability studies of the upper airway morphology influence.

Results and Discussion
In this section, after validating the CFPD model, results and discussion of the numerical results focus on the following topics: (1) Mainstream and secondary flow patterns at multiple cross-sections (see Section 3.2. for details), i.e., standard deviation of the velocity. (2) Length and deviation of laryngeal jet cores from the centerlines.
(3) Local turbulence kinetic energy (TKE) distribution and averaged TKE at cross-sections AA to LL listed in Table 1.
Morphological parameters such as circularity and constriction ratio that affect airflow patterns are identified as well.

Model Validations
Due to the importance of accurately predicting transitional and turbulent flows in human upper airways, airflow field simulations with the SST transitional model had to be validated. So, the turbulent airflow results were compared with 3-D in vitro data of velocity contours in a subject-specific human respiratory system (i.e., Geometry as shown in Figure 1a), using Magnetic Resonance Velocimetry (MRV) of the identical physical model, i.e., Geometry A [11].
The numerical setup for the comparison was the same as for the experimental measurements [11]. For example, the inflow rate was 3.78 L/min of water, and the outlet condition was 20% of the total inflow for each lobe. Contours of non-dimensional local velocity magnitudes

Results and Discussion
In this section, after validating the CFPD model, results and discussion of the numerical results focus on the following topics: (1) Mainstream and secondary flow patterns at multiple cross-sections (see Section 3.2. for details), i.e., standard deviation of the velocity. (2) Length and deviation of laryngeal jet cores from the centerlines.
(3) Local turbulence kinetic energy (TKE) distribution and averaged TKE at cross-sections AA′ to LL′ listed in Table 1.
Morphological parameters such as circularity and constriction ratio that affect airflow patterns are identified as well.

Model Validations
Due to the importance of accurately predicting transitional and turbulent flows in human upper airways, airflow field simulations with the SST transitional model had to be validated. So, the turbulent airflow results were compared with 3-D in vitro data of velocity contours in a subjectspecific human respiratory system (i.e., Geometry as shown in Figure 1a), using Magnetic Resonance Velocimetry (MRV) of the identical physical model, i.e., Geometry A [11].
The numerical setup for the comparison was the same as for the experimental measurements [11]. For example, the inflow rate was 3.78 L/min of water, and the outlet condition was 20% of the total inflow for each lobe. Contours of non-dimensional local velocity magnitudes ‖ ‖ at the sagittal plane y = 0 (Figure 4a) are shown in Figure 4a-c. Additionally, iso-surfaces, visualizing the epiglottal and laryngeal jet cores (‖ ‖ = 1.5), are also compared in Figure 5a,b. Good agreements can be observed in both cases.

Secondary Flow Patterns in Cross-Sections and Sagittal Planes
Local velocity distributions and secondary flow vectors at different cross-sections and sagittal planes (see Table 1) in the three geometries for Q in = 37 L/min are shown in

Secondary Flow Patterns in Cross-Sections and Sagittal Planes
Local velocity distributions and secondary flow vectors at different cross-sections and sagittal planes (see Table 1   the contour in the sagittal plane in Figure 7). Recirculation is generated at the upper trachea. However, no such flow patterns are observed in Geometries A and C. 3. Although the lung airway trees employed in the three geometries are the same, mass flow distributions to the left and right lobes are different because of geometric differences among upper airways. For example, Dean's flows can be observed at KK′ and LL′. However, the local mass flow rate distributions are different, which will potentially influence the fraction of nanoparticles moving into different lung branches.

2.
A laryngeal jet core is generated in each geometry at the glottis associated with recirculation zones in the trachea.
Distinguished local airflow behavior can be summarized as follows: 1.
There are two counter-rotating vortices following the mainstream at both sides in the oral cavities in Geometries A and B, due to the non-circular geometry effect (see AA and BB in Figures 6-8). However, no similar structures exist in the idealized upper airway (Geometry C) with perfect circular outlines at those cross-sections.

2.
Due to the large constriction ratio and low circularity at the glottis in Geometry B, high flow resistance around the glottis region breaks up the jet core and spreads it towards the wall (see the contour in the sagittal plane in Figure 7). Recirculation is generated at the upper trachea. However, no such flow patterns are observed in Geometries A and C.

3.
Although the lung airway trees employed in the three geometries are the same, mass flow distributions to the left and right lobes are different because of geometric differences among upper airways. For example, Dean's flows can be observed at KK and LL . However, the local mass flow rate distributions are different, which will potentially influence the fraction of nanoparticles moving into different lung branches.
In summary, Geometry C is not able to quantitatively maintain the complex secondary flow patterns in subject-specific geometries. Nevertheless, it is still a computationally efficient geometry to qualitatively mimic the geometric characteristics and hence the resulting major flow patterns. In summary, Geometry C is not able to quantitatively maintain the complex secondary flow patterns in subject-specific geometries. Nevertheless, it is still a computationally efficient geometry to qualitatively mimic the geometric characteristics and hence the resulting major flow patterns.

Laryngeal Jet Core Structures
To better visualize the laryngeal core structures for = ‖ ‖/ = 1.5, Figure 9a-c show the

Laryngeal Jet Core Structures
To better visualize the laryngeal core structures forÛ = → u /U in = 1.5, Figure 9a-c show the iso-surfaces in three geometries for two steady-state inhalation flow rates, i.e., Q in = 37 L/min and Q in = 75 L/min. Specifically, U in is the average inlet velocity at the mouth.
In contrast to the conclusion drawn by Choi et al. [2], i.e., "at a higher flow rate, the jet core is found to be shorter and detaches from the tracheal wall", no decrease in jet-core length was detected between Q in = 37 L/min and Q in = 75 L/min in Geometries B and C (see Figure 10). However, in Geometry A (see Figure 9a), jet core length reduction was found, i.e., higher flow rate induced a shorter length of the laryngeal jet core due to the energy and momentum dissipations when the flow is passing the complex laryngeal region. Specifically, higher inhalation flow rates will lead to stronger impaction, which may or may not induce the redistribution of the jet core region. Additionally, jet-core impaction and breakup depend not only on the inhalation flow rate but also on the jet core orientation with the trachea centerline. impaction, which may or may not induce the redistribution of the jet core region. Additionally, jetcore impaction and breakup depend not only on the inhalation flow rate but also on the jet core orientation with the trachea centerline. Indeed, Figures 6-10 indicate that laryngeal jet cores deviate from the trachea centerlines with highly distinguished patterns in all three geometries for different inhalation flow rates. For example, the jet core impacts the front right trachea in Geometry C, the front left trachea in Geometry B, and the back left trachea in Geometry A, respectively. No common characteristics were discovered, due to the highly different geometric characteristics among the three geometries. Furthermore, Figure 10 indicates that the constriction ratio has an impact on the length of the jet core. Specifically, a small constriction ratio generates a stronger jet core with a longer core length (Geometry C vs. Geometry A).  Because nanoparticles follow the airflow streams in human respiratory systems, it is reasonable to expect that the impingement location and the length of the laryngeal jet core can significantly influence the local deposition of nanoparticles. By guiding patients to operate a drug inhaler with different inhalation flow rates, local deposition of pulmonary nano drugs can be improved for treating diseases in different lung regions. impaction, which may or may not induce the redistribution of the jet core region. Additionally, jetcore impaction and breakup depend not only on the inhalation flow rate but also on the jet core orientation with the trachea centerline. Indeed, Figures 6-10 indicate that laryngeal jet cores deviate from the trachea centerlines with highly distinguished patterns in all three geometries for different inhalation flow rates. For example, the jet core impacts the front right trachea in Geometry C, the front left trachea in Geometry B, and the back left trachea in Geometry A, respectively. No common characteristics were discovered, due to the highly different geometric characteristics among the three geometries. Furthermore, Figure 10 indicates that the constriction ratio has an impact on the length of the jet core. Specifically, a small constriction ratio generates a stronger jet core with a longer core length (Geometry C vs. Geometry A).  Because nanoparticles follow the airflow streams in human respiratory systems, it is reasonable to expect that the impingement location and the length of the laryngeal jet core can significantly influence the local deposition of nanoparticles. By guiding patients to operate a drug inhaler with different inhalation flow rates, local deposition of pulmonary nano drugs can be improved for treating diseases in different lung regions. Indeed, Figures 6-10 indicate that laryngeal jet cores deviate from the trachea centerlines with highly distinguished patterns in all three geometries for different inhalation flow rates. For example, the jet core impacts the front right trachea in Geometry C, the front left trachea in Geometry B, and the back left trachea in Geometry A, respectively. No common characteristics were discovered, due to the highly different geometric characteristics among the three geometries. Furthermore, Figure 10 indicates that the constriction ratio has an impact on the length of the jet core. Specifically, a small constriction ratio generates a stronger jet core with a longer core length (Geometry C vs. Geometry A).
Because nanoparticles follow the airflow streams in human respiratory systems, it is reasonable to expect that the impingement location and the length of the laryngeal jet core can significantly influence the local deposition of nanoparticles. By guiding patients to operate a drug inhaler with different inhalation flow rates, local deposition of pulmonary nano drugs can be improved for treating diseases in different lung regions.
Another interesting observation is that, despite the axisymmetric shape of the idealized upper airway of Geometry C, the laryngeal cores shown in Figures 8 and 9c are not axisymmetric. This is due to the influence of the reverse flow in the trachea, which conveys the non-axisymmetric geometric effect of the lower lung airways towards the upstream region.
Since the conclusions are diversified between our studies and existing papers [2], more subjects should be included in the future to identify geometric and operational parameters that significantly influence the laryngeal jet core structure. Another interesting observation is that, despite the axisymmetric shape of the idealized upper airway of Geometry C, the laryngeal cores shown in Figures 8 and 9c are not axisymmetric. This is due to the influence of the reverse flow in the trachea, which conveys the non-axisymmetric geometric effect of the lower lung airways towards the upstream region.

Turbulence Kinetic Energy (TKE)
Since the conclusions are diversified between our studies and existing papers [2], more subjects should be included in the future to identify geometric and operational parameters that significantly influence the laryngeal jet core structure.  Additionally, Figure 12 shows the averaged TKE of cross-sections AA′ to JJ′ (see Figures 6-8 for cross-section locations). The averaged TKE abruptly increases after the constriction of the glottis. It can be observed that the constriction ratio has a significant impact on the higher TKE in the trachea. Indeed, Geometries B and C have similar trends of TKE after the airflows pass the glottis contractions (see Figure 12). In contrast, with a higher constriction ratio, the increase of TKE in Geometry A after the glottis is less steep than that in the other two. The comparison between Geometries B and C also Additionally, Figure 12 shows the averaged TKE of cross-sections AA to JJ (see  for cross-section locations). The averaged TKE abruptly increases after the constriction of the glottis. It can be observed that the constriction ratio has a significant impact on the higher TKE in the trachea. Indeed, Geometries B and C have similar trends of TKE after the airflows pass the glottis contractions (see Figure 12). In contrast, with a higher constriction ratio, the increase of TKE in Geometry A after the glottis is less steep than that in the other two. The comparison between Geometries B and C also implies that the glottis contraction serves as a transitional flow regulator due to the forced mixing effect when airflow passes the contraction region. Although the averaged TKE at cross-sections differs in the two geometries, their downstream-glottis trends are similar. However, more geometries need to be included in the inter-subject variability study to test this "glottis regulator" concept.

Conclusions
In this study, numerical inter-subject variability studies were performed for airflows, considering drug inhalation flow rates. It can be concluded that the upper airway morphology has a significant influence on actual flow patterns. This will potentially lead to different deposition patterns of nanoparticles. It was also found that the jet core length does not increase with the increase of inhalation flow rates in all upper airways.
The idealized human upper airway geometry (Geometry C) may not be able to fully represent the naturally asymmetrical flow patterns in subject-specific human respiratory systems. However, Geometry C can be used for initial studies because it qualitatively represents basic flow patterns in human respiratory systems. However, the use of Geometry C may lead to lower local nanoparticle depositions when compared to subject-specific lung-airways.

Future Work
Considering the complexity of multiple geometric parameters, it is necessary to include additional airway geometries representing different population groups to generate statistically robust conclusions, i.e., CFPD simulation results with "error bars." The proposed framework of in silico subject variability study is shown in Figure 13, to pave the way for the development of a virtual human system as the individualized digital twin for personalized pulmonary disease treatment planning.

Conclusions
In this study, numerical inter-subject variability studies were performed for airflows, considering drug inhalation flow rates. It can be concluded that the upper airway morphology has a significant influence on actual flow patterns. This will potentially lead to different deposition patterns of nanoparticles. It was also found that the jet core length does not increase with the increase of inhalation flow rates in all upper airways.
The idealized human upper airway geometry (Geometry C) may not be able to fully represent the naturally asymmetrical flow patterns in subject-specific human respiratory systems. However, Geometry C can be used for initial studies because it qualitatively represents basic flow patterns in human respiratory systems. However, the use of Geometry C may lead to lower local nanoparticle depositions when compared to subject-specific lung-airways.

Future Work
Considering the complexity of multiple geometric parameters, it is necessary to include additional airway geometries representing different population groups to generate statistically robust conclusions, i.e., CFPD simulation results with "error bars." The proposed framework of in silico subject variability study is shown in Figure 13, to pave the way for the development of a virtual human system as the individualized digital twin for personalized pulmonary disease treatment planning.