How microplastics are transported and deposited in realistic upper airways?

Microplastics are tiny plastic particles and have an impact on public health and environmental hazards, which is a global concern. Tiny microplastic debris is generated from the degradation of plastic products, consumer products, tire wear, and industrial breakdown. Millions of tons of these microplastic particles have been found in seas, atmospheric air, and soil.¹ Global microplastic production is surging, and the density of microplastics in the air is increasing significantly. It is evident that the tiny respirable microplastics can transport in atmospheric air. For the first time, in 2022, studies found microplastics in human airways, which raises the concern of serious respiratory health hazards.² Microplastics usually contain toxic pollutants and chemicals. The toxic substances from the inhaled microplastics could potentially cause serious health hazards in the respiratory system. Research shows humans might inhale about 16.2 bits of plastic every hour, which is equivalent to plastics used to make a credit card in an entire week.³ To date, a range of studies have analyzed the movement and behavior of microplastics in atmospheric air.^4,5 However, the transport behavior, migration, and toxicity of microplastics in human airways are never conducted (Fig. 1).

In the last few years, scientists have found microplastics in drinking water.⁶ The World Health Organization has called for more attention to microplastics.⁷ Most previous studies consider microplastics to be ingested by humans orally through water or food. However, inhalation also allows microplastics to reach human tissues.⁸ Recently, microplastic particles have been tracked in the lower airways for the first time in history.^9,10 The presence of microplastics in lower airways increases the long-term microplastic exposure concerns for respiratory health. Microplastics are tiny enough to enter the bronchioles and trap in different regions of the airway branches. The accumulation of deposited microplastics in the airways may adversely impact the human respiratory system. Computational fluid dynamics (CFD) has been employed to investigate the pollutants particle behavior in the human airway for decades.^11–15 Compared with in-vivo or in-vitro experiments, CFD simulation has proven to be an efficient method to predict the trajectories of inhaled particles without risk to human beings and help physicians learn the impact of deposited particles. There are several techniques that have been used to create the human lung airways for CFD simulations. These included the 3D airway model based on Weibel's dimensions and the image-based airway models such as CT and MR images. The CT and MR imaging technologies have been widely used for the CFD simulations to investigate the airflow characteristic, particle transport, and deposition due to their advantages that provide the realistic models under the asymmetric branching patterns.¹⁶ Some studies focused on the numerical simulation for the particle transport and deposition under various sections of the human respiratory system, including the nasal cavity through the upper airways,¹⁷ oral airway through the upper airways,^18–20 and both nasal cavity and oral airway through the extrathoracic airways.²¹ Some studies also focused on the specific region for the different purposes, such as the nasal flow therapy,²² drug delivery for human maxillary sinus,^23–25 nasal irrigation delivery for the patients after the functional endoscopic sinus surgery,²⁶ and aerosol transport through stenosis upper airway.^27,28 However, to date, accurate knowledge of microplastic transport in airways is lacking in the literature.

This study aims to analyze the microplastic transport dynamics in extrathoracic airways for the first time. A CFD-DPM model is used to analyze the fluid and particle motion dynamics in airways.

This study employed 30-year male DiCom images and constructed the three-dimensional upper airway model. The upper airway consists of nasal cavities, mouth-throat, pharynx, epiglottis, and trachea. The healthy airway model is used to analyze microplastic transport. The upper airway model and selected locations for flow analysis are presented in Fig. 2. Figure 3(a) presents the front view of the airway model, and Fig. 3(b) shows the unstructured mesh elements. The nasal cavities and mouth-throat anatomy are highly complex and contain a strong change of curvatures with uneven surfaces. Therefore, unstructured tetrahedral elements are generated, and ANSYS meshing is used for the grid generation. The grid independency test for the lung model from the present study was performed from the previous study.³¹ Inflation layers are used for the airway wall, and a proper grid refinement is performed for the present model. The final model comprises 4.3 × 10⁶ elements.

