The effect of alveolar mixing on particle retention and deposition investigated by a dynamic single-path model

Abstract

A dynamic, single-path model is developed for dry powder transport in the lungs. The model differentiates between particle behavior in the respiratory bronchioles and alveolar ducts on one hand and inside the alveoli on the other. In particular, it considers the alveolar volume of each generation as a mixing chamber. Air inflow to the alveoli is calculated by accounting for the deformation of airways during breathing. Particle dispersion along the respiratory tract is taken into account and mechanistic deposition rates are developed for the alveoli. Deposition by Brownian diffusion is modeled by a concentration boundary layer, whose thickness varies inversely with the intensity of mixing. The plausibility of the assumption of alveolar mixing is tested indirectly by comparison of model predictions with benchmark data of the exhaled concentration profile and of the pulmonary deposition of continuously inhaled aerosols. The observed agreement lends support to the hypothesis that alveolar mixing represents fundamental physics of the breathing process. It also supports the suggestion that alveolar mixing provides an additional axial dispersion mechanism in the acinus, which is independent of particle size and is active at zero gravity.

Copyright © 2020 American Association for Aerosol Research

EDITOR:

• Warren Finlay

1. Introduction

Prediction of the characteristics of particle transport and of the spatio-temporal variation of deposition in the respiratory tract has been a central research goal of the past decades, motivated by efforts to assess the consequences of air-born pollutants and to optimize the location of inhaled drug delivery (Heyder 2004; Kleinstreuer and Zhang 2010; Longest and Holbrook 2012).

Key parameters that affect aerosol transport and deposition are related to the particle itself (size, density, shape, charge, hygroscopicity), the breathing pattern (inhaled volume, flow rate and its time-variation during inhalation and exhalation) and the geometrical and dynamical characteristics of the respiratory tract, i.e., airway diameters and lengths at functional residual capacity, as well as their temporal variation, and size and distribution of alveoli.

Particle deposition along conducting airways has been typically modeled in the idealized tubular geometry and reliable correlations are available (Finlay 2001). Alternatively, three-dimensional simulations–based on the actual shape of upper airways as recorded by tomography and/or acinar models derived from morphometrical measurements–are increasingly becoming available (Khajeh-Hosseini-Dalasm and Longest 2015; Kolanjiyil and Kleinstreuer 2016; Hofemeier et al. 2018). However, because the computational load is excessive, such approaches usually concentrate on one region of the lung, or provide global predictions in combination with a model of the rest of the lung that is not simulated in detail.

Deposition in the pulmonary region, i.e., along the walls of the respiratory bronchioles and inside the alveoli, appears to be understood less satisfactorily and many open questions remain. Preliminary efforts to deal with alveolar volumes as extensions of the airways of the respective generation (Egan and Nixon 1985; Darquenne and Paiva 1994) were met with some success. A major improvement in the modeling strategy was the consideration by Choi and Kim (2007) of a transport scheme that differentiates between airway and alveolar volumes and includes interchange of particles between the two.

More than two decades ago, it was first proposed that the flow inside an alveolus may result in chaotic mixing despite the very small velocities involved (Tsuda, Henry, and Butler 1995; Henry, Butler, and Tsuda 2002). Chaotic mixing was argued to be the outcome of the interaction of axial flow due to expansion/contraction with recirculation due to shear imposed by the airway flow. This idea was subsequently investigated by three-dimensional simulations and experiments (Tsuda et al. 2002, Tsuda, Henry, and Butler 2008, Tsuda, Laine-Pearson, and Hydon 2011; Sznitman et al. 2009; Henry et al. 2012), and was more recently enriched by a complementary mixing mechanism (Darquenne and Prisk 2003; Sznitman et al. 2009; Ma and Darquenne 2012; Fishler et al. 2017) based on streamline crossing by the particles in the presence of high velocity gradients. The idea of chaotic mixing was exploited in a simple filtration model that assumed pure convection along airways in combination with perfect mixing inside the alveolar volume of each generation (Georgakakou et al. 2016), and which–despite its simplicity–produced reasonable results (Gourgoulianis et al. 2017). It is noted that that such flow-based mixing is predominantly confined to the proximal generations of the acinus (e.g., Fishler, Mulligan, and Sznitman 2013).

The present work focuses on the pulmonary region and interrogates the dominant patterns of mixing and their effect on particle retention and deposition. An Eulerian, single-path deposition model (Taulbee and Yu 1975; Egan et al. 1989; Darquenne and Paiva 1994; Choi and Kim 2007) is used, which assumes symmetrically dichotomous airways. This class of models can handle time-varying deposition during realistic breathing cycles, in contrast to filtration–or Lagrangian–models (Kim et al. 1983; Martonen 1993), which apply only to steady flow conditions.

A key aspect of the present dynamic model of particle transport and deposition lies in the differentiation between the dynamics of bronchiolar airways and alveolar ducts on the one hand and alveolar volumes on the other. Indeed, the possibility that particle concentrations in respiratory airways and in alveoli follow fundamentally different dynamics has not attracted much attention in the literature, apart from the work of Choi and Kim (2007).

Flow in the ducts is described by one-dimensional convection with superimposed axial dispersion (Darquenne and Paiva 1994), whereas the alveolar volumes of each generation are cumulatively considered as an efficient mixing chamber (Georgakakou et al. 2016). The results indicate the consequences of the assumption of perfect alveolar mixing and–through comparison with benchmark data–serve as an indirect test of the plausibility of this assumption.

Predictions by the model are compared to two categories of independent data: (i) Exhaled concentration profiles of a continuously inhaled aerosol (Choi and Kim 2007). Exhalation data, and in particular the long-time tail, represent a sensitive indicator of the retention pattern and release rate of particles from the acinus. (ii) Pulmonary deposition of a continuously inhaled aerosol, as estimated (after appropriate correction, ICRP 1994) by the slow-cleared fraction (Heyder et al. 1986). These data span a wide range of particle sizes and breathing patterns and thus represent a stringent test of any deposition model.

It is demonstrated that predictions are in satisfactory agreement with all the above data. Based on this favorable comparison, it is argued that the differentiation between the dynamics of open-ended ducts (conducting/respiratory bronchioles, alveolar ducts) and closed volumes (alveoli), and in particular the hypothesis of good mixing in the latter, may represent fundamental physics of the breathing process.

2. Methods

2.1. Modeling of time-varying lung geometry and air flow

The respiratory tract is modeled as a symmetrically bifurcating network of the trachea and N = 23 generations of flexible cylindrical airways, with the diameter, $d_{0,i},$ and length, $l_{0,i},$ of each generation $i$ at functional residual capacity, given by one of the standard available models (Finlay 2001). The same model provides reference values of $\mathit{FRC}$ and the total alveolar volume per generation,$V_{a0,i}.$ All alveoli are assumed to have the same size and expand/contract at the same rate. The fraction of the airway wall covered by alveoli, $f_i,$ is an increasing function of generation number, starting from the first alveolated generation.

