Particulate Deposition in a Patient With Tracheal Stenosis

The rate of respiration, airway geometry, particle sizes, and the presence or absence of airway disease are all factors that significantly influence the accumulation of aerosol or particles in the respiratory tract of humans; as a result, these factors also affect the therapeutic efficacy of medications that are inhaled. Airway diseases may modify deposition locations, which eventually impacts the effectiveness of inhaled drugs. Other factors including the thickness of the mucus layer, humidity, and mechanisms of mucus clearance could also impact the absorption rate and deposition of the inhaled drug. In recent decades, the effectiveness of inhaled particulate and aerosolized medications as a drug administration route has been explored extensively. Among the most critical components affecting the deposition of particles or inhaled medications include the tracheobronchial tree geometry, the respiration rate, and the particle diameter of the aerosolized medication. Airway obstruction results in a reduction of airway lumen, which consequently alters the flow characteristics and deposition pattern of the particles that are inhaled [1].

Ruffin et al. [2–4] explored the deposition of histamine aerosol in the human respiratory system and deduced that the presence of increasing airway blockages, depositions were more efficient in larger airways as compared to the distal airways distributed throughout the lung. This result implied that the larger airways represent the primary location of histamine receptors.

Labiris and Dolovich [1] have demonstrated that the drug delivery efficiency and consequent efficacy of therapy may decline if the pharmaceutical aerosols are directed to an unintended region of the airways. In patients with asthma, inflammatory cells such as lymphocytes are homogeneously dispersed throughout both small and large airways [5]; therefore, they proposed that drug delivery efficiency might increase and be more valuable if the inhaled anti-inflammatory medication could be dispersed to the entire tracheobronchial tree.

In the present investigation, numerical investigation of the effect of tracheal stenosis on air flow characteristics and particle depositions is presented. Tracheal stenosis, the most commonly occurring injury in the trachea, is often triggered by intubation. Tracheal intubation is a procedure where a flexible plastic tube or stent is placed into the trachea. Granulation formation resulting from the steady deterioration of cartilage could occur at either end of the stent [6].

Airway stenosis results in increased flow resistance due to obstruction and thus augments breathing problems especially in patients with pneumonia and collagen vascular diseases. It has resulted in increased death occurrences [7,8].

Computational fluid dynamics (CFD) has been used by researchers to evaluate airflow characteristic and particle depositions in human airways. Many groups including Brouns et al. [8], Zhang and Kleinstreuer [9], and Luo et al. [10,11] have explored tracheal stenosis in the context of the idealized/smooth shaped model. Since geometry possesses a significant effect on particle deposition and aerosolized medication, this makes simulations of these events using computed tomography (CT)-based models necessary. Investigation of flow characteristics of a realistic three-dimensional (3D) model including the lower generations of the airways (downstream of the tracheal stenosis) is needed to identify aerosol/particle deposition and transport patterns in both the trachea and branching sections. The inclusion of the branching airways was not present in several previous investigations [6,8,9,12,13], with few exceptions [10,14]. Luo et al. [10] examined the obstruction of airways in different generations of a Weibel-based (idealized) model, while Xi et al. [14] examined the effects of laryngeal jet inclusion on the behavior of inhaled aerosols deposition in an approximated model up to the sixth generation (G6) branching of the airways. It is concluded that the most effective position for the deposition of aerosolized antibiotics or therapeutic particulates is greatly dependent on the infection type or disease being treated.

It is vital to be familiar with the particle depositions' efficiencies in drug delivery, whether they are toxic particles or therapeutic aerosols. The location and distribution uniformity of the deposited particles play major roles in the efficacy of the therapy. The details of the aforementioned processes are yet to be fully elucidated, especially in patients with diseased airways and are thus the concentration of the current study.

Multidimensional numerical models including a 3D CFD coupled with zero-dimensional [15,16] or one-dimensional (1D) [17–22] models can be used to simulate the deeper generations together with the larger 3D airways. Tawahi et al. [19,22] have modeled the full conducting airway system by using a volume-filling algorithm. This approach is used for 3D and 1D boundary coupling and has also been utilized with the CT-based geometry [21,22].

Malve et al. [23] simulated the human trachea using fluid–structural interaction under poststent, normal breathing, and coughing conditions; airway branches and different generations were not included in this study. In subsequent investigations [24,25], they have added a 1D outflow boundary conditions model to their 3D tracheal model to predict the impedance and pressure wave forms. The 3D–1D impedance boundary condition was simulated for evaluation of tracheal stenosis [25], including a single bifurcation before and after prosthesis implantation [23]. Their results showed that tracheal muscle deflection is prevented by the existence of the stent. Malve et al. [23] indicated that the presence of local recirculatory flow might play a role in mucus accumulation on the stent's top surface.

