Aerosol Particles in Lungs: Theoretical Modeling of Deposition and Mucociliary Clearance

Abstract

A theoretical model of the movement of aerosols in the lungs is proposed. The model is based upon the transport equations taking into account the aerosol inertial deposition processes. Particles move along curvilinear trajectories in self-twisted vortex air flows. Deposition occurs over several cycles of inhalation and exhalation. This mechanism works for particles with a diameter greater than 1–2 microns. All particles with a diameter of 4 microns and more are captured in the respiratory tract before the terminal bronchioles. Only particles with a size of less than 2 microns can penetrate into the respiratory parts of the lungs, but the cleaning coefficient for them is close to unity. The lung cleaning model describes the limiting capabilities of the mucociliarу system. The importance of taking into account the temporal characteristics of the mucociliary escalation of dust deposited in the lungs has been demonstrated. The existence of a mode of accumulation of particles in the lungs, due to a lack of cleaning time during periodic dust exposure, has been established.

1. Introduction

Rhinitis and pulmonary pathology frequently occur together. The epidemiologic, pathophysiologic, and clinical data are so compelling that the concept of “one airway, one disease” is accepted [1]. Several mechanisms have been held responsible for the interaction between the upper and lower airways (see Review [2]). A good example is the high prevalence of chronic rhinosinusitis in patients with asthma. Recent data provide ample evidence of existence of a systemic pathway between allergic rhinitis and asthma. It has been found that the severity of asthma directly correlates with the severity of sinus disease [3, 4]. Patients with chronic rhinosinusitis and asthma constitute the most severe form of unified respiratory tract disease, which is characterized by older age of the patients, greater duration of nasal symptoms, extent of sinus radiological changes, more prominent systemic inflammation markers, greater bronchial obstruction, and incidence of perennial allergic rhinitis [5].

Treating asthma with inhaled steroid therapy appears to be important in treating refractory chronic rhinosinusitis associated with asthma. In [6], the results of topical inhalation treatment of acute bacterial rhinosinusitis and standard systemic antibiotic therapy are compared. It has been found that topical inhalation therapy of acute bacterial rhinosinusitis may provide better treatment options, because systemic antibiotics can be associated with different adverse effects.

Aerosol in lungs are the task the solution of which is necessary in various areas. For example, in medical applications, aerosol therapy allows the medication to be delivered directly to those parts of the lungs where it should act.

In this chapter, we consider both deposition and clearance kinetics throughout the respiratory tract. Aerosol particles can penetrate the human lungs deep into the smallest respiratory tracts and parts of the lungs during inhalation. The effects, caused by inhaled particles, depend on the site at which they deposit within the respiratory system. Information concerning particle deposition in human lungs is of special interest because it is useful for a quantitative risk evaluation from exposure to aerosol particles. It is also vital to the effective administration of pharmaceutical aerosols by inhalation to targeted regions of the respiratory tracts. The deposition of particles in the lungs for given exposure conditions depends on parameters such as the airflow field, aerosol particle properties, breathing pattern, and geometric airway characteristics. There have been many investigations on this subject. Some recent reviews have been carried out and reported in [7, 8, 9]. In Section 3 of this chapter, we describe the aerosol inertial capture processes. Particles move along curvilinear trajectories in self-twisted vortex air flows. Deposition occurs over several cycles of inhalation and exhalation. The result is the distribution of deposition of inhaled aerosol particles in the airway generations in the respiratory tract of the lung.

Once deposited, specific clearance mechanisms will become effective, leading to a clearance of the deposited substances or a retention in different compartments of the respiratory tract. The main clearance mechanism for insoluble particles deposited in the conducting airways is via the mucociliary escalator. For the calculation of retention curves, this bronchial clearance model was applied to predicted deposition patterns for different particle sizes and breathing conditions. The importance of taking into account the temporal characteristics of the mucociliary escalation of aerosol particles deposited in the lungs has been demonstrated.

2. Model of human respiratory system

Morphometric models were initiated by E. Weibel [10] who assigned lengths and diameters to a symmetrical tracheobronchial tree (TBT) based on measurements of human airway casts. Weibel’s airway tree model assumes dichotomous branching and creates a fractal structure, which is physiologically consistent in the amount of air delivered to the terminal airways and the spatial arrangement of branches within the lung. A generation is defined as all airways occurring at a specific number of bifurcations from the trachea.

The number N (i) of the branches of the ith generation and the diameter D (i) are given by the formulas

$$\mathrm{N}\left(\mathrm{i}\right)=2^{\mathrm{i}} \mathrm{;} \mathrm{D}\left(\mathrm{i}\right)=\mathrm{D}\left(0\right)\ast 2^{-\mathrm{i}/3}$$ (1) E1

The airway diameter decreases distally, but because of the increasing number of tubes, the total cross section increases and the air velocity decreases. The law of change of air velocity follows from the requirement of conservation airflow before and after dividing branches of TBT:

$$\mathrm{V}\left(\mathrm{i}\right)=\mathrm{V}\left(0\right)\ast 2^{-\mathrm{i}/3}$$ (2) E2