In order to take into account the variation of airway diameter during breathing, we resort to a generally accepted model of airway elasticity (Lambert et al. 1982) and use it in combination with the compliance properties of the lung. More specifically, the airway diameter is expressed as a function of intramural pressure, ΔP, by the following two hyperbolae, matched at $\Delta P=0:$ (1) $$\frac{d_i^2}{d_{i,\mathit{max}}^2}={\alpha }_0{\left(1-\frac{\Delta P}{P_1}\right)}^{-n_1}\mathrm{ },\mathrm{ }\Delta P\le 0$$(1) (2) $$\frac{d_i^2}{d_{i,\mathit{max}}^2}=1-\left(1-{\alpha }_0\right){\left(1-\frac{\Delta P}{P_2}\right)}^{-n_2}\mathrm{ },\mathrm{ }\Delta P>0$$(2)

Here,$P_1={\alpha }_0{n}_1/{\alpha }_0^{\prime }, P_2=-n_2\left(1-{\alpha }_0\right)/{\alpha }_0^{\prime }, \mathrm{ }{d}_{i,\mathit{max}}$ is the maximum airway diameter of generation i (i.e., at the end of full inspiration) and ${n_1, n}_2,{\alpha }_0,{\alpha }_0^{\prime }$ are the elasticity parameters of each generation (Lambert et al. 1982). Note that the original data have been presently extrapolated by using the parameters of generation 16 also for the airway generations 17–23.

The intramural pressure at time t is related to the instantaneous lung volume,$V_L(t),$ according to the following standard expression of lung elasticity (Pollak 1998) (3) $$\Delta P\left(t\right)=\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{C_{L0}}\mathrm{ }ln\left(\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{V_{L,\mathit{max}}-V_L\left(t\right)}\right)$$(3) where $C_{L0}$ is the lung compliance at $\Delta P=0$ and $V_{L,\mathit{max}},V_{L,\mathit{min}}$ are defined in terms of TLC and RV ($V_{L,\mathit{max}}=\mathit{TLC}+0.1L, V_{L,\mathit{min}}=RV-0.1L$) so that Equation (3)(3) $$\Delta P\left(t\right)=\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{C_{L0}}\mathrm{ }ln\left(\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{V_{L,\mathit{max}}-V_L\left(t\right)}\right)$$(3) does not produce singularities in the breathing range. Here, quasi-equilibrium is assumed and also pressure drop along the airways is neglected, i.e., $\Delta P$ is defined by taking the intramural pressure equal to the alveolar recoil pressure according to Equation (3)(3) $$\Delta P\left(t\right)=\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{C_{L0}}\mathrm{ }ln\left(\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{V_{L,\mathit{max}}-V_L\left(t\right)}\right)$$(3) . The maximum airway diameters, $d_{i,\mathit{max}},$ are calculated from the values at FRC of the model used, by substituting $V_L\left(t\right)=\mathit{FRC}$ in Equation (3)(3) $$\Delta P\left(t\right)=\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{C_{L0}}\mathrm{ }ln\left(\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{V_{L,\mathit{max}}-V_L\left(t\right)}\right)$$(3) and then using the value of $\Delta P$ in Equation (2)(2) $$\frac{d_i^2}{d_{i,\mathit{max}}^2}=1-\left(1-{\alpha }_0\right){\left(1-\frac{\Delta P}{P_2}\right)}^{-n_2}\mathrm{ },\mathrm{ }\Delta P>0$$(2) , together with $d_i=d_{0,i}.$

The present model admits arbitrary, generally time-varying, inhalation and exhalation flow rates, $Q_{0, \mathit{inh}}(t)\ge 0$ for ${0<t<t}_{\mathit{inh}}$ and $Q_{0, \mathit{exh}}(t)\le 0$ for ${0<t<t}_{\mathit{exh}},$ which are defined at the mouth. Thus, the instantaneous lung volume, needed for Equation (3)(3) $$\Delta P\left(t\right)=\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{C_{L0}}\mathrm{ }ln\left(\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{V_{L,\mathit{max}}-V_L\left(t\right)}\right)$$(3) , is calculated by simple mass balances of entering and leaving air. For example, in a breathing cycle starting from $\mathit{FRC},$ the lung volume during inhalation and exhalation is given by the following Equations (4)(4) $$V_{L,\mathit{inh}}(t)=\mathit{FRC}+\underset{0}{\overset{t}{\int }}{Q}_{0,\mathrm{\hspace{0.17em}}\mathit{inh}}\left(\tau \right)d\tau$$(4) and (5)(5) $$V_{L,\mathit{exh}}(t)=V_{L,\mathit{inh}}(t_{\mathit{inh}})+\underset{0}{\overset{t}{\int }}{Q}_{0,\hspace{0.17em}\mathit{exh}}\left(\tau \right)d\tau$$(5) , respectively: (4) $$V_{L,\mathit{inh}}(t)=\mathit{FRC}+\underset{0}{\overset{t}{\int }}{Q}_{0,\mathrm{\hspace{0.17em}}\mathit{inh}}\left(\tau \right)d\tau$$(4) (5) $$V_{L,\mathit{exh}}(t)=V_{L,\mathit{inh}}(t_{\mathit{inh}})+\underset{0}{\overset{t}{\int }}{Q}_{0,\hspace{0.17em}\mathit{exh}}\left(\tau \right)d\tau$$(5)

Next, the time-varying characteristics of the alveolar volumes are considered. Inhalation is used as an example, but the exact same reasoning and equations apply for exhalation. Air flow during inhalation is driven by the expansion of both the airways and the alveoli. By combining Equations (1)–(5), we obtain the instantaneous values of airway diameters at every time instant. Thus, we may calculate the part of air flow rate, $Q_{d, i},$ that is accommodated at each generation because of the diameter expansion (6) $$Q_{d,i}\left(t\right)=\frac{d}{dt}\left(2^i{l}_{0,i}\frac{\pi d_i^2}{4}\right)$$(6)

The difference between the inflow rate at the mouth and the sum of the contributions of all the airway expansion rates constitutes the total expansion rate of all the alveolar volumes. Thus, the air flow rate that enters the alveolar volume of generation $i$ is calculated as: (7) $$Q_{a,i}\left(t\right)=\left(\frac{V_{a0,i}}{\sum _{j=0}^N{V}_{a0,j}}\right)\left(Q_0-\sum _{j=0}^N{Q}_{d,j}\right)$$(7)