Air movement during inspiration is the result of downward movement of the diaphragm, which increases the volume of the thoracic cavity, thus lowering the air pressure in the lungs, in comparison to the atmospheric pressure. Modeling a pressure-driven flow is one of the challenges in defining the boundary conditions for CFD simulations of human airways, since in vivo measurements are not practically possible in all generations of the bronchial tree.

Yang et al. [11] designated a constant inlet mass flow rate in their idealized model with steady pressure at the outlets of contracted airways in different branches. They demonstrated that the airway constrictions result in a jet-like behavior, and at some points, the amplification of flow momentum, and the skewed velocity profile resulting from the jet-like profile, does not possess a sufficient length to recover prior to the oncoming branching section. In the work of Xi et al. [14], the inlet velocity and mass flow distribution rate were used for boundary conditions.

Mylavarapu et al. [13] designated pressure boundary conditions at the inlet and outlet to make the system pressure driven. Although their model included the Larynx and Oropharynx areas, they did not include different branching sections after the trachea.

Chen et al. [12] defined the boundary conditions as the widely employed zero pressure outlet and an assigned velocity inlet. Although these boundary conditions may be suitable in the healthy trachea and when lower regions are not included, it may not properly characterize the patient with an obstruction that causes a considerable pressure drop in the system, or models in which the lower generations of the airways are integrated.

Walters et al. [26] proposed a simple approach for simulating flow in the terminal and truncated outlets of a human airway, based on steady-state simulation with stochastically coupled pressure boundary conditions. The method uses a prescribed time-dependent flow during the breath cycle as input, and it does not require prior prescribed distal mass flow or pressure distribution. They compared their results with corresponding results for fully resolved airway geometries where they concluded that reduced geometry with appropriate boundary conditions could be an alternative to fully 3D simulations, throughout all 16 generations of the conducting zone.

In this investigation, CT-based model along with patient-specific outlet boundary conditions is used in order to make the flow pressure driven. The effect of tracheal stenosis pre- and postintervention on particle deposition in the lower generations of the airways is presented.

For the current investigation, the model (Fig. 1.) used is a CT-based model of the airways from a patient with tracheal stenosis (the patient has a blockage only on the lateral sides of the trachea with no posterior and anterior changes when compared to a healthy subject without the flow limiting segment). De-identified CT scan images consisting of 127 slices in the axial position, with a slice thickness of 3 mm and pixel size of 0.976 mm, were imported into the 3D-Doctor software, and the 3D models were generated using a variety of segmentation methods [27]. The model was subsequently imported into starccm+ for CFD analysis (Fig. 1). The final model includes up to the eight generations of the airways. Comparisons of flow characteristics, particle deposition, and pressure variation before and after the virtual intervention, when the stenosis is removed, were presented. The reduction in the size of the tracheal lumen due to stenosis was manually modified to the average tracheal diameter, representing the patient's ideal postsurgical outcome, without the flow-limiting segment. Inhalation is an effective method of drug delivery; however, precise quantities of the therapeutic medication must be deposited after the oropharyngeal region [1]; the models in this investigation include the branching airways downstream of this area.

Here, P = Pressure, t = time, and BP_max = maximum branch pressure found from the initial step.

As a result of the structure of the throat, oral cavity, or various obstructions, the flow may be laminar, becomes transitional, and subsequently turbulent. Thus, accurate turbulence modeling is necessary to describe various flow regimes [29]. The critical Reynolds number for steady flow through a straight pipe with smooth inner walls and a circular cross section is commonly used to approximate the flow pattern between these conditions. However, care needs to be employed when predicting flow conditions based solely on the Reynolds number, due to the fact that obstructions, skewed airways, and cartilaginous rings can trigger turbulence in larger airways [30]; thus, turbulence may occur in the trachea despite the fact that the critical Reynolds number has not been reached [31].

Implicit unsteady simulation for 5 s of respiration was implemented using sinusoidal-type flow behavior during normal tidal volume (500 ml) for a single breath. The flow governing equations are discretized, using second-order finite volume schemes, on the computational domain and a second-order implicit discretization was used for the time integration.

Equations (1) and (2) terms: ρ is the density, u are the velocity components, k is the kinetic energy of turbulence, x is the Cartesian coordinate, P is the production of kinetic energy of turbulence, σ and β are the turbulence-model coefficients, β*and γ are turbulence-model coefficients that are modified for with low Reynolds numbers, ω is the specific dissipation of turbulence kinetic energy, μ is the dynamic viscosity, and F₁ is the auxiliary function in the turbulence model.