The trachea (0th generation) has D(0) = 1.44 cm and V(0) = 290 cm/s. The branches of the 18th generation correspond to terminal bronchioles. The radius R(i) = D(i)/2, and length L(i) of each branch is assumed to be equal to three times the diameter.

In the adopted model, the time tp for the passage of the air channel is the same for all generations:

$${\mathrm{t}}_{\mathrm{p}}=\mathrm{L}\left(\mathrm{i}\right)/\mathrm{V}\left(\mathrm{i}\right)=3\ast \mathrm{D}\left(\mathrm{i}\right)/\mathrm{V}\left(\mathrm{i}\right)=3\ast \mathrm{D}\left(0\right)/\mathrm{V}\left(0\right)\approx 0.015sec.$$ (3) E3

We can compare this time with the characteristic inhalation time t_in = 1.25 s (we assume that the breathing frequency is 12 1/min, and for each period, there are four identical phase durations: inhale, pause, exhale, and pause). The ratio t_in/t_p ≈ 100 can be considered as the possible number of TBT generations that air could pass during inhalation. As was to be expected, this number is significantly greater than the total number of generations (≈ 25) of TBT. The passage of the airways to the level of the alveoli takes no more than 10% of the total inhalation time.

2.1 Air flow characteristics

An important dimensionless quantity in air mechanics is the Reynolds number (Re). It is used in the scaling of similar but different sized flow situations [11]. The value of Re is the ratio of inertial forces to viscous forces within an air which is subjected to relative internal movement due to different air velocities, which is known as a boundary layer in the case of a bounding surface such as the interior of an airway. In the air channel of TBT, the Reynolds number can be defined as

$$Re\left(\mathrm{i}\right)=\mathrm{V}\left(\mathrm{i}\right)\ast \mathrm{R}\left(\mathrm{i}\right)/\mathrm{\nu }$$ (4) E4

where ν is the air kinematic viscosity. The Reynolds number is important for different situations where a fluid is in relative motion to a surface. Matching the Reynolds numbers in different airways is sufficient to guarantee similitude of airflows. Table 1 shows the Reynolds numbers characterizing the air flow in the respiratory tract.

i=	0	2	4	6	8	10	12	14	16	18	20
Re =	1392.0	552.4	219.2	87.0	34.5	13.7	5.4	2.2	0.9	0.3	0.1

Table 1.

Reynolds numbers in the branches of the TBT.

It can be seen that in almost all branches of the TBT, the flow is laminar, however with Re > 1.

The behavior of particles suspended in an air flow is characterized [11] by Stokes number (Stk). It is defined as the ratio of the characteristic time of an aerosol particle to a characteristic time of the flow. In the air channel of TBT, the Stokes number can be defined as

$$\mathrm{Stk}=\left(1/18\right)\ast \left({\mathrm{\rho }}_{\mathrm{p}}/{\mathrm{\rho }}_{\mathrm{a}}\right)\ast \left({{\mathrm{d}}_{\mathrm{p}}}^2/\mathrm{\nu }\right)\ast \left[\mathrm{V}\left(\mathrm{i}\right)/\mathrm{R}\left(\mathrm{i}\right)\right]$$ (5) E5

where ρ_p and ρ_a are the particle and air densities and d_p is the particle diameter. Value R(i)/V(i) is the characteristic time of the flow. A particle with a low Stokes number follows fluid streamlines, while a particle with a large Stokes number is dominated by its inertia and continues along its initial trajectory. For Stk ≫ 1, particles will detach from a flow especially where the flow decelerates abruptly. For Stk ≪ 1, particles follow fluid streamlines closely.

In the accepted TBT model, the ratio of air velocity to diameter is the same for all generations; therefore, the Stokes number is the same in all branches, and it changes depending on the particle diameter. The corresponding results are shown in Table 2.

dp, μm	1	2	4	6	10	20	40
Stk =	0.0015	0.0060	0.024	0.054	0.149	0.60	2.39

Table 2.

Stokes numbers for particles of different diameters.

It can be seen that particles with a diameter of up to 10 μm (Stk ≪ 1) are involved in the movement of air. Only particles with a large diameter can “break away” from the flow and, for example, settle onto the wall directly due to the inertial mechanism.

2.2 Mucociliary escalator

The kinetics of mechanical particle transport from the TBT via the mucociliary escalator to the larynx cannot be measured directly. Average tracheal mucus velocities in healthy nonsmokers, as measured by noninvasive radiological techniques, ranged from 4 to 6 mm/min consequently; both ICRP [12] and NCRP [13] committees adopted a value of 5.5 mm/min. Because of the very limited experimental data on mucus clearance velocities in human bronchial airways, mucus velocity in other bronchial airways has to be calculated, using assumptions about TBT morphology and properties of the mucus flow in these airways. The law of change of mucus velocity U(i) follows from the requirement of conservation flows before and after dividing branches of TBT:

$$2\ast \mathrm{\pi }\ast \mathrm{D}\left(\mathrm{i}\right)\ast \mathrm{U}\left(\mathrm{i}\right)=\mathrm{\pi }\ast \mathrm{D}\left(\mathrm{i}-1\right)\ast \mathrm{U}\left(\mathrm{i}-1\right)$$ (6) E6

The diameter of each airway is determined by (1). So, we have

$$\mathrm{U}\left(\mathrm{i}\right)=\mathrm{U}\left(0\right)\ast 2-\text{2i}/3$$ (7) E7

This model uses an empirically derived expression for mucociliary transport rates in TBT.

For the calculation of retention characteristics, this bronchial clearance model was applied to predicted deposition patterns for different particle sizes and breathing conditions.

3. Inertial capture

Inertial capture of aerosol particles from air flow is illustrated by Figure 1.

Figure 1.

Trajectories of particles entering the air flow in a straight pipe with round cross section at distances r_с from the center. A cylindrical coordinate system (r,z) is used. The airway radius is R; length is L.

Particles entering the channel near the center fly through it during the time of the t_p ≈ 0.015 sec (see above Section 2.1.). Particles entering the channel near the wall move more slowly due to the effect of deceleration of the air flow in the boundary layer. During inhalation, they may not reach the end of the channel. They turn back during exhalation.

An important circumstance is the air swirling flow at the point of bifurcation of the branches of TBT [7]. Numerical studies on the bifurcation flow of lung models consistently identified different types of secondary flow pattern such as swirl-flow pattern or double vortex, also known as Dean flow. At inspiration, the rate of swirling flow at the inlet to the branch of TBT is 20–25% of the axial velocity.

In such a flow, the centrifugal force shifts the particles to the wall. Here, the particle is additionally inhibited. After several cycles of “inhale-exhale,” the particle will deposit on the wall. The capture of particles in some branch of the TBT occurs if the transit time of this branch is longer than the inhalation time tin≈ 1.25 s (see Section 1). The transit time increases with increasing distance from the channel axis. The particles that fly into the airway at a distance r > rc from axis are captured. This condition determines the probability P (i) of particle capture in the ith branch of the TBT:

$$\mathrm{P}\left(\mathrm{i}\right)=\left({\mathrm{\pi R}}^2-{{\mathrm{\pi r}}_{\mathrm{c}}}^2\right)/{\mathrm{\pi R}}^2=1-{{\mathrm{r}}_{\mathrm{c}}}^2/{\mathrm{R}}^2$$ (8) E8

This mechanism is the most effective even for small aerosol particles.

3.1 Air flow characteristics

There are numerous research studies focusing on airflow in the respiratory system (see Review [14]). The results of these studies provide an almost complete picture of the airflow in the lungs. It is easy to recognize the “classical” flow patterns in smooth cylindrical tubes (see [11]).

The axial velocity distribution at the entrance can be considered uniform. Further downstream, a Prandtl boundary layer forms along the walls. Its thickness is gradually increasing. In sufficiently long pipes at a distance Le from the entrance, the boundary layer fills the entire cross section. It establishes a Poiseuille flow with a parabolic distribution of the axial velocity along the radius r. The steady-state flow does not depend on the velocity distribution at the inlet. The distance at which the Poiseuille flow is formed is called the input length L_e. For L_e we have

$${\mathrm{L}}_{\mathrm{e}}\approx 0.1\ast \mathrm{D}\ast {\mathrm{R}}_{\mathrm{e}}$$ (9) E9

For a small (compared to the channel radius) thickness of the boundary layer, the axial air velocity V_z is described by the relation (Blasius solution):

$${\mathrm{V}}_{\mathrm{z}}={\mathrm{V}}_{\text{in}}\ast \mathrm{f}\left(\mathrm{\zeta }\right)$$ (10) E10

where V_in is the velocity far from the wall, which can be equated to the velocity at the entrance to the channel, and ζ is the self-similar variable:

$$\mathrm{\zeta }=\left(\mathrm{R}-\mathrm{r}\right)\ast {\left[{\mathrm{V}}_{\text{in}}/\mathrm{z}\mathrm{\nu }\right]}^{1/2}$$ (11) E11

The function f (ζ) is the result of a numerical solution of the continuity and Navier–Stokes equations after the transition to the self-similar variable ζ. The function f (ζ) initially grows linearly with increasing ζ and then approaches to 1. The scale of change is ζ₀ ≈ 3.

The swirling of the air flows was found in the numerical simulation of currents in the TBT channels [15, 16]. The distribution of the azimuthal velocity V_φ along the radius in the cross section of the channel is well interpolated by the function

$${\mathrm{V}}_{\mathrm{\phi }}=\mathrm{S}\ast {\mathrm{V}}_{\text{in}}\ast \left(\mathrm{r}/\mathrm{R}\right)\ast \left(1-\mathrm{r}/\mathrm{R}\right)$$ (12) E12

This function describes the real decrease in azimuthal velocity as it approaches the axis and the channel wall. The intensity of the vortex (swirl number) S is defined as the ratio of the average over the cross-sectional azimuthal velocity V_φ to axial veloсity V_z. At the entrance to the air channel of the TBD, the value of S = S0 ≈ 0.15–0.25 [15], and then S (z) decreases exponentially. The characteristic decay distance is close to the input length L_e [17].