The first parenthesis in Equation (7)(7) $$Q_{a,i}\left(t\right)=\left(\frac{V_{a0,i}}{\sum _{j=0}^N{V}_{a0,j}}\right)\left(Q_0-\sum _{j=0}^N{Q}_{d,j}\right)$$(7) expresses the assumption that all alveoli expand at the same rate.

Based on the above, the remaining flow and anatomic parameters may be calculated. In particular, the air flow rate entering the airways of generation $i$ is: (8) $$Q_{i,\mathit{inlet}}(t)=Q_0-\sum _{j=0}^{i-1}\left(Q_{d,j}+Q_{a,j}\right)$$(8) and its variation along the airway is taken linear in the distance, x, from the entrance (9) $$Q_i(x,t)=Q_{i,\hspace{0.17em}\mathit{inlet}}-\frac{x}{l_k}\left(Q_{d,i}+Q_{a,i}\right)$$(9)

Also, the instantaneous diameter, $D_a,$ of the alveoli (taken as spherical cavities) may be calculated by dividing the total alveolar volume by the total number of alveoli. The instantaneous alveolar volume of generation $i$ is given by the expression (10) $$V_{a,i}(t)=V_{a0,i}+\underset{0}{\overset{t}{\int }}{Q}_{a,\hspace{0.17em}i}\left(\tau \right)d\tau$$(10)

2.2. Modeling of spatio-temporal variation of particle concentration

The concentration, $c(x,t),$ of particles along the airways is described by a one-dimensional differential mass balance that includes convection and axial dispersion, as well as source and sink terms that correspond to mass exchange with the alveoli and to deposition. The concentration of particles inside the alveoli of each generation, $c_{a,i}(t),$ is taken as spatially uniform, and thus the respective mass balances depend only on time. More specifically, the following equations are formulated for generation $i$ during inhalation: (11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) (11b) $$\frac{d}{dt}\left(V_{a,i}{c}_{a,i}\right)=\frac{Q_{a,i}}{l_{0,i}}\underset{0}{\overset{l_{0,i}}{\int }}\mathit{cdx}+f_i\underset{0}{\overset{l_{0,i}}{\int }}{L}_{i}dx-L_{a,i}$$(11b)

The respective mass balances during exhalation are: (12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) (12b) $$\frac{d}{dt}\left(V_{a,i}{c}_{a,i}\right)=-Q_{a,i}{c}_{a,i}+f_i\underset{0}{\overset{l_{0,i}}{\int }}{L}_{i}dx-L_{a,i}$$(12b)

Note that $A_i$ is the total airway flow cross section and $V_{a,i}$ the total alveolar volume of generation $i.$ Also, $f_i$ is the fraction of airway wall that has been replaced by alveolar openings and $b$ is the fraction of the area of the spherical alveolus that is missing because of the same opening. Term $L_i$ is the deposition rate per unit length of the cylindrical airways of generation $i$ and $L_{a,i}$ the deposition rate in the spherical alveolar volumes of generation $i.$

The terms on the rhs of Equations (11a)(11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) and (12a)(12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) express respectively the contribution of axial convection, axial dispersion, radial convection to/from the alveolar volume, deposition along the airway and deposition on the alveolar walls. The terms on the rhs of Equations (11b)(11b) $$\frac{d}{dt}\left(V_{a,i}{c}_{a,i}\right)=\frac{Q_{a,i}}{l_{0,i}}\underset{0}{\overset{l_{0,i}}{\int }}\mathit{cdx}+f_i\underset{0}{\overset{l_{0,i}}{\int }}{L}_{i}dx-L_{a,i}$$(11b) and (12b)(12b) $$\frac{d}{dt}\left(V_{a,i}{c}_{a,i}\right)=-Q_{a,i}{c}_{a,i}+f_i\underset{0}{\overset{l_{0,i}}{\int }}{L}_{i}dx-L_{a,i}$$(12b) express respectively convective exchange between airway and alveoli, and interchange caused by deposition inside the airway and inside the alveoli.

Though the exact calculation of deposition terms will be described in the following section, it is important to observe here that these terms have non-trivial contributions to the mass balances. Of the total deposition rate, $L_i,$ per unit length of a cylindrical airway, only $\left(1-f_i\right)L_i$ actually arrives on the wall, while $f_i{L}_i$ enters the alveolar volume through the openings. Similarly, out of the total deposition rate $L_{a,i}$ inside spherical cavities, only $(1-b)L_{a,i}$ actually arrives at the alveolar walls, while $b{L}_{a,i}$ returns to the airway through the openings. A schematic representation of the various particle fluxes is given in Figure 1.

Figure 1. Schematic of the various particle fluxes between airways and alveoli.

The above analysis may be summarized by stating that particles move between the airways and the alveolar volumes of each generation by two mechanisms whose contribution is additive, (a) by convection with the air flow entering or leaving the alveoli and (b) by gravitational settling and/or Brownian diffusion on the fraction of the airway wall that is absent due to the alveolar openings. This summation is equivalent to the classical formulation of mass flux of a chemical species in a mixture as the sum of a convective and a molecular component.

According to the present model, dispersion along the respiratory tract is accounted for by an explicit and an implicit contribution. The former is dominant in the conducting airways and gradually diminishes distally. The latter is a direct consequence of alveolar mixing and thus is active only in the acinus. This line of thought is in accord with the results of recent stochastic modeling of particle transport in the lung (Hofmann, Pawlak, and Sturm 2008), which showed that axial diffusion in conducting airways and mixing in alveoli constitute the dominant contributions to dispersion.

The explicit contribution is the effective diffusivity, $\mathcal{D}$ in Equations (11a)(11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) and (12a)(12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) , which is taken to be the combined effect of Brownian diffusion and convective mixing. Following Choi and Kim (2007), we write (13) $$\mathcal{D}={\mathcal{D}}_{Br}+{\mathcal{D}}_{\mathit{mix}}=\frac{C_{C}kT}{3\pi \mu d_p}+0.5u_i{d}_i=\frac{C_{C}kT}{3\pi \mu d_p}+\frac{2Q_i/2^i}{\pi d_i}$$(13)