The particles that are transported in the unsteady flow simulation are injected uniformly at the inlet surface with the equivalent inlet air velocity. One-way coupling for air and particles was assumed to emphasize the fluid forces on various sized particle dispersion and deposition [35]. The Lagrangian phase model was employed to replicate transport and depositions of 1, 2.5, 4, 6, 8, and 10 μm size. The selection of sizes was determined by the Environmental Protection Agency testing standards of PM2.5 and PM10 [36]. The selection of 1 and 2.5 was chosen to determine the difference in fine particle deposition (PM2.5). The coarse particle range of PM2.5 to PM10 has a significant variation in size; therefore 4, 6, and 8 were selected. Solid round particles with a density of 1200 kg/m³ were utilized with an injection rate of 1.3 g/mm³/s for all particles throughout the inhalation process. Other studies have demonstrated the importance of modeling the Brownian motion [37,38] for nanoparticles. However, the present investigation does not include modeling of the Brownian motion due to the dominance of the flow inertial on micron particle depositions in large and medium airways.

The particle transport is governed by Newton's second law of motion where the acting forces are the drag and gravity forces and the near wall stresses. The following models [39] are chosen for capturing particle transport and deposition in the human respiratory system:

Here, $mp$ is the mass of the particle, $vp$ is the velocity of the particle, $Fg$ is the gravity force, $Fp$ is the pressure gradient force, and $FD$ is the drag force acting on the particle. One-way coupling is assumed, meaning that the continuous phase influences the particles not the other way around. The effect of droplet evaporation and secondary breakups was not taken into consideration.

where $∇Pstatic$ is the static pressure gradient in the continuous phase, and $Vp$ is the particle volume.

Normally, the particles do not rebound in the actual human trachea as a result of the airway's mucus layer. Therefore, a stick wall boundary condition was utilized. Both inspiratory and expiratory processes were analyzed, due to the back pressures produced by exhalation, which can cause deposition before the first generation, downstream of the tracheal stenosis. The importance of analyzing both inspiration and expiration as opposed to just inspiration is that the majority of particles are deposited during the inspiration cycle; however, suspended particles after inspiration may be deposited in the airways or expelled from the body during expiration, as shown experimentally by Usmani et al. [40].

A mesh independency test was executed with 0.6, 1.4, 2.5, and 6 × 10⁶ mesh cells, with particle depositions (for 10 μm) of 48.2%, 46.9%, 47%, and 47%, respectively. 1.4–2.4 ×10⁶ mesh cells have been demonstrated [28] to be adequate with respect to the deposit of particles. Over 2.4 × 10⁶ prismatic and polyhedral layer meshes added near the wall with a y⁺ < 1 were chosen for the present simulations.

The CFD pressure data validation was completed by comparing the corresponding results for a rigid model of the identical airway up to three generations. The model was generated via utilization of additive layer manufacturing. Figure 2(a) illustrates the experimental setup while Fig. 2(b) is a schematic of the experiment. The outlet of each branch was linked to a pressure transducer with T-branch tubes, and each branch was fastened to a solenoid valve in which the flow volume could be controlled by a variety of valve gate openings. For each lobe of the lung, the valve outlets were attached to a manifold; the manifold was then attached to a piston-cylinder system using a computer-controlled motor, where a multitude flow conditions could be reproduced. The experiments were conducted with a sinusoidal flow (unsteady), which represented a natural breathing cycle of inhalation and exhalation. The experimental and numerical results for the outlet pressure were compared to validate flow characteristics of the simulation.

Figure 3 shows contours of the mean velocity at the sagittal plane during the inhalation process for pre- and post-intervention. Figure 4 shows the mean velocity and the vorticity magnitudes during the inhalation at this plane for different time-steps. As the air flows through the obstructed region, there is flow acceleration with an area of recirculation detected downstream from the stenosis. However, regions of flow separation and recirculation are greatly decreased following the interventions. Wilcox [41] showed that an obstruction resulting from the stenosis would cause flow separation, which results in a transition of flow from laminar state to a turbulent state and finally reattachment. Zang and Kleinstreuer [9] detected that the flow is mostly laminar near the intake; however, it could potentially become turbulent as a result of the stenosis and then relaminarize in subsequent generations of the airways as there is a reduction in the Reynolds number. Based on the current results, downstream from the stenosis, there are areas of large recirculation and Coanda-type wall re-attachment. The instability of the jet, a characteristic generally denoted as the Coanda effect, is the consequence of large region of recirculation, which induces a strong Laryngeal jet, in addition to the secondary flows and very skewed velocity profiles both downstream and upstream of the stenosis during expiration phase [9,14].