3.2 Model of the capture of aerosol particles

This data is sufficient to make an estimate of the probability of the capture in the ith branch of the TBT. The motion of aerosol particles along the z axis is determined by the speed [Eq. (10)]. In the radial direction, the aerosol particle moves under the action of centrifugal force:

$$3\ast \mathrm{\pi }\ast {\mathrm{d}}_{\mathrm{p}}\ast {\mathrm{\rho }}_{\mathrm{a}}\ast \mathrm{\nu }\ast \left(\mathrm{dr}/\mathrm{dt}\right)=\mathrm{m}\ast {{\mathrm{V}}_{\mathrm{\phi }}}^2/\mathrm{r}$$ (13) E13

Time does not explicitly enter into Eqs. (10) and (13). They define (r, z) trajectory. In dimensionless variables x = z/R and y = r/R, this trajectory is described by the formula

$$\mathrm{y}={\mathrm{y}}_1/\left\{{\mathrm{y}}_1+\left(1-{\mathrm{y}}_1\right)\ast exp.\left[-\left(2/3\right)\ast {\mathrm{Г}}^2\ast \mathrm{Stk}\ast {\mathrm{\zeta }}_{\mathrm{o}}\ast {\mathrm{x}}^{3/2}/{Re}^{1/2}\right]\right\}$$ (14) E14

Here, y₁ = r₁/R is the dimensionless radius of the particle entrance. The full trajectory will be obtained after twisting (r, z) trajectories around the axis of the channel in accordance with the rotation of the particles with the azimuthal velocity V_φ along with the flow.

From Eq. (14), it follows that if Re ≫ 1 (initial generation of TBT) and/or Stk ≪ 1 (small particle size), then у ≈ у₁ all the way. In this instance, the particles do not shift to the wall. The shift is noticeable in those air channels where

$$\mathrm{\lambda }=\left(2/3\right)\ast {\mathrm{Г}}^2\ast \mathrm{Stk}\ast {\mathrm{\zeta }}_{\mathrm{o}}\ast {{\mathrm{x}}_{\mathrm{с}}}^{3/2}/{Re}^{1/2}>1$$ (15) E15

Even if λ ≪ 1, always y < 1, that is, direct deposition of particles due to their centrifugal demolition to the channel walls does not occur. However, the removal of particles into the zone of slow near-wall flow can significantly increase the time of flight of the particle through the channel. The definition of the transit time can be obtained, for example, by integrating the equation of motion for x (τ = t*V_in/R). The estimate of the channel transit time τ_c, taking into account the removal of a particle into the zone of a slow near-wall flow, has the form:

$${\mathrm{\tau }}_{\mathrm{с}}=\left\{\mathrm{\lambda }+\left[exp\left(\mathrm{\lambda }\right)-1\right]\ast {\mathrm{у}}_1/\left(1-{\mathrm{y}}_1\right)\right\}/\mathrm{Stk}\ast {\mathrm{Г}}^2$$ (16) E16

As λ increases to large values, the time span τc grows exponentially and can reach values greater than the inspiratory time. This can be considered the effect of trapping a particle in the corresponding channel. According to Eq. (16), the value of τ_c is different for particles entering the channel at different distances y₁ from the center of the channel. Only particles are captured that enter the channel at distances greater than that determined from the equation

$${\mathrm{у}}_1/\left(1-{\mathrm{y}}_1\right)=\left[{\mathrm{Г}}^2\ast \mathrm{Stk}\ast {\mathrm{\tau }}_{\text{in}}-\mathrm{\lambda }\right]/\left[exp\left(\mathrm{\lambda }\right)-1\right]$$ (17) E17

The dimensionless inspiratory time τ_in is defined in Section 2.1; it is ≈ 250.

The right-hand side of Eq. (17) varies with the generation number and with the characteristics (mass and size) of aerosol particles. Changes in the entrance radius of the capture of particles at the channel entrance lead to changes (according to Eq. (8)) of the probability of their settling in the ith generation.

3.3 Results

The results of the calculation of the probability of the capture of particles of different diameters in the generation of the respiratory tract TBT are presented in Figure 2.

Figure 2.

Dependence of the probability P of the capture of a particle (ordinates) with diameter d in the respiratory tracts of the ith generation (abscissa axis) of TBT.

It can be seen that the air cleaning mechanism due to the capture of aerosols in the airways is quite effective even for particles with a diameter of 1–2 μm. Particles with a diameter of 4 μm and more are almost completely captured in channels with a generation number < 19, that is, to respiratory branches of TBT. Particles of large diameters (>10 μm) are captured in the upper respiratory tract, that is, the trachea, the zonal extrapulmonary bronchi, and the intrapulmonary subsegmental bronchi.

The probability of capture (formula Eq.(8)) determines the proportion of captured particles from the total number of particles entering the channel. It is possible to determine the probability Pa (i) of the particles reaching the ith generation:

$$\mathit{Ра}\left(i\right)=\prod _{k=0}^{k=i-1}\left[1-P\left(k\right)\right]$$ (18) E18

The probability of reaching the deep parts of the lungs is almost zero for particles with sizes larger than 10 microns. Particles with sizes of 4–6 microns can penetrate into the lungs only up to the terminal bronchioles. Only particles with a size of less than 2 microns can penetrate into the respiratory branches of the TBT, but the cleaning coefficient for them is close to 1. These data also indicate a high efficiency of the capture of particles as a mechanism for cleaning air from aerosols.

Finally, we present data on the probability of trapping particles in the airway of the TBT. Denoting this value as Р_с (i), we can calculate it using the obvious formula

$${\mathrm{P}}_{\mathrm{c}}\left(\mathrm{i}\right)={\mathrm{P}}_{\mathrm{a}}\left(\mathrm{i}\right)\ast \mathrm{P}\left(\mathrm{i}\right)$$ (19) E19

The results of the calculation of P_c (i) are shown in Figure 3.

Figure 3.

The distribution of trapped particles of different diameters d by the numbers of generation of TBT.

The latter result gives a clear picture of the areas of the lungs in which the capture of aerosol particles of different diameters occurs. It is seen that even the smallest particles are mainly trapped in the bronchioles to the respiratory part. Large particles (with a diameter of 6 microns and more) do not penetrate into the terminal bronchioles at all.

4. Model of mucociliary clearance of the lungs

As a rule, the modification of the model of E. Weibel leads to more realistic conclusions about the work of the lungs. In addition, it allows us to describe a number of effects that are impossible in the “classical” TBT model. There are, however, negative aspects of such a modification. Asymmetric statistical models are analyzed only by numerical methods. At the same time, however, the visibility of the results and the possibility of generalization the conclusions are lost. Therefore, for the initial estimates, it is advisable to limit ourselves to the simpler initial morphometric model of E. Weibel.

4.1 Formulation of the problem

We introduce adequate variables to describe the dynamics of mucociliary cleaning of the lungs. We will describe the flow j (z, t) of the mucociliary escalation by the product of the velocity U(z) and the density n (z, t). The latter is the number of aerosol particles per unit length of the corresponding generation of TBT. Нere and later, t is time, and z is the coordinate along the TBT generations. It is counted from the larynx. Here, the advantage of considering generations instead of TBT branches is used explicitly. Generation is the entire set of TBT branches with the same division number. If you follow the flow of particles along the branches of the TBT, it will have gaps (vary approximately twice) at the points of division. If the flow is attributed to generation, it will be continuous. The continuity equation will be

$$\frac{\mathrm{\partial n}}{\partial }-\frac{\mathrm{\partial Un}}{\mathrm{\partial z}}=\mathrm{q}\left(\mathrm{x},\mathrm{t}\right)$$ (20) E20

Here, q (z, t) denotes the density of aerosol deposition on the inner surface of the TBT (per unit generation length and per unit time). Since precipitation occurs faster (in a few breath-exhalation cycles) than cleaning the lungs, sedimentation density can be factorized:

$$\mathrm{q}\left(\mathrm{z},\mathrm{t}\right)=\mathrm{Q}\left(\mathrm{z}\right)\ast \mathrm{H}\left(\mathrm{t}\right)$$ (21) E21

Here, the function H (t) describes the temporal dependence of the entry of aerosol particles into the lungs, and Q (z) is the distribution of the trapped particles along the TBT generations. The latter function is directly related to the probability of capture P_i found in Section 3. To determine the relationship, let us establish the correspondence of the z coordinate to the generation number i. Choosing the origin (z = 0) at the beginning of the trachea, we get that the end of the ith generation corresponds to the coordinate

$$\mathrm{z}={\mathrm{z}}_{\mathrm{i}}\equiv \sum _{\mathrm{k}=0}^{\mathrm{k}=\mathrm{i}}{\mathrm{l}}_{\mathrm{k}}$$ (22) E22

Here, l_k is the length of the kth generation. In the model of E. Weibel, it changes with the number k as well as the diameter: l_k = l₀*2^−k/3. The geometric progression (Eq. (22)) is summed. For z_i we have

$${\mathrm{z}}_{\mathrm{i}}={\mathrm{l}}_0\ast \frac{1-2^{-\left(\mathrm{i}+1\right)/3}}{1-2^{-1/3}}=\mathrm{Z}\ast \left(1-2^{-\frac{\mathrm{i}+1}{3}}\right)$$ (23) E23

Here, Z = l₀/ (1–2^−1/3) ≈ 29.1 cm is the total length of the TBT generations. The distribution of trapped particles Q (z) can be determined through the air exchange rate in the lungs W and the concentration of particles R in inhaled air:

$$\mathrm{Q}\left(\mathrm{z}\right)=\mathrm{W}\ast \mathrm{R}\ast {\mathrm{P}}_{\mathrm{i}}/{\mathrm{l}}_{\mathrm{i}}$$ (24) E24

We can express the mucociliary velocity U in terms of the z coordinate along the TBT generations:

$$\mathrm{U}\left(\mathrm{z}\right)={\mathrm{U}}_0\ast {\left(1-\mathrm{z}/\mathrm{Z}\right)}^2$$ (25) E25

where U0 ≈ 2.4 * 10⁻² cm/s is the velocity at the entrance to the TBT (near the larynx). With such values of the scale of length and velocity, the time scale is T = Z/U0 ≈ 1.2*10³ sec = 0.337 hour. The relation (Eq. (25)) is valid only at the bifurcation points, but in the future, it is also advisable to use it at intermediate values of z.

Instead of the density n (z, t), we introduce a new variable—the stream j (z, t) = n (z, t) * U (z), and instead of the z coordinate, a new variable

$$\mathrm{\vartheta }\left(\mathrm{x}\right)={\int }_0^{\mathrm{x}}\mathrm{dx}/\mathrm{U}\left(\mathrm{x}\right)$$ (26) E26

The meaning of this variable is the time required for the escalation of particles from depth z to the TBT. For the dependence U (z) given by formula (Eq. (25)), we obtain

$$\frac{\mathrm{\vartheta }}{\mathrm{T}}=\frac{\mathrm{z}/\mathrm{Z}}{1-\mathrm{z}/\mathrm{Z}};\frac{\mathrm{z}}{\mathrm{Z}}=\frac{\mathrm{\vartheta }/\mathrm{T}}{1+\mathrm{\vartheta }/\mathrm{T}}$$ (27) E27

From Eq. (20) we obtain the equation for the flow:

$$\frac{\mathrm{\partial j}}{\mathrm{\partial t}}-\frac{\mathrm{\partial j}}{\mathrm{\partial \vartheta }}=\mathrm{q}\left(\mathrm{\vartheta },\mathrm{t}\right)\ast \mathrm{U}\left(\mathrm{\vartheta }\right)$$ (28) E28

Its solution for the natural boundary conditions (j → 0 as ϑ → ∞ or t → - ∞) has the form:

$$\mathrm{j}\left(\mathrm{\vartheta },\mathrm{t}\right)={\int }_0^{\infty }\mathrm{q}\left(\mathrm{\vartheta }+\mathrm{u},\mathrm{t}-\mathrm{u}\right)\ast \mathrm{U}\left(\mathrm{\vartheta }+\mathrm{u}\right)\ast \mathrm{du}$$ (29) E29

Solution (Eq. (29)) of the problem of the dynamics of mucociliary clearance of the lungs allows us to estimate the parameters of escalation in several practically important cases.

4.2. Examples of solutions

Let us consider, for example, the dynamics of excretion of particles inhaled in a short time T0. In this case, the time dependence of the settling of particles can be chosen in the form H (t) = T₀*δ (t), where δ (t) is the Dirac impulse delta function.

The presence of the delta function allows us to determine the integral (29) and record the outflow of TBT (with z = 0) of the flow of settled particles in the form

$$\mathrm{j}\left(\mathrm{z}=0,\mathrm{t}\right)={\mathrm{T}}_0\mathrm{WR}{\left[\mathrm{U}\left(\mathrm{z}\right)\ast {\mathrm{P}}_{\mathrm{i}}/{\mathrm{l}}_{\mathrm{i}}\right]}_{\mathrm{z}=\mathrm{z}\left(\mathrm{t}\right)}$$ (30) E30

The content of square brackets in this ratio is taken for z, which varies with time according to Eq. (27). The results of calculations of the time dependence of the removal of particles of different diameters are shown in Figure 4.

Figure 4.

Dependence on time t (hour) of the particle flux at the exit of TBT.

As was to be expected, particles of large diameters (10–20 μm), which are captured in proximal sections of the TBT (generation with number i = 6–11), are removed in the first hours (5–9 hours). This can be explained by the high rate of mucociliary escalation (U = 10⁻³–10⁻⁴ cm/s). The maximum is clearly visible in the time dependence of the flow. Small particles (d = 1–4 μm) are removed much more slowly (1 day or more), and the time dependence of the flux decreases monotonously.

Integration over time of the effluent (30) gives us

$${\int }_0^{\infty }\mathrm{j}\left(\mathrm{o},\mathrm{t}\right)\mathrm{dt}={\mathrm{T}}_0\mathrm{WR}{\int }_0^{\mathrm{Z}}\left({\mathrm{P}}_{\mathrm{i}}/{\mathrm{l}}_{\mathrm{i}}\right)\mathrm{dz}=>\left(\mathit{dz}={\mathrm{l}}_{\mathrm{i}}\mathrm{di}\right)=>{\mathrm{T}}_0\mathrm{WR}\sum {\mathrm{P}}_{\mathrm{i}}={\mathrm{T}}_0\mathrm{W}$$ (31) E31

The last equality indicates the preservation of the number of aerosol particles: the total number of particles released (left side of Eq. (31)) is equal to the number of inhalation particles (right side of Eq. (31)).

We get other results in the case of prolonged (many hours) inhalation of dusty air. This situation is realized in the case of industrial dust exposure, which continues during the work shift. We assume that such an effect begins at t = 0 and stops at t = t0 = 8 hours. In this case, the function H (t) in the relation (Eq. (21)) can be given by the formula