The first term on the rhs of Equation (13)(13) $$\mathcal{D}={\mathcal{D}}_{Br}+{\mathcal{D}}_{\mathit{mix}}=\frac{C_{C}kT}{3\pi \mu d_p}+0.5u_i{d}_i=\frac{C_{C}kT}{3\pi \mu d_p}+\frac{2Q_i/2^i}{\pi d_i}$$(13) is the classical expression for particles of diameter $d_p$ in a fluid of absolute temperature $T$ and viscosity $\mu ,$ including the Cunninghum correction factor, $C_C.$ The second term is derived from the popular expression (Taulbee and Yu 1975; Ultman and Thomas 1979) ${\mathcal{D}}_{\mathit{mix}}=0.167u_i{l}_i,$ where the mean velocity is calculated from the instantaneous values of flow rate and airway cross-section and the length is substituted by the diameter using the approximate rule $l_i/d_i\cong 3.$ The above slight modification, which is reminiscent of a formula proposed by Scherer et al. (1975) was chosen in order to take consistently into account the temporal variation of diameter and flow rate. However, it turns out that it gives indistinguishable results from the original form.

Current understanding of dispersion in the lungs focuses on the role of bifurcations (Sarangapani and Wexler 1999) and on the local variability of particle paths caused by the coupling of cross-stream velocity differences and gravitational settling (Darquenne and Prisk 2003). The expression for convective mixing in Equation (13)(13) $$\mathcal{D}={\mathcal{D}}_{Br}+{\mathcal{D}}_{\mathit{mix}}=\frac{C_{C}kT}{3\pi \mu d_p}+0.5u_i{d}_i=\frac{C_{C}kT}{3\pi \mu d_p}+\frac{2Q_i/2^i}{\pi d_i}$$(13) is presently applied as an empirical correlation, which accounts only for these effects. The gravity mechanism predicts that dispersion should vary with particle size, and indeed such an effect has been observed (Darquenne and Prisk 2003). However, experiments at zero-gravity (Darquenne, West, and Prisk 1998, 1999) indicate that significant dispersion remains, which is independent of particle properties.

It is presently argued that the assumption of alveolar mixing provides an additional, implicit dispersion mechanism, which is active in the acinus and could potentially explain the above zero-gravity dispersion. This mechanism is demonstrated by the toy-example of Figure 2. (The figure is purely schematic and aims at differentiating between conduits with an inlet and an outlet–such as bronchiolar airways and alveolar ducts–on the one hand and volumes with a single opening–such as alveoli–on the other.)

Figure 2. Sketch of particle transport into and out of alveoli during (a) inhalation and (b) exhalation of a bolus.

In order to single-out the contribution of alveolar mixing to axial dispersion, we consider a bolus of particles with radially uniform concentration, which travels along the duct with negligible axial dispersion. During inhalation (Figure 2a), the bolus moves forward and consecutively introduces particles into the expanding alveoli. The bolus reaches a certain penetration depth by the end of inhalation and starts moving back during exhalation (Figure 2b). Now all alveoli start to contract simultaneously and–as a result of mixing–release particles both before and after the passage of the returning bolus. Thus, despite the neglect of axial dispersion in the duct, the exhaled concentration profile is smoothed out by the pattern of particle release from the alveolar volumes. Only if the inflow and outflow of air masses in the alveoli is perfectly ordered and fully reversible, do the inhaled particles exit with the exact timing to reconstruct the original bolus.

2.3. Modeling of airway and alveolar deposition

Deposition probabilities along the airways are calculated by appropriate correlations, which provide the contributions of inertial impaction$P_{I,i}$ (Cai and Yu 1988), gravitational settling $P_{G,i}$ (Heyder and Gebhart 1977) and Brownian diffusion $P_{D,i}$ (Ingham 1975). The total probability of deposition, $P_{T,i},$ in an airway of generation $i$ is given by (14) $$P_{T,i}=P_{I,i}+P_{G,i}+P_{D,i}-P_{I,i}{P}_{G,i}-P_{I,i}{P}_{D,i}-P_{G,i}{P}_{D,i}+P_{I,i}{P}_{G,i}{P}_{D,i}$$(14)

However, the dimensionless numbers that appear in the above correlations are presently calculated at each time step, using the instantaneous values of air flow rate (spatial mean along the airway) and airway diameter. This procedure implies a quasi-steady assumption, i.e., parameters are changing with time slowly enough for the expressions, rigorously valid for steady conditions, to be applicable. Combining the above, the deposition rate per unit length of airway that appears in Equations (11(11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) ) and (12(12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) ) is given by (15) $$L_i=\frac{P_{T,i}{Q}_{i}c}{l_i}$$(15) where $c$ is the local value of concentration and $Q_{i}c$ is the local flow rate of particles through the cross-section at $x.$

Extrathorasic deposition is not presently modeled, instead the value of the deposited fraction is enforced exactly equal to the measurements by Heyder et al. (1986). This choice avoids any bias to the predicted pulmonary deposition resulting from the modeling of extrathoracic deposition. The latter is known to be highly uncertain due to intersubject (and even intrasubject) variability (Grgic et al. 2004).

Alveolar deposition is modeled taking into account the instantaneous volumes of the alveolar cavities. For a settling velocity $u_s=C_C{\rho }_{p}g{d}_p^2/18\mu ,$ the particles that deposit by gravity in a spherical alveolus of diameter $D_a$ during a differential time interval $dt$ are contained in the volume $u_{s}dt(\pi D_a^2/4).$ Thus the deposition rate by gravity inside the alveoli for the entire generation $i$ is, (16) $$n_{G}dt=\left(u_{s}dt\right)\left(\pi D_a^2/4\right)c_{a,i}{N}_{a,i}=\frac{3u_s}{{2D}_a}\left(V_{a,i}{c}_{a,i}\right)dt$$(16) where $N_{a,i}$ is the number of alveoli in generation $i.$

Diffusional deposition is modeled as governed by a boundary layer of thickness $\delta$ close to the wall, across which particle concentration falls from the uniform value $c_{a,i}$ in the bulk of the alveolar volume down to zero (Georgakakou et al. 2016). Therefore, the deposition rate by Brownian diffusion inside the alveoli for the entire generation $i$ is, (17) $$n_{D}dt={\mathcal{D}}_{Br}\frac{c_{a,i}}{\delta }\left(\pi D_a^2\right)N_{a,i}dt={\mathcal{D}}_{Br}\frac{6\Delta }{D_a^2}\left(V_{a,i}{c}_{a,i}\right)dt$$(17) where ${\Delta }^{-1}$ is the dimensionless boundary layer thickness normalized with the alveolar diameter (18) $${\Delta }^{-1}=\frac{\delta }{D_a}$$(18)

In the frame of the alveolar mixing assumption, the thickness $\delta$ is determined by the competition between stirring, which tends to unify particle concentration, and diffusional deposition, which tends to deplete of particles the region close to the wall. Thus, the stronger the mixing, the smaller the value of $\delta .$ In this sense, parameter $\Delta$ may be physically understood as an indicator of the intensity of alveolar mixing.