The behavior of the air flow is very different throughout the exhalation cycle. As exhalation occurs and air flows from different airway branches merge together, higher levels of turbulence are detected as compared to the flow conditions observed during the inhalation following the intervention.

A drop in static pressure, as well as a maximal wall shear stress, is detected proximally to the stenosis in the trachea [7,13]. Figure 5 displays the changes in wall shear stress pre- and post-interventions at the peak of inhalation. The region of stenosis has the greatest value of wall shear stress in the pretreatment trachea. Analyzing immediately downstream of the stenosis, the area of increased velocity and recirculation exhibits a large wall shear stress in comparison to the same area post-intervention (which has a dramatic reduction in wall shear stress). In the branch areas for pre- and post-intervention, the wall shear stress values are unchanged.

The Coanda flow has resulted in increased turbulence at the stenosis and immediately downstream. Aerosol/particle deposition is generally augmented as a result of increased turbulence in the respiratory tract [14].

The air flow is unsteady due to the complex airway structure in the upper airways. The presence of an obstruction alters the distribution of pressure, which increases airway resistance and may give rise to breathing difficulties. Figures 6 and 7 display the drop in pressure within the left and right branches at various breathing volumes for pre- and post-interentions. Significant pressure drops are detected upstream of the stenosis at a location where the effects are amplified with increased air volume. Drops in pressure can also be observed in other regions, further downstream at the lower generations, in both pre- and post-tracheal stenosis intervention.

The experimental setup for validation of boundary conditions is shown in Fig. 2. Figure 8 displays the relationship of the normalized pressure results from the CFD simulation with the two cases of inspiratory–expiratory flow (sinusoidal air flow). The data of the pressure for every outlet/lobe are presented at the third generation. Based on the primary results and careful examination of the airway model, it was deduced that the surface roughness influenced the pressure differences. In the CFD model, the surface was assumed to be smooth while for the experimental model, there was a significant roughness throughout the contour of the interior walls due to the quality of the 3D printed material. The model was treated by an acetone bath to remove these imperfections and to smooth the surface. The comparison of the initial and cleaned model for sinusoidal flow shows that five branches experienced 5–98% improvement, three branches experienced 34–42% decline, and one branch is unchanged. An ongoing study is focused on experimenting with a variety of different methods and materials, to generate a smoother surface, between correlation between CFD and experimental models.

Stokes number represents the influence of inertial effects during particle's trajectory [29], and can describe the probability of particle deposition via impaction. Results show as the size of the particle increases, the probability of particle deposition also increases, and as the Stokes number increases, the particle deposition due to inertial impact increases. In comparison of surface deposition between a stenosis and a post-intervention trachea, the surface deposition is relatively higher at corresponding Stokes number in a stenosis trachea.

Figure 10 shows the distribution of particles within the flow domain. Results show increased particle depositions in a patient with stenosis, with solid 10 μm particles (PM-10) having a higher deposition than 2.5 μm particles (PM-2.5). PM-10 deposition pattern downstream is locally concentrated after the stenosis while in post-intervention, it is uniformly dispersed in the trachea. PM 2.5 has a high local concentration at the stenosis with no impact in post-intervention at the same location. Flow reversal in between inspiration–expiration phases may influence particle deposition which requires further investigations.

Different airway generations are predicted to experience different levels of particle depositions, which is critical for developing an effective patient-specific drug delivery method. Figures 11 and 12 demonstrate the percentage of total deposition across different generations, at different air volumes. Labiris and Dolovich [1] observed that aerosolized medication is typically deposited more central as a result of the inertial impact, and that it is deposited in a more uniform pattern in a healthy subject; however, it is also greatly dependent on other factors such as the rate of respiration and the size of the particle [27,28]. The majority of smaller particles, such as the 2.5 μm particles, are deposited in the deeper generations of the downstream airways up to the eighth generation, and tracheal stenosis does not have a crucial impact on the pattern of deposition. Particles that are deposited in the larger airways are mostly deposited between generations 4 and 6.

Large particles (PM10) are mostly deposited in the upper region of the respiratory tract, and the tracheal stenosis has a considerable effect on their deposition pattern. In the state of stenosis, most of the particles are deposited downstream of the stenosis in the G0 (trachea)–G1 (first Generation) and G3–G5; however, following removal of the stenosis, the deposition is more pronounced in G3–G5. Overall, most of the particles tend to reside in G4 and G5 for every case.