$$\mathrm{H}\left(\mathrm{t}\right)=\mathrm{\eta }\left(\mathrm{t}\right)\ast \mathrm{\eta }\left({\mathrm{t}}_0-\mathrm{t}\right)$$ (32) E32

Here, η (t) is the unit Dirac function, which is zero for negative values of the argument and one for positive.

In this case, the solution (Eq. (29)) of Eq. (28) can be written as

$$\mathrm{j}\left(\mathrm{\vartheta },\mathrm{t}\right)=\left\{\begin{array}{cc}0& \text{if}\mathrm{t}<0\\ {\int }_{\mathrm{\vartheta }}^{\mathrm{\vartheta }+\mathrm{t}}\mathrm{Q}\left({\mathrm{\vartheta }}^{\prime }\right)\mathrm{v}\left({\mathrm{\vartheta }}^{\prime }\right){\mathrm{d\vartheta }}^{\prime }& \text{if}{\mathrm{t}}_0>t>0\\ {\int }_{\mathrm{\vartheta }+\mathrm{t}-{\mathrm{t}}_0}^{\mathrm{\vartheta }+\mathrm{t}}\mathrm{Q}\left({\mathrm{\vartheta }}^{\prime }\right)\mathrm{v}\left({\mathrm{\vartheta }}^{\prime }\right){\mathrm{d\vartheta }}^{\prime }& \text{if}\mathrm{t}>{\mathrm{t}}_0\end{array}\right\}$$ (33) E33

Interesting is the observed value of the flow of precipitated particles at the exit from the TBT, with ϑ = 0. If we return from the variable ϑ to the coordinate z, we get

$$\mathrm{j}\left(0,\mathrm{t}>0\right)={\int }_{{\mathrm{z}}_1}^{{\mathrm{z}}_2}\mathrm{Q}\left({\mathrm{z}}^{\prime }\right){\mathrm{dz}}^{\prime }$$ (34) E34

where

$${\mathrm{z}}_2\left(\mathrm{t}\right)={\mathrm{U}}_0\ast \mathrm{t}/\left[1+{\mathrm{U}}_0\ast \mathrm{t}/\mathrm{Z}\right];{\mathrm{z}}_1\left(\mathrm{t}\right)=0\text{for}\mathrm{t}<{\mathrm{t}}_0\text{and}{\mathrm{z}}_1\left(\mathrm{t}\right)={\mathrm{z}}_2\left(\mathrm{t}-{\mathrm{t}}_0\right)\text{for}\mathrm{t}>{\mathrm{t}}_{0.}$$

Formula (34) clearly demonstrates where the deposited particles are drawn from at a certain point in time.