In Georgakakou et al. (2016), $6\Delta =10$ was chosen in order to fit the benchmark deposition data of Heyder et al. (1986) at the small particle-diameter end. However, there is no a priori reason to assign the same value of Δ for all alveolar generations, and actually the discussion in Section 4.2 argues in favor of an improved choice. Thus in the following we assign as ${\Delta }_i$ the value of $\Delta$ for generation i.

Combining Equations (16)(16) $$n_{G}dt=\left(u_{s}dt\right)\left(\pi D_a^2/4\right)c_{a,i}{N}_{a,i}=\frac{3u_s}{{2D}_a}\left(V_{a,i}{c}_{a,i}\right)dt$$(16) and (17)(17) $$n_{D}dt={\mathcal{D}}_{Br}\frac{c_{a,i}}{\delta }\left(\pi D_a^2\right)N_{a,i}dt={\mathcal{D}}_{Br}\frac{6\Delta }{D_a^2}\left(V_{a,i}{c}_{a,i}\right)dt$$(17) , the deposition rate inside the alveolar volume of generation $i,$ which appears in the mass balance Equations (11(11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) ) and (12)(12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) is given by (19) $$L_{a,i}=\left(\frac{3u_s}{{2D}_a}+\frac{6 {\mathcal{D}}_{Br}{\Delta }_i}{D_a^2}\right)\left(V_{a,i}{c}_{a,i}\right)$$(19)

A fundamental advantage of Equation (19)(19) $$L_{a,i}=\left(\frac{3u_s}{{2D}_a}+\frac{6 {\mathcal{D}}_{Br}{\Delta }_i}{D_a^2}\right)\left(V_{a,i}{c}_{a,i}\right)$$(19) is that it provides the instantaneous deposition rate and is thus appropriate for a differential mass balance like Equations (11(11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) ) and (12(12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) ).

2.4. Method of solution

The system of Equations (11a,b)(11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) and (12a,b)(12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) is augmented by the initial condition that sets all concentrations equal to zero, $c_{a,i}\left(0\right)=c\left(x,0\right)=0,$ and by boundary conditions at the two ends, $x=0$ and $x=L$ of the respiratory tract. More specifically, during inhalation and for continuous particle breathing (20) $$c\left(0,t\right)=1 , \frac{\partial c}{\partial x}\left(L,t\right)=0$$(20)

During exhalation, the conditions imposed at the two ends are (21) $$\frac{\partial c}{\partial x}\left(0,t\right)=\frac{\partial c}{\partial x}\left(L,t\right)=0$$(21)

The above boundary conditions are standard for this problem (Darquenne and Paiva 1994).

The system of Equations (11a(11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) ,b) and (12a(12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) ,b) and boundary conditions Equations (20)(20) $$c\left(0,t\right)=1 , \frac{\partial c}{\partial x}\left(L,t\right)=0$$(20) and (21)(21) $$\frac{\partial c}{\partial x}\left(0,t\right)=\frac{\partial c}{\partial x}\left(L,t\right)=0$$(21) is solved numerically by a finite-volume scheme, using fully implicit time discretization, upwinding for the convection terms and central differences for the diffusion terms (Patankar 1980). Each airway is discretized by N = 10 grid points and inhalation/exhalation is divided into M = 1000 time steps each. This is the choice in Choi and Kim (2007) and our tests also confirm that the results do not change with further grid and time step refinement.

It is noted that the finite-volume method that is applied for the airways conserves mass by definition, a property which is very useful in the present application. However, it is vital that Equations (11b)(11b) $$\frac{d}{dt}\left(V_{a,i}{c}_{a,i}\right)=\frac{Q_{a,i}}{l_{0,i}}\underset{0}{\overset{l_{0,i}}{\int }}\mathit{cdx}+f_i\underset{0}{\overset{l_{0,i}}{\int }}{L}_{i}dx-L_{a,i}$$(11b) and (12b)(12b) $$\frac{d}{dt}\left(V_{a,i}{c}_{a,i}\right)=-Q_{a,i}{c}_{a,i}+f_i\underset{0}{\overset{l_{0,i}}{\int }}{L}_{i}dx-L_{a,i}$$(12b) are also integrated fully implicitly in time, and that the source/sink terms that appear in the two equation sets–and represent exchange of particles between airways and alveoli–are computed consistently, so that mass is also conserved between different lung volumes.

3. Results

The following parameter values are used in the computations. Airway dimensions at $\mathit{FRC}=3 L$ are taken according to the model of Finlay (2001), which considers the 15^th generation as the first of respiratory bronchioles. Following Pollak (1998), the values used in Equation (3)(3) $$\Delta P\left(t\right)=\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{C_{L0}}\mathrm{ }ln\left(\frac{V_{L,\mathit{max}}-V_{L,\mathit{min}}}{V_{L,\mathit{max}}-V_L\left(t\right)}\right)$$(3) are $V_{L,\mathit{max}}=7.5 L, V_{L,\mathit{min}}=1.4 L, C_{L0}=0.0057 L/Pa.$ The fractional coverage of airway walls by alveolar openings is taken $f_{1-14}=0,f_{15}=0.2,f_{16}=0.4,\hspace{0.17em}{f}_{17}=0.7,\hspace{0.17em}{f}_{18-23}=1$ (Tsuda, Henry, and Butler 2008). Each alveolar opening is taken to subtract 20% of the spherical alveolar wall (Weibel 1963), thus $b=0.2$ in Equations (11a)(11a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)-\frac{Q_{a,i}c}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(11a) and (12a)(12a) $$\frac{\partial }{\partial t}\left(A_{i}c\right)=-\frac{\partial }{\partial x}\left(Q_{i}c\right)+\frac{\partial }{\partial x}\left(\mathcal{D}{A}_i\frac{\partial c}{\partial x}\right)+\frac{Q_{a,i}{c}_{a,i}}{l_{0,i}}-L_i+b\frac{L_{a,i}}{l_{0,i}}$$(12a) . Finally, the variation of parameter Δ per generation is taken as: ${\Delta }_{15-23}=(3000,1500,900,400,200,100,30,10,0)$ The physical relevance of the ${\Delta }_i$ values and the sensitivity of the predictions to their parametric variation will be considered in Discussion.