The numerical model is validated with the literature^32,33 (Fig. 4). The overall deposition rate is calculated against the impaction parameters and compared with the published literature. The deposition efficiency (DE) is compared against of impaction parameter [(d_p² (μm²) Q (cm³/s)] for the nasal cavity. The present model's finding shows a similar trend in the deposition rate with the available literature.

Auxiliary planes and lines have been created to better represent the fluid flow characteristics at different parts of the model. The first plane is a coronal plane at the nasal cavities, and the last plane is a transverse plane located at the end of the trachea. In between these two planes, there are other lines and planes that encompass the whole model, including the larynx, pharynx, and the main trachea itself, as indicated in Fig. 2 (see for the location of selected lines and planes).

Velocity contours and vectors on different planes for the flow rates of 7.5 and 30 l/min are depicted in Fig. 5. The comparison of velocity on different planes indicates that the local distortions in the geometry of airways substantially affect momentum distribution and change the velocity contour patterns. At planes 1 and 2, the velocity is highly affected by local complications of the airways and is unevenly distributed on these planes. At plane 1, the velocity at the center of the coronal plane is higher, and there is a velocity gradient throughout the plane, though the gradient is more intense at the flow rate of 30 l/min compared with the flow rate of 7.5 l/min. The nasal cavity has a complicated geometry with tiny passages and more complicated morphology than other parts of the airways, which results in lower velocity in the vicinity of the walls and higher velocity at the core of the passage, which is associated with the velocity gradient observed on planes 1 and 2. Despite the nasal cavity, there is not a significant velocity gradient in the planes located downstream of the flow inside the airways. For instance, the distribution of momentum is more uniform on planes 4 and 5 located in the trachea. The velocity vectors are also included on the planes for velocity contours, demonstrating the direction of velocity on different planes. Plane 2 is located at the section after the pharynx, which means the flow should pass a curved path to reach this region. The presence of deflections on the fluid flow path generates another flow pattern perpendicular to the main flow called secondary flows. Among all other planes, the most secondary solid flows have been formed on plane 2 compared to other planes downstream of the airways (Fig. 5). The velocity magnitude at plane 5 near the tracheal exit is significantly higher at 30 l/min compared to 7.5 l/min. The flow field is uniform because of the less asymmetric structure of the trachea at this section.

The velocity profiles on lines 1–4 are depicted in Fig. 6 for two flow rates of 7.5 and 30 l/min. The velocity profiles support the data inferred by the velocity contours in Fig. 5. The velocity profiles on lines 1 and 2 are more asymmetric. Line 1 is located in the nasopharynx, being highly influenced by the complexity of the structure of the nasal cavities, and line 2 is located immediately after the larynx. The velocity is not distributed uniformly on this line. Lines 3 and 4 (locate in the trachea) are smoother than lines 1 and 2. These lines are highly affected by the complexity of the geometry in those regions. While comparison of the velocity profile for the two flow rates shows a significant difference in all lines.

Figure 7 shows the wall shear stress at various parts of the model for different flow rates. Wall shear stress is generated due to the velocity gradient on the layer of fluid adjacent to the wall, and the higher wall shear is observed due to the higher velocity gradient. The wall shear stress is the least in the nasal cavity and the larynx regions. The highest amount of wall shear stress occurs on the pharynx and the trachea. The velocity increment in the oropharynx and the nasopharynx, and the wall shear stress increases in these regions. In the larynx region, there is a great reduction in the wall shear; however, it increases in the trachea.

Figure 8 presents the wall shear stress at different planes. According to Fig. 8, the highest amount of shear stress occurs on plane 2, possibly due to the pharynx effect. After plane 2, planes 4, 3, and 5 have a higher shear stress rate. As mentioned previously, the shear rate is proportional to the velocity gradient. A higher inlet velocity generates a more significant velocity gradient, as shown in Figs. 7 and 8. The shear stress rate for the 30 l/min is higher than the 7.5 l/min conditions.

Figure 9 indicates pressure contours in different parts of the airways from the nasal cavity down to the trachea. The pressure is higher in the nasal cavity and reduces while going down in the airways because of the friction with walls and loss of energy. Pressure loss inside the airways could occur due to the simultaneous effect of wall friction and the complexity of the geometry. As the flow passes through a complicated geometry with tiny passages in some regions, flow movement would be accompanied by consecutive high expansion and contradictions, which affect the pressure. A significant change in pressure occurs in the regions with high variation in the cross section of the airways. One such instance is the pharynx, where both the oropharynx and nasopharynx passages join together, followed by a reduction in flow cross section. Accordingly, there is a pressure drop at the end of the oropharynx and nasopharynx. There is an expansion in the airways cross section, which results in a pressure increase. The larynx remains the same effect on the flow pressure and another sudden drop in the domain.

Figure 10 indicates this phenomenon more clearly. The pressure drops rapidly when moving from the first plane on the nasal cavity to subsequent planes on the flow path, and the pressure on plane 5 is the lowest observed. The pressure drop occurs more intensely at the flow rate of 30 l/min due to higher friction between the flow and the walls of the airways. Figure 10 represents this finding more quantitatively by quantitatively showing these results. The maximum pressure reduces gradually on successive planes along the flow path.

Figure 11 depicts the turbulence intensity in different regions for 7.5 and 30 l/min conditions. Turbulence intensity is the ratio of the velocity fluctuation deviation of the flow to the mean flow velocity and is a good measure to show the turbulence level. The turbulence intensity is the highest in the pharynx and the trachea, while there is a reduction in the turbulence intensity in the larynx region. A comparison of turbulence intensity in two flow rates of 7.5 and 30 l/min shows that flow rate is directly influential. The turbulence intensity in a higher flow rate is greater than in a lower flow rate.

Figure 12 shows the turbulent intensity on different sample planes, as indicated in Fig. 2 for different flow rates. The maximum turbulence intensity occurs on plane 2, which is followed by its rapid reduction on consecutive planes along the flow direction for 30 l/min. A similar trend is observed for 7.5 l/min with a slight difference on plane 4, having a little higher turbulence intensity. The turbulent intensity is generally related to the level of turbulent airflow flowing through the airways. The higher velocity magnitude will generate a higher turbulent level. This is obvious from the velocity contour from Fig. 5 that the velocity magnitude at P2 (plane 2) is higher than other areas (other planes). Furthermore, the vortex (refer to the direction of the vectors) is also found from this plane for both flow rates, while there is no vortex from other areas. Complex lung geometry is one of the main reasons for generating turbulent flow. Therefore, these reasons lead to a higher turbulent intensity at plan 2 than other planes.

Despite the turbulence intensity, turbulence kinetic energy is another parameter that represents the turbulence rate of the flow. The turbulence kinetic energy is higher in the pharynx and the trachea, with a significant reduction in the larynx region (Fig. 13).

A similar trend in the turbulence kinetic energy is observed on consecutive planes in the flow direction in Fig. 14. The turbulence kinetic energy is highest on plane 2, located at the pharynx, and rapidly reduces on other planes until plane 5, which has the lowest amount of turbulence kinetic energy compared with other planes. The flow rate is influential here, as the level of turbulence kinetic energy is higher for 30 l/min compared with the 7.5 l/min cases.

Flow streamlines can enhance our understanding of the flow field inside the airways and adequately characterize the flow behavior inside the complicated respiratory system passages. Velocity streamlines are depicted in Fig. 15 for two flow rates of 7.5 and 30 l/min. The highest velocity occurs at the inlet of the nasal cavity, the pharynx, and the trachea. The morphology of the airways substantially influences the streamlines, as in the regions where the airways have more complexity, and the streamlines have an intense disorder due to the deflections in the flow path. This is particularly recognizable in the nasal cavities, where the presence of complex and small passages results in the deflection of the streamlines and circulation of the airflow inside the maxillary sinuses. The circulation of the flow is more intense at the flow rate of 30 l/min [Fig. 15(b)]. The streamlines are highly affected by the local geometry of the model in the larynx, which is shown in a magnified view to better distinguish the streamlines in this region. The flow passing through the larynx gets disturbed abruptly due to the contradiction of the flow cross section and increment in the flow velocity. The local complexity of the geometry imposes disturbances to the flow. A higher disorder accompanies the higher flow rate in the flow streamlines. The highest velocity occurs in the nasopharynx, although the velocity magnitude is higher in the flow rate of 30 l/min.

Total deposition efficiencies of spherical, cylindrical, and tetrahedral microplastics with different diameters are presented in Fig. 16 for 7.5 and 30 l/min with various microplastics sizes, including 1.6, 2.56, and 5.56 μm, respectively. For 7.5 l/min, the deposition of 2.56 and 5.56-μm cylindrical microplastics is higher, while the 1.6-μm cylindrical particles have a lower deposition rate than others. At this flow rate, the deposition of spherical and tetrahedral microplastics is identical for all particle sizes. The deposition of 1.6 and 2.56-μm spherical and tetrahedral particles is quite similar, and the deposition rate increases for the 5.56-μm particles, while the cylindrical particles have rapid growth of deposition when increasing their diameter from 1.6 to 2.56 μm and from 2.56 to 5.56 μm. At a flow rate of 30 l/min, however, a different trend is observed. Figure 13(b) indicates that only for the particle diameter of 2.56 μm, the cylindrical microplastics have a higher deposition. Among the 1.6 and 5.56 μm particles, the spherical particles are higher than others. The slope of the deposition curve for the spherical particles is relatively constant by increasing the diameter. However, the slope of the deposition rate for the cylindrical and tetrahedral particles reduces by increasing the particle diameter from 2.56 to 5.56 μm, which means that the deposition of tetrahedral and cylindrical particles is more dependent on particle diameter at smaller sizes. The comparison of Figs. 16(a) and 16(b) demonstrates that the deposition of microplastics with different shapes and sizes is higher for similar particles at the flow rate of 30 l/min. However, the overall deposition efficiency of 7.5 l/min is higher than the overall deposition efficiency of 30 l/min. This is applied to all microplastic sizes. One of the main reasons for a higher deposition efficiency is the residence time to pass through the upper airway region from a lower flow rate, which is different from a higher flow rate. The lower flow rate leads microplastic to spend a longer time passing through the upper airway region compared to a higher flow rate. Furthermore, gravitational sedimentation and Brownian diffusion play an important role in addition to the inertial impaction effect. Brownian diffusion plays a more critical role at lower flow rates and gets weakened by increasing flow rates. Brownian diffusion occurs due to the collision of the particles with the air molecules and is more effective at lower flow rates.

Figure 17 indicates the deposition patterns for microplastics with different shapes at a flow rate of 7.5 l/min. There are several main hotspots that have the greatest amount of microplastic deposition. These hotspots are the nostrils, oropharynx, nasopharynx, and trachea. The oropharynx has the lowest concentration of deposited microplastics at a flow rate of 7.5 l/min. The dominance of inertial impaction is one of the particle deposition mechanisms. As particles have inertia when passing through the airways, when there is a change in the direction of the flow, the particles with high inertia may not be able to adjust their path with the fluid flow and consequently impact the wall. Despite the spherical and tetrahedral microplastics, the results are different from the cylindrical-shaped microplastics. In Fig. 17(c), at a flow rate of 7.5 l/min, only the particles with a diameter of 1.6 μm are trapped in the airways.

Figure 18 represents the distribution of microplastic particle deposited in the mouth-throat model at a flow rate of 30 l/min. In the flow rate of 30 l/min, the flow velocity is higher, and the particles have higher inertia compared with the flow rate of 7.5 l/min. Moreover, the particles have to pass through a steep curved path to enter the trachea, which is an accumulation of microplastics at the end of the oropharynx when the flow rate increases. This phenomenon is observed for particles with different shapes. The distribution of microplastics throughout the nasal cavity is more uniform at the flow rate of 7.5 l/min compared to the flow rate of 30 l/min, at which most of the particles have accumulated at the end of the oropharynx, and few particles are deposited on the nasal cavity walls. Figures 17 and 18 compared the deposition patterns of different shapes of microplastics at different flow rates. Figure 17(a) shows that most of the spherical microplastics with 1.6-μm diameter could pass the nasal cavity and the trachea and move toward the lower airways for 7.5 l/min. On the other hand, the 5.6-μm spherical particles get trapped in the upper airways at 30 l/min. However, this conclusion about the cylindrical and tetrahedral-shaped particles cannot be drawn. The overall deposition pattern demonstrates that most particles are deposited on the walls of the airways along the lung domain for 30 l/min. However, for the flow rate of 7.5 l/min, the particles are mostly trapped on the wall along the lung airways, especially on the oropharynx.

Figure 19 depicts the microplastic deposition hotspots for different shapes of microplastic. These figures provide a deeper understanding of how microplastic shapes influence the transport and deposition of microplastics in the airways and which regions are most prone to receiving these hazardous particles. Despite the nasal cavity, which is the most significant hotspot for the deposition of microplastics, the trachea is also another region where the deposition of particles is remarkable. In Fig. 19, particles accumulate downstream of the trachea, located between y = −0.06 and y = −0.1. In this region, with the presence of curvature in the trachea, the particles get trapped for not being able to follow the flow path properly.

Figures 20–22 depict the deposition hotspots for 1.6, 2.56, and 5.56-μm microplastics at the flow rate of 7.5 l/min. The deposition rate of particles of different shapes and sizes is the highest in the nostrils and the nasal cavity. This particular region is located between y = 0.04 and y = 0.07. The nasal cavity has a very complicated geometry, making it difficult for particles to follow the airflow path and deposit on its walls. However, the deposition rate in the upper airways for 1.6-micron tetrahedral microplastics at the flow rate of 7.5 l/min is slightly lower in the region between y = 0.04 and y = 0.07 (Fig. 20). A similar trend is reported for the 2.56 and 5.56-μm particles about the main deposition hotspots. In addition, most of the particles are deposited at the oropharynx region because of the dominance of inertial impaction mechanisms.

The deposition hotspots for 1.6, 2.56, and 5.56-μm microplastics at a flow rate of 30 l/min are shown in Figures 23–25. At the region between y = 0.07 and y = 0.09, which is associated with the upper parts and the top of the nasal, different amounts of deposition occur depending on the shape and size of particles. This difference in deposition is particularly more evident in Figs. 23–25 where a greater deposition rate occurs. However, when the particle size increases, the particles can barely reach the top wall of the nasal cavity, and y = 0.1 is free of particles or receives the least amount of particles at 30 l/min as demonstrated in Figures 23 and 24. Figures 23–25 present the deposition hotspots for the flow rate of 30 l/min. The deposition of spherical particles is lower when the particle size is smaller compared to the tetrahedral and cylindrical particles. Figures 23 and 24 (flow rate of 30 l/min) show that the deposition of spherical particles (1.6 and 2.56 μm particles) is much less than the tetrahedral and cylindrical particles, particularly in the region between y = 0.06 and y = 0.08. In this region, the particles enter the nasal cavity at the middle section of the air passages and have a further distance from the top and bottom walls of the nasal cavity. Consequently, they can more easily find their path toward the lower airways at the flow rate of 30 l/min. However, for larger spherical particles (5.56 μm), this phenomenon does not occur because the larger particles have higher inertia and easily get distracted from their path and impact the wall of the nasal cavity. Along with the inertia, the asymmetric shape of the microplastic also plays an important role in the overall deposition pattern. The rotational and translational movement of the non-spherical particle also influences the overall transport behavior. For Fig. 25, the amount of deposited 5.56-μm spherical particles in the region between y = 0.06 and y = 0.08 is higher.

Exposure to hazardous particles is one of the potential concerns of human health, especially during human inhalation. The main characteristic of tiny particles can transport to the deeper lung airways compared to the larger particle sizes. Microplastic exposure has become a health concern as it is recently found in the deep lung airway. The comprehensive investigation of microplastic transport and deposition is the key factor in increasing the understanding of its behavior. The present study analyzes the effect of microplastic transport within the upper airways under the three shapes of microplastics, including spherical, tetrahedral, and cylindrical shapes with 7.5 and 30 l/min and microplastic sizes as 1.6, 2.56, and 5.56 μm, respectively. The key findings are as follows:

• The nasal cavity area has a lower flow velocity compared to other areas. However, the flow fluctuates in this area, while the flow becomes uniform after passing the trachea instead. A higher flow velocity is found after passing the trachea for all flow rates. The flow velocity is affected by the lung geometry. The complex region leads to a higher flow velocity resulting in a higher level of turbulent flow. This leads to the difference of wall shear stress and turbulence intensity, which are related to the lung geometry and flow velocity. The wall shear stress at the nasal cavity and larynx regions is lower than in other regions. Higher turbulence intensity and kinetic energy are found at the complex parts, including the pharynx and the trachea.

• A high pressure is observed at the nasal cavity, mouth-throat, and trachea. The significant change in pressure is the result of a high variation in the cross section of the airways. Low pressure is found in the area where it is less complex.

• A higher microplastic deposition rate is observed in the nasal cavity compared to the trachea and other areas. The flow rates, shape, and size of microplastics influence the overall deposition pattern. A higher flow rate (30 l/min in the present study) leads to a lower deposition efficiency for all microplastic shapes.

• The highest deposition efficiency from the flow rate at 30 l/min is less than 13%, while there is around 18.5% of the deposition efficiency from the flow rate at 7.5 l/min. The lower flow rate (7.5 l/min) leads microplastic particles to spend a longer residence time passing through the upper airway region. The gravitational sedimentation and Brownian diffusion are also the main factors in addition to the inertial impaction effect. At 7.5 l/min, the deposition efficiency from tetrahedral and spherical shapes are similar and lower than the cylindrical shape. This is applied to all microplastic sizes. However, the deposition efficiency from spherical shapes is higher than other shapes at a 30 l/min flow rate. This is applied only for the larger microplastic size at 5.56 μm. Other smaller sizes have similar deposition efficiency.

• The airflow and particle transport are considered for an inhalation condition in the present study.

• The current model includes only the nasal cavity and oral airway as the inlets with the trachea section as the outlet.

• The dynamic wall motion or airway deformation is not considered in the present study.

• An open outlet condition and constant pressure are applied at the outlet.

• The effect of different environmental conditions such as humidity and temperature is not considered in the present study.

The results of this study could be useful for the health risk assessment under the exposure to hazardous particles topic. However, the comprehensive investigation of microplastic and its characteristic will be further studied to enhance the understanding of its effect on the human respiratory system. This includes the effect of microplastic with various shapes under the consideration of various microplastic sizes and flow rates. Moreover, the complex lung geometry with upper and lower lung generations will be further considered under the healthy and diseased airways conditions. The effect of environmental conditions, such as humidity and temperature, will also be considered in the future study.

The authors highly acknowledge the UTS high-performance computational facility.

The authors have no conflicts to disclose.

Mohammad S. Islam: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Md. Mizanur Rahman: Formal analysis (equal); Visualization (equal); Writing – review & editing (equal). Puchanee Larpruenrudee: Data curation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Akbar Arsalanloo: Formal analysis (equal); Visualization (equal); Writing – original draft (equal). Hamidreza Mortazavy Beni: Methodology (equal); Writing – review & editing (equal). Md. Ariful Islam: Resources (equal); Visualization (equal). Yuantong Gu: Project administration (equal); Supervision (equal); Writing – review & editing (equal). Emilie Sauret: Formal analysis (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.