Dependences on time t are depicted graphically in Figure 5. Along the abscissa, the time (in hours) after the start of inhalation of aerosol particles (10 μm in diameter) is plotted. Inhalation is assumed to occur over a period of time from t = 0 to t₀ = 8 hours. The boundaries z₁ and z₂ of the integration range are represented by the corresponding curves. The sedimentation density distribution of particles along the z coordinate (the integrand function Q (z) in the formula (Eq. (34)) is shown in the same graph. It can be seen that in the first hours after the start of inhalation, particles are removed from the TBT sections that are at a distance from the beginning to ≈ 0.8 * Z. A small part of the aerosol particles is captured here. Over time, the length of the cleaning section increases. After 5 hours, the removal of particles from a depth of ≈ 0.94 * Z begins. From here, a significant part of the trapped aerosol particles is removed. The maximum elimination occurs at t ≈ 10 hours after the start. After this, the boundaries of the dust removal area are shifted further into TBT, where relatively little dust accumulates. By 15–20 hours after the start of inhalation (7–12 hours after the end), almost all captured particles are removed from TBT.

Figure 5.

The dependence on time t (hour) of the boundaries of removal from the TBT trapped particles with a diameter of 10 mkm.

Thus, breaks lasting 7-12 hours eliminate the risk of accumulation of aerosol particles with a diameter of 10 mkm in the lungs. It is assumed, of course, that the mass of the trapped aerosol particles is not too large, so that its removal does not exceed the mass transfer capabilities during mucociliary escalation.

Cleaning the lungs from smaller dust particles occurs differently. The most probable capture depth is ≈ 0.98 * Z for particles with g = 4 μm. In these generations of TBT, the rate of mucociliary escalation is 10 times less than in the trachea. Cleaning time is much longer.

The results of direct calculations of the flow at the TBT outlet confirm these conclusions (see Figure 6).

Figure 6.

The time dependence of the elimination of aerosol particles trapped in 8 hours of work.

It is seen that particles with d = 20 μm practically do not stay in the lungs; they come out of TBT as they arrive, of course, if the mass of trapped particles does not exceed the mass transfer capacity during mucociliary escalation. Particles with d < 4 μm are deposited in the deep sections of TBT. They do not have time to get out of the lungs during the break. Their share exceeds 30%. This indicates the possibility of accumulation in the lungs of trapped particles of small diameter. For almost complete removal of such particles, it is necessary not less than 2 days. More time is needed to remove smaller particles.

The fact of accumulation of small particles during periodic inhalation of dusty air (during working hours), alternating with periods of rest when mucociliary clearance of the lungs occurs, raises the question of the average number of particles characterizing the result of such accumulation. Using solution (Eq. (33)) of Eq. (20), which describes the dynamics of mucociliary cleaning of the lungs, we can determine the density n (z, t) of the captured particles n(z, t) = j (z, t) / U (z). And after integrating by TBT generation, we have the total amount of N (t) particles in the lungs:

$$\mathrm{N}\left(\mathrm{t}\right)={\int }_0^{\mathrm{Z}}\mathrm{dz}\frac{\mathrm{j}\left(\mathrm{z},\mathrm{t}\right)}{\mathrm{U}\left(\mathrm{z}\right)}={\int }_0^{\infty }\mathrm{d\vartheta j}\left(\mathrm{\vartheta },\mathrm{t}\right)={\int }_0^{\infty }\mathrm{d\vartheta }{\int }_0^{\infty }\mathrm{d\tau q}\left(\mathrm{\vartheta }+\mathrm{\tau },\mathrm{t}-\mathrm{\tau }\right)\mathrm{v}\left(\mathrm{\vartheta }+\mathrm{\tau }\right)$$ (35) E35

Obvious substitutions of order and integration variables simplify this expression to the form:

$$\mathrm{N}\left(\mathrm{t}\right)={\int }_0^{\infty }\mathrm{d\vartheta U}\left(\mathrm{\vartheta }\right)\mathrm{Q}\left(\mathrm{\vartheta }\right){\int }_0^{\mathrm{\vartheta }}\mathrm{d\tau H}\left(\mathrm{t}-\mathrm{\tau }\right)$$ (36) E36

Now, the time dependence of the capture of particles H (t) is a sequence of rectangular pulses of the form (Eq. (32)), duration t0 (working time), and periodicity t₁ (1 day). A rigorous computation of the double integral (Eq. (36)) is possible, but not interesting. It is clear that this is a certain function of time, which varies with the period of the day around a certain mean value <N>. It can be estimated quite accurately. Indeed, the calculation of the average < N > in time from N (t) reduces to replacing in the integral H (t-τ) by <H (t-τ)>. The last value does not depend on τ; therefore, from Eq. (36) we can get

$$<N>=<H>\ast {\int }_0^{\infty }\mathrm{d\vartheta v}\left(\mathrm{\vartheta }\right)\mathrm{Q}\left(\mathrm{\vartheta }\right)\mathrm{\vartheta }=<H>\ast {\int }_0^{\mathrm{X}}\mathrm{dxQ}\left(\mathrm{x}\right)\mathrm{\vartheta }\left(\mathrm{x}\right)$$ (37) E37

As in formulas (Eq. (26–27)), the meaning of the variable ϑ (z) is the time required for the escalation of particles from depth z to the beginning of the TBT. The same technique that was used in the transition in formula (Eq. (21)) from the integral to the sum over the generation numbers leads to the expression

$$<N>=<H>\ast W\ast R\ast <\vartheta >$$ (38) E38

Here, the average exit time < ϑ > is determined by the probability distribution of the capture of particles of the corresponding diameter:

$$<\vartheta >=\sum {\mathrm{P}}_{\mathrm{i}}\ast {\mathrm{\vartheta }}_{\mathrm{i}}$$ (39) E39

The relation (Eq. (38)) is quite clearly interpreted: the average number of aerosol particles in the lungs is equal to their number in the volume of air, which is determined by pulmonary air exchange over the average time of exiting <ϑ > .

The obtained values of the average time of exit are in complete agreement with the results on the time dependence of the output stream for various inhalation regimes given above. The average exit time for particles of different diameters is given in Table 3.

d, mkm	1	2	4	6	10	20
<ϑ>,hours	27.7	23.7	12.3	7.4	3.6	1.1

Table 3.

The average exit time for particles of different diameters.

5. Conclusions

We investigated the deposition, retention, and clearance of inhaled aerosol particles in the lungs. Each of these processes is intensively investigated in numerous papers. As a rule, for their analysis, complex computer simulation of air and hydrodynamic flows in the branches of TBT is used.

The phenomenological model we have constructed is a simplified picture of real processes. It is used to study the key aspects of aerosol behavior in the lungs. This is a self-consistent and closed set of assumptions about the process, reflecting the basic, though not all, properties of a real phenomenon. Simulation is a research method with several goals:

• Structuring existing knowledge, giving it a certain shape, turning a data set into some information design

• The use of already accumulated information to determine the priority directions of its detailed elaboration, ranking of new information

The model reflects the unity of the main functional structures and, at the same time, the specifics of the processes occurring in them. Despite the descriptive nature of the proposed modeling, its meaning is that it gives the plot accuracy and compositional clarity, which are not inherent in real objects of modeling.

Conflict of interest

The author does not have a conflict of interest.