The first set of results are predictions of particle concentration in exhaled air, which are compared to literature data. Choi and Kim (2007) measured the temporal variation of particle concentration in the mouth during expiration. In the experiment, air was inhaled from $\mathit{FRC}=3000 \mathrm{ml}$ up to a tidal volume of $V_T=500 ml,$ with constant flow rate $Q_0=250 ml/s.$ Then, it was exhaled at the same flow rate and down to the reserve volume (RV). The inhaled air was an aerosol containing particles with diameters $d_p=1, 3$ or $5 \mu m.$ Figure 3a is a reproduction of the above data (points) for all three particle sizes, plotted as mouth concentration versus expired air volume. Dotted lines are predictions according to the model of Choi and Kim (2007) and continuous lines are predictions based on the present work. The role of convective mixing in the present predictions is interrogated in Figure 3b by repeating the computation with incorrect values of ${\mathcal{D}}_{\mathit{mix}}.$ Continuous lines are the correct results, whereas dashed lines correspond to ${\mathcal{D}}_{\mathit{mix}}=0$ and dotted lines to ${\mathcal{D}}_{\mathit{mix}}^{\prime }=2{\mathcal{D}}_{\mathit{mix}}.$

Figure 3. Measurements and predictions of particle concentration at the mouth during exhalation. Points are data from Choi and Kim (2007) for 1 μm (○), 3 μm (△) and 5 μm (□) particles. (a) Dashed lines are predictions by a model of the above authors and continuous lines are predictions by the present model. The improvement of the present prediction amounts to 70-80% reduction in the mean relative error with data. (b) The measurements are compared with modified predictions of the present model. Dashed lines correspond to the neglect of axial dispersion and dotted lines to doubling its normal value.

The second set of results are predictions of deposition, which are compared to literature data. Heyder et al. (1986) measured the slow-cleared and fast-cleared fraction of monodisperse, unit-density spheres with diameters in the range $0.05\mathrm{\mu }\mathrm{m}⩽d_p⩽15\mathrm{\mu }\mathrm{m}.$ They used three different breathing patterns that may schematically be classified as shallow: $V_T=500 {\mathrm{cm}}^3,Q_o=250{\mathrm{cm}}^3/s,$ intermediate: $V_T=1000 {\mathrm{cm}}^3, Q_o=250{\mathrm{cm}}^3/s$ and deep: $V_T=1500 {\mathrm{cm}}^3, Q_o=750{\mathrm{cm}}^3/\mathrm{s}.$ Their data may be converted to pulmonary (P) and tracheo-bronchial (TB) deposition following a recommendation by ICRP (1994) that TB deposition also has a slow-cleared fraction, ${\phi }_{sl},$ which is estimated as ${\phi }_{sl}=0.5\hspace{0.17em}\text{exp}[-0.63(d_p-2.5)]$ for $d_p⩾2.5 \mathrm{\mu }\mathrm{m}$ and ${\phi }_{sl}=0.5$ for $d_p<2.5 \mathrm{\mu }\mathrm{m}.$

Predictions of deposition per generation for representative particle sizes are shown in Figures 4a–c, corresponding respectively to the shallow, intermediate and deep breathing pattern. However, these predictions cannot be directly compared with the aforementioned measurements, as the latter offer only cummulative results.

Figure 4. Predictions of the deposited fraction per generation for various particle sizes, and for the (a) shallow, (b) intermediate, and (c) deep breathing pattern.

Thus, predictions are summed up to give the slow-cleared fraction, i.e., the pulmonary and the appropriate part of the tracheo-bronchial deposition, and are compared to data in Figures 5a–c (corresponding again to the shallow, intermediate and deep breathing pattern). Circles are the measurements by Heyder et al. (1986) and continuous lines the results of the present model. The other lines are predictions by the simple filtration model of Georgakakou et al. (2016). Dashed lines are the original results in their article, which were based on a uniform value of parameter ${\Delta }_i=10/6.$ Dotted lines are the updated predictions of that model, using the present values of parameter ${\Delta }_i.$

Figure 5. Pulmonary deposition during (a) shallow, (b) intermediate, and (c) deep breathing pattern. Measurements by Heyder et al. (1986) (points) are compared to predictions by the present model (bold lines), by the original model of Georgakakou et al. (2016) (dashed lines) and by the model of Georgakakou et al. (2016) using the presently selected values of ${\Delta }_i$ (dotted lines). The present prediction reduces the mean error by 42% in (a), 65% in (b), and 87% in (c) compared to the original model of Georgakakou et al. Use of the latter with the present ${\Delta }_i$ leads to large overprediction.

4. Discussion

4.1. Trends in particle trapping and deposition

Attention is first drawn to the predictions of particle concentration in exhaled air. It is evident from Figure 3a that, beyond the first $200 ml$ of expiration, the model is in good agreement with the data. The relative error, $\sum \mathit{abs}\left[\left(c_{\mathit{pred}}-c_{\mathit{dat}}\right)/c_{\mathit{dat}}\right],$ is for the present model 20–30% of the error by the model of Choi and Kim (2007). The temporal variation of concentration in the exhaled air–and in particular the longtime tail–is expected to reflect the rate of release of particles trapped in the alveoli. Thus, the observed agreement argues in favor of the assumption of alveolar mixing, on which the present model is based. However, it is admitted that other factors–such as ventilation distribution between dependent and non-dependent regions of the lung and airway branching asymmetry–are likely to also affect exhaled concentration profiles, and thus the quality of the agreement may to some extent be fortuitous.

It is noted that the above agreement does not depend on the choice of the only adjustable parameters in our model, i.e., the dimensionless boundary layer thicknesses in the alveoli, ${\Delta }_i^{-1},$ as this choice affects only the small-diameter end of the particle spectrum $\left(d_p<1 \mathrm{\mu }\mathrm{m}\right).$ On the contrary, the prediction depends on the correct modeling of axial dispersion, as is indicated by Figure 3b. Differences observed with the erroneous selections of ${\mathcal{D}}_{\mathit{mix}}$ are more pronounced for $d_p=1 \mathrm{\mu }\mathrm{m},$ small for $d_p=3 \mathrm{\mu }\mathrm{m}$ and negligible for $d_p=5 \mathrm{\mu }\mathrm{m}.$ As particles with $d_p=1 \mathrm{\mu }\mathrm{m}$ deposit the least, this behavior is attributed to the larger concentrations that float along the airways.

The importance of the inclusion of airway deformation in the model was interrogated by additional computations, holding airway dimensions constant at FRC values. In all cases, the exhaled particle concentrations decreased, with maximum deviations 15–30% around the middle of the exhalation maneuver and reduction of total exhaled particles by 5–8% (deviations increased with particle size). This behavior is understood by considering that, with rigid airways, the same inhaled volume results in larger alveolar expansion. Thus, more particles are trapped and fewer are exhaled.

An interesting deduction from the model is that the exhaled concentration profiles depend very weakly on the extent of alveolar deposition. Indeed, setting the latter artificially equal to zero (which amounts to roughly 25% reduction in total deposion) increases total exhaled particles by less than 1%. This is a confirmation of similar remarks in the literature (Darquenne and Paiva 1994), cautioning against the use of exhaled concentration profiles as indicators of deposition, because of the significant fraction of particles that may remain trapped in the alveoli.