Drug delivery optimization varies for different diseases. Labiris and Dolovich and Ruffin et al. [1–4] emphasized the necessity of this subject. Obstruction of the airway in lower generations can redirect the inhaled air to the unblocked/healthy areas, which significantly influences the position where the particles/aerosols are deposited [1]; however, the characteristics of the air flow with tracheal obstructions are different.

Particle deposition during inhalation commonly results from inertial impaction as the inhaled air splits into diverse branches. Particle deposition and transport are also affected by secondary flows. This behavior is enhanced in the event of stenosis. In particular, regions of recirculation below the stenosis would enhance particle depositions and particle residence time. During exhalation, particle deposition behavior is similar to inhalation. However, particles also sedimented below the stenosis (main bronchi) due to the blockage effect and increased pressure during exhalation. The secondary flows, velocity contours of axial plane, and their effects on particle depositions were demonstrated in authors' previous investigation [28].

Computational fluid dynamics-based “virtual surgery” [42] may also be used as an additional method for planning and strategizing surgeries. There are uncertainties in CFD simulations of biological systems, which can be categorized into modeling, numerical, experimental, and biological behavior uncertainties. Modeling uncertainties consist of automated 3D geometry conversion from CT scan, turbulence models, and boundary conditions. In this study, manual segmentation and semi-automatic conversions reduce the errors associated with the conversion of CT images to 3D models. Comparison of appropriate turbulence modeling is essential for realistic flow behavior and has been studied by other investigators [29–31]. Two-step boundary conditions mimic the realistic patient-specific flow behavior to each lobe of the lung. Although the pharyngeal regions may alter particle deposition patterns in the trachea for a healthy subject, the object of this study was to emphasize the effect of pre- and post-tracheal stenosis intervention on particle depositions downstream in lower generation of the airways. The effect of pharyngeal on total depositions for smaller particles is minimal in the trachea, but significant for larger particles. The ratio of particle deposition in the lower generations for 3D models that exclude nasal cavity and pharyngeal regions compared to a complete upper airway model is the authors' current ongoing investigations [43]. Mesh refinement and grid dependency analysis were performed to reduce discretization errors. Experimental measurement uncertainties such as fabrication, material removal, and surface roughness contributed to differences in flow patterns between the numerical and experimental analyses. It should be noted that wall compliance of the airway might play a significant role in predicting accurate flow characteristics, thus requiring further investigation.

Computational fluid dynamics analysis was performed using a patient-specific trachea with stenosis, and following the removal of the stenosis (post-intervention condition). The study also investigated the consequence of stenosis on particle deposition, to determine the method by which inhaled particles transport and deposit in the presence and absence of stenosis. Investigations were performed for particles of 1, 2.5, 4, 6, 8, and 10 μm. Results indicate the existence of the stenosis inflicts a substantial pressure force on the flow of air that is inhaled, resulting in major particles deposition downstream and upstream of the stenosis. Downstream of the stenosis, the flow is initially turbulent, and there are regions of secondary flow with significant flow separation, which has resulted in higher particles deposition.

The CFD results were validated using a 3D printed airway model up to the third generation. The experimental investigation outcomes were presented for two cases of constant and sinusoidal flows. Relatively good agreements were obtained between the CFD and experimental pressure results; the large variations found at selective lobes were ascribed to flaws in the experimental 3D printed model. Experimental correlation was substantially improved after the model was treated by an acetone bath to remove imperfections and to smooth the surface.

The CT images were provided by Long Beach Veterans Administration Hospital (LBVA), pulmonary division, for which the authors are deeply appreciative.

$Ap$ =

particle projected area

BP_max =

maximum branch pressure found from the initial step

$Cd$ =

drag coefficient of the particle

D =

characteristic length scale

$Dp,dp$ =

particle diameter

$FD$ =

drag force acting on the particle

$Fg$ =

gravity force

$Fp$ =

pressure gradient force

F₁ =

auxiliary functions in turbulence model

k =

kinetic energy of turbulence

$mp$ =

mass of the particle

P =

pressure

P_k =

production of kinetic energy of turbulence

Re =

Reynolds number

St =

Stokes number

t =

time

U =

peak velocity of inhalation

u_j =

the velocity components

$vp$ =

velocity of the particle

$vS$ =

particle slip velocity

$Vp$ =

particle volume

$xj$ =

Cartesian Coordinates

β and σ =

turbulence-model coefficients

β^*and γ =

turbulence-model coefficients that are modified for with low Reynolds numbers

μ =

dynamic viscosity

ρ =

density

$τv$ =

momentum relaxation time-scale

ω =

specific dissipation of turbulence kinetic energy

$∇Pstatic$ =

continuous phase static pressure gradient