Next, attention is focused on predictions of deposition and their comparison with data in Figures 5a–c. The results of the present model (bold lines) are an improvement over those of the Lagrangian model of Georgakakou et al. (dashed lines), particularly at the small and intermediate diameter range $\left(d_p<1 \mathrm{\mu }\mathrm{m}\right),$ which is affected by the choice of parameters ${\Delta }_i.$ In this range, the mean error between predictions and data is for the present model 58%, 35% and 13% of that of the older one, for the shallow, intermediate and deep breathing pattern respectively. The improvement is also reflected in the accurate prediction of the slope of the deposition curve at the fine-particle end and of the curvature around the minimum.

The dotted lines in Figures 5a–c are derived from the Lagrangian model of Georgakakou et al. (2016) using the present distribution of ${\Delta }_i$ values. The strong over-prediction of the data observed is presently understood as a consequence of the neglect of airway deformation during the breathing maneuver. If airway diameters are held constant, then only the alveoli expand during inhalation, resulting in significantly higher air and particle inflow. Once the number of particles that enter the alveolar volumes is overpredicted, agreement with data necessitates artificially low deposition rates. This explains the small value of $\Delta$ found by Georgakakou et al. (2016) to fit best the data. It is concluded that accurate prediction of the air volumes that actually reach the alveoli is of vital importance for the quality of the model.

Finally, a note is in order for the apparently counter-intuitive outcome (both in the data and the model) that the intermediate breathing pattern results in higher pulmonary deposition than the deep breathing one. This behavior is explained if it is noted that the former pattern has twice the inhalation time of the latter, thus affording more time for gravitational settling in the acinus. Also, the higher velocities during deep breathing, result in stronger inertial deposition along the upper airways and thus decreased particle availability deep in the lung.

4.2. Implications of the assumption of perfect mixing

The central assumption of the present model, i.e., that the particle concentrations in the alveolar volumes of each generation may be treated as spatially uniform, has direct consequences on both particle transport and deposition. With respect to transport, it dictates that particles that have entered the alveoli are released during exhalation at a very slow rate, and thus remain trapped for long times even if they do not deposit. This behavior (kinematic irreversibility) is contrasted to the reversible movement of air in and out of the alveoli, which would result in the removal during exhalation of all the particles that had entered the alveoli and did not deposit on their wall.

It is also useful for the discussion that follows to recall that, according to the present model, particles move between the airways and the alveolar volumes of each generation by two mechanisms whose contribution is additive, (a) by convection with the air flow entering or leaving the alveoli and (b) by gravitational settling and/or Brownian diffusion on the fraction of the airway wall that is absent due to the alveolar openings. (This summation is equivalent to the classical formulation of mass flux of a component in a mixture as the sum of a convective and a molecular contribution.) In particular, particles enter the alveoli by mechanism (a) during inhalation and by mechanism (b) during the entire breathing cycle.

With respect to alveolar deposition, the assumption of spatially uniform particle concentration per alveolar generation is the basis of the simple expression, Equation (19)(19) $$L_{a,i}=\left(\frac{3u_s}{{2D}_a}+\frac{6 {\mathcal{D}}_{Br}{\Delta }_i}{D_a^2}\right)\left(V_{a,i}{c}_{a,i}\right)$$(19) . However, the significance of the assumption is different for each of the two pulmonary deposition mechanisms, i.e., gravitational settling and Brownian diffusion.

More specifically, gravitational deposition is correctly predicted by Equation (16)(16) $$n_{G}dt=\left(u_{s}dt\right)\left(\pi D_a^2/4\right)c_{a,i}{N}_{a,i}=\frac{3u_s}{{2D}_a}\left(V_{a,i}{c}_{a,i}\right)dt$$(16) , even when the particle distribution is non-uniform. For example, a higher density close to the wall would result in larger or smaller deposition depending on the orientation of the alveolus with respect to gravity, but the overall effect of a large number of randomly oriented alveoli would be equivalent to that resulting from uniform particle distribution. Thus, it is argued that gravitational deposition inside the alveoli is influenced by the rate of particle entrance and exit but is largely insensitive to the detailed spatial distribution.

It is further noted that, for particles with $d_p\ge 2\text{ μm},$ the characteristic gravitational time in Equation (16)(16) $$n_{G}dt=\left(u_{s}dt\right)\left(\pi D_a^2/4\right)c_{a,i}{N}_{a,i}=\frac{3u_s}{{2D}_a}\left(V_{a,i}{c}_{a,i}\right)dt$$(16) , ${\tau }_{gr}=\left({2D}_a/3u_s\right),$ is less than one second, and thus all the particles that enter the alveoli during inhalation (and do not move out during exhalation) have enough time to deposit. Also, and because of the high deposition probabilities of large particles in the proximal airways, the further influx by mechanism (b) is not significant. Thus, the prediction of gravitational deposition in the alveoli appears relatively straightforward. In particular, it depends to a large extent on the correct prediction of tracheo-bronchial deposition, as the latter determines the number of particles that actually reach the alveolated airways.

Prediction of the alveolar deposition for particles affected by Brownian diffusion is more uncertain because it depends on the extent of alveolar mixing per generation (which in the present model is expressed by parameters ${\Delta }_i$). Small particles in the considered range ($0.05 \mathrm{\mu }\mathrm{m}\le d_p\le 1\text{ μm}$) have two important characteristics: (i) They deposit far less in the tracheo-bronchial section, so they are contained in higher concentrations in the air that reaches the acinus. (ii) Their deposition probability in the distal airways increases strongly with generation. Thus, the fraction of particles that enter alveoli by deposition on the wall openings (mechanism (b) that is operative during the entire breathing cycle) grows rapidly in importance distally.

In order to understand the deposition behavior of these small particles, we consider the timescale for depletion of the alveolus by deposition due to Brownian diffusion, ${\tau }_{\mathit{diff}}={D_a}^2/\left(6{\Delta }_i{\mathcal{D}}_{Br}\right),$ which is derived from Equation (17)(17) $$n_{D}dt={\mathcal{D}}_{Br}\frac{c_{a,i}}{\delta }\left(\pi D_a^2\right)N_{a,i}dt={\mathcal{D}}_{Br}\frac{6\Delta }{D_a^2}\left(V_{a,i}{c}_{a,i}\right)dt$$(17) . Values of ${\Delta }_i\sim {10}^3$ lead to ${\tau }_{\mathit{diff}}\sim 0.2 \mathrm{s},$ indicating that particles which enter the alveoli have enough time to deposit almost completely and the actual value chosen for ${\Delta }_i$ is of limited importance. This appears to be the case for the first three alveolated generations, where particle inflow by inhalation (mechanism (a)) is dominant.

However, prediction of deposition for the next generations becomes increasingly more sensitive to the value of ${\Delta }_i.$ This is understood in terms of the increasing influence of particle inflow to alveoli through deposition on the openings of the airway wall (mechanism (b) that operates during the entire breathing cycle). The particles that enter during exhalation have less time to deposit and their fate depends more strongly on the intensity of mixing. In particular, no mixing occurs in the last generation, and thus deposition of small particles there is expected to be negligibly small.

Selection of appropriate values for ${\Delta }_i$ was presently done on a purely empirical basis, i.e., by fitting predictions to data of pulmonary deposition. However, there appears to be direct correlation between ${\Delta }_i$ and parametric conditions that lead to chaotic mixing. More specifically, the original literature (Tsuda, Henry, and Butler 1995, Tsuda et al. 2002, 2008) indicates that the onset of chaotic mixing requires large values of the ratio $Q_i/Q_{a,i}.$ Using the estimates provided by Tsuda, Henry, and Butler (2008), the order of magnitude of this ratio from the 15^th to the 23^rd generation is (3000, 1500, 900, 400, 200, 100, 40, 15, 3), to be compared with the above discussed values of ${\Delta }_i.$ Though there is no fundamental justification for equating the dimensionless boundary layer thickness to the ratio $Q_{a,i}/Q_i,$ the similarity of scales is notable.

The present model of diffusional deposition in the alveoli is based on a double role of Brownian diffusivity. First, it defines the lengthscale, $x_{\mathit{diff}},$ below which diffusion smears out the closely spaced concentration gradients that are created in the alveoli by the stretching and folding deformations of chaotic flow (Tsuda, Laine-Pearson, and Hydon 2011), leading to a uniform concentration in the core of the alveolus. This lengthscale may be estimated as $x_{\mathit{diff}}\sim \sqrt{{\mathcal{D}}_{Br}t},$ which–for a particle of $0.1–1 \mathrm{\mu }\mathrm{m}$ and typical duration of a breathing cycle–amounts to $x_{\mathit{diff}}\sim 10–30 \mathrm{\mu }\mathrm{m}.$ Second, Brownian diffusivity determines deposition as a boundary-layer effect from this idealized, uniform concentration in the bulk, toward the wall where particle concentration is by definition zero (any particle that touches the wall deposits instantaneously).

4.3. Implications of the modeling of lung morphology and its temporal variation

Modern morphometric studies, based on the characterization of microscopic sections by stereology (Ochs et al. 2004; Weibel, Hsia, and Ochs 2007; Vasilescu et al. 2012), have provided valuable new information on the detailed lung micro-structure and in particular the topology of alveolar walls and openings. These findings support the simplified model used in the present work (which consists of open-ended ducts in one-dimensional flow, surrounded by uniform peripheral chambers in good mixing) in two respects: First, they confirm that the alveolar entrance rings, formed by the strong fiber tracts in the free edges of the alveolar septa, constitute a virtual duct. Second, they indicate that the size of a single alveolus varies little around a mean of $4.2x{10}^6 {\text{μm}}^3\mathrm{ }(D_a\approx 200 \mathrm{\mu }\mathrm{m}),$ and it is mainly the number of alveoli that varies with different lung volumes.

Airway elasticity has presently been taken into account, motivated by the need to estimate correctly the flow rate of air that actually reaches the alveolar volumes. The elasticity model employed is based on intramural pressure, and extrapolates measurements taken in the central and conducting airways of excised human lungs (Lambert et al. 1982) in order to estimate the expansion of alveolated ducts. There are other, equally legitimate modeling choices, such as for example isometric volume scaling (Choi and Kim 2007).

However, the reasonable criticism–that the extrapolation of elasticity characteristics from the conducting airways to the acinus is open to question–raises a more fundamental issue: Namely, what is the contribution of ducts and alveoli to volume change during expansion and contraction of the acinus, and how do the alveoli change in shape. Potential answers were discussed long ago (Greaves, Hildebrandt, and Hoppin 1986), and the question was highlighted more recently as “… may be the last frontier of gross anatomy” (Smaldone and Mitzner 2012).

Though the above issue is clearly beyond the scope of the present work, it is useful to clarify that our model corresponds to isotropic expansion and contraction of alveoli, i.e., it is valid for internal parenchymal surface area that varies with lung volume raised to the 0.67 power (For isotropic expansion, $S\sim d^2, V\sim d^3\to S\sim V^{2/3}$). Isotropic variation without change in the number of active alveoli (i.e., no recruitment) has been questioned by measurements of the apparent diffusivity of ³He and of the acinar deposition of monodispersed aerosols during breathhold (Smaldone, Mitzner, and Itoh 1983; Hajari et al. 2012), but is supported by data from postmortem anatomy of frozen lungs (Glazier et al. 1967; Soutiere and Mitzner 2004).

5. Conclusions

A dynamic, single-path model has been developed for the prediction of particle transport and deposition deep in the lung. The model takes into account convection and dispersion along the air ducts, as well as particle exchange between ducts and alveoli. In particular, it models alveolar volumes of each generation as independent mixing chambers. In this respect, it is demonstrated that the assumption of alveolar mixing introduces an implicit axial dispersion mechanism that is unique to the acinus and independent of particle size and gravity.

Air inflow to (and outflow from) alveoli is calculated by accounting for the deformation of airways during breathing. Mechanistic models are derived for deposition inside the alveoli. In particular, deposition by Brownian diffusion is modeled by a concentration boundary layer, whose dimensionless thickness, ${{\Delta }_i}^{-1},$ varies inversely with the intensity of mixing.

The values of ${\Delta }_i$ per generation were tailored so as to optimize the prediction of benchmark data of pulmonary deposition (Heyder et al. 1986). The agreement is satisfactory for all breathing patterns tested and extends over the entire particle diameter range measured ($0.05 \mathrm{\mu }\mathrm{m}-15 \mathrm{\mu }\mathrm{m}$). It is notable that the selected values of ${\Delta }_i$ are found to correlate with the ratio, $Q_i/Q_{a,i},$ of airway flow to alveolar inflow. High values of this ratio are known to result in chaotic mixing (Tsuda, Henry, and Butler 1995, Tsuda et al. 2002).

The predictions of the present model are found to reproduce independent measurements of particle concentration in the mouth during pronounced exhalation of a continuously inhaled aerosol (Choi and Kim 2007). Because such measurements are expected to depend strongly on the release rate of particles retained in the acinus, it is argued that the agreement lends support to the assumption of alveolar mixing that is central to the model.

Funding

Partial financial support of this research by a donation from ELPEN Pharmaceutical is gratefully acknowledged.