Environ. Sci. Technol. 2001, 35, 3727-3734
Integrated Exposure and Dose Modeling and Analysis System. 3. Deposition of Inhaled Particles in the Human Respiratory Tract M I H A L I S L A Z A R I D I S , †,‡ D A V I D M . B R O D A Y , †,§ ØYSTEIN HOV,‡ AND P A N O S G . G E O R G O P O U L O S * ,† Environmental and Occupational Health Sciences Institute, Rutgers University and UMDNJ - R.W. Johnson Medical School, 170 Frelinghuysen Road, Piscataway, New Jersey 08854, Norwegian Institute for Air Research, P.O. Box 100, Instituttveien 18, N-2007 Kjeller, Norway, and Agricultural Engineering Technion - Israel Institute of Technology, Haifa 32000, Israel
Detailed information on the composition-resolved size distribution of particulate matter deposited along the human respiratory tract can help linking epidemiological, toxicological, and pathological studies and thus potentially improve the understanding of the origin of pulmonary disorders induced by respirable pathogens. For this purpose, a new mechanistic dosimetry model describing the dynamics of respirable particles in the human airways was developed. Model predictions of transport and fate of inhaled aerosols are based on solutions of the aerosol general dynamic equation, which describes changes in particle size and mass distributions resulting from processes such as nucleation, condensation, coagulation, gas phase chemical reaction, and deposition. To compensate for approximating the three-dimensional problem by considering only axial variations along the airways, boundary layer effects are introduced via appropriate dimensionless transport parameters. The architecture of the human lung is described by Weibel’s simple regular dichotomous model. An important advantage of the present approach is that it allows testing the significance of intersubject lung morphology and ventilation variability for particle deposition and dose calculations. The model predicts the evolution of size and composition distributions of inhaled particles and the deposition profile along the human lower respiratory tract: in general, model predictions are in qualitative and quantitative agreement with tracheobronchial and alveolar deposition data.
Introduction A conceptual and theoretical framework for a modular integrated exposure and dose modeling and analysis system (EDMAS) has been described (1), implemented, and tested for (i) multipathway/multiroute exposures to benzene vola* Corresponding author phone: (732)445-0159; fax: (732)445-0915; e-mail:
[email protected]. † Rutgers University and UMDNJ. ‡ Norwegian Institute for Air Research. § Agricultural Engineering Technion - Israel Institute of Technology. 10.1021/es001545w CCC: $20.00 Published on Web 08/11/2001
2001 American Chemical Society
tilized from tap water (2) and (ii) inhalation of ozone (3). The EDMAS approach was further extended to predict indoor air concentration of gaseous pollutants and associated fine particulate matter (PM) due to temporal variations in ambient pollutant concentrations and as a result of indoor sources and sinks (4). In particular, PM formation and removal indoors, resulting from direct emissions (indoor sources, particle resuspension), deposition, and physicochemical transformation processes, were considered. Apart from transport and mixing, processes accounted for include heterogeneous (reaction on surfaces) and homogeneous (formation of secondary organic and inorganic particles, nucleation) aerosol and gas-phase chemistry, condensation, and coagulation. The present dosimetry model represents the third contribution to the integrated paradigm, where lung deposition of inhaled PM with specific size-resolved chemical characteristics is predicted. Combined with the previous work, this enables study of different “real-world” exposuredose scenarios for PM. The relationship between exposure to airborne PM and resulting health effects is the subject of an ever-increasing number of studies (5-11). While epidemiological studies suggest an association between ambient PM concentrations and increased human morbidity and mortality (12, 13), toxicological studies have recently begun to provide potential biological explanations for this observed association (14). Assessing health risks from exposure to particulate matter includes consideration of particle-specific mechanisms of chemical disposition, toxicant-target interactions, and tissue response. Including such factors as regional variability in deposition and clearance rates by gender, age, and state of health (14, 15) makes exposure assessment studies very complex endeavors. The development of mechanistic dosimetry models is an important step in understanding exposure-dose-response relationships for ambient PM and can certainly help in assessing human health risks from inhalation of different toxicants in occupational and residential environments. Such mechanistic models can utilize data reflecting intra- and intersubject biological variability by employing statistical distributions for model parameters, enabling prediction of the biological dose distribution in target individuals and specific subpopulations. Numerous mathematical models for predicting particle deposition in the human airways have been developed over the years. In general, available practical (i.e. “routinely operational”) models can be classified into three major groups. The first group follows an approach based on the concept of applying compartmental modeling to the human lung anatomy. Numerous compartmental models have been proposed, differing in the representation of the tracheobronchial tree, the breathing physiology and the resulting air flow, and the expressions used for calculating particle deposition efficiencies (16-18). In the second approach, the human lung is modeled as a chamber shaped like a trumpet with a variable crosssectional area (19-22). This approach is not intended to provide a realistic morphometric description of the lung but rather to enable efficient implementation of one-dimensional numerical schemes for solving the spatiotemporal concentration field. The third approach uses a combination of theoretical and empirical expressions to predict particle deposition in the human airways (23). Such is the approach taken by International Commission on Radiological Protection ICRP66 recommended model (15). Since this model is inherently VOL. 35, NO. 18, 2001 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
3727
semiempirical, it has a limited applicability to nonstandard cases that are not included in the database from which the empirical expressions were derived. Indeed, the ICRP66 has a questionable predicting power when aerosol particles cannot be considered inert and stable, i.e., for particles that do not preserve their composition and size during the calculation. Moreover, ignoring the particle physicochemical properties and referring to generic PM is inconsistent with current atmospheric and indoor aerosol transport and fate models, where chemical reaction and particle growth result in time-dependent size and composition distributions. In fact, the lack of physicochemical attributes of aerosol particles found in outdoor, indoor, or occupational environments is the major shortcoming of most previous attempts to predict particle transport and fate in the human respiratory tract. For example, not accounting for properties such as solubility and hygroscopicity may have a profound effect on the predicted deposition profile (24). The present model was developed in order to overcome certain limitations of existing lung deposition models and to provide a more robust quantitative description of particle dynamics in the human respiratory tract while at the same time maintaining a simple description of respiratory physiology. The model is in fact one-dimensional and as such only simplified velocity profiles are considered. For small radial temperature gradient, e.g. when temperature difference e10 °C prevails between the walls and the air, a onedimensional flow within a tubular enclosure is a valid model with tolerable errors (25-27). Nevertheless, boundary-layer type analysis was implemented to improve the predicted size-distribution profiles of both the persistent and deposited aerosol particles. Thus, this model transcends previous models in which (a) boundary layer effects were not considered (20-22), (b) empirical expressions were used for describing particle transport and deposition, and (c) gasphase processes were not included (17, 28). The mechanistic approach and the flexible nature of the present model are applied to study the relative importance of different physicochemical processes and the effect of different lung morphologies and breathing conditions on aerosol dynamics and deposition in the human airways. Predicted deposition profiles are generally in agreement with available experimental data.
is expressed in terms of the total deposition velocity of particles in the size class k, vd,k, the differential volume element ∆V, the effective deposition area Ad (the projection of ∆V normal to the flow direction), and the air velocity u. For the simplified plug flow considered here, the left-hand side (LHS) of eq 1 can be written as
Aerosol Dynamics in the Human Respiratory Tract
where Fa ) ∑jcjMw,j is the air density. The LHS of eq 5 equals the sum over all gaseous species depletion sinks and formation sources. For steady-state conditions, the first term on the RHS of eq 5 is zero. Furthermore, since airway bifurcations are treated as discontinuities in both the crosssectional area and the air velocity, eq 5 reads
Different simultaneous processes control the transport and fate of particles and gases in the human respiratory tract. For example, size changes due to both condensation and agglomeration can significantly affect the dynamics and deposition of inhaled particles (14). Likewise, at high enough concentrations of vapor precursors secondary particle formation may occur in the upper airways, leading to changes in the particle size and composition distributions. Assuming steady-state conditions (29), the aerosol general dynamic equation (GDE) in a one-dimensional form along the flow direction is
( )
dnk 1 dnk ) Jkδ(k - k*) + dx u dx
coag
+
( ) dnk dx
grow
-
vd,kAd n u∆V k (1)
where nk is the number concentration of particles in the size class k (having a radius rk). The first term on the right-hand side (RHS) of eq 1 corresponds to particle formation by nucleation, where Jk is the nucleation rate and Dirac’s delta function, δ(k - k*), indicates the size of the nucleating droplet. The other terms on the RHS of eq 1 represent size changes due to particle coagulation, particle growth by condensation and chemical reactions, and particle removal by deposition, respectively. In particular, the differential particle removal 3728
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 35, NO. 18, 2001
[
dnk 1 ∂nk 1 ∂ ) + (Aunk) dx u ∂t A ∂x
]
(2)
where time is related to axial position of the aerosol bolus via the air velocity u. The mass distribution of different species is calculated by solving the condensed phase species equation
[
]
d(Fp,iΥij) d(Fp,iΥij) ) dx dx
+
coag
[
]
d(Fp,iΥij) dx
gtp
-
vd,iAd (F Υ ) u∆V i ij (3)
where Fp,i and Υij are the particle average density and the fraction of species j at the ith grid point (bin), respectively. The first term on the RHS of eq 3 corresponds to particle coagulation, the second term refers to gas-to-particle transformation due to nucleation and condensation growth, and the last term represents composition-resolved deposition. The rate of change of the molar concentration of the jth gaseous phase species, cj, is calculated from the gas-phase species conservation equation
[ ]
dcj dcj ) dx dx
-
form
[ ] dcj dx
-
ncr
vd,jAd c u∆V j
(4)
The first term on the RHS of eq 4, the species formation rate, can be calculated from reaction kinetics or from local gasphase chemical equilibrium; the second term stands for the rate of vapor depletion due to nucleation, condensation and chemical reactions; and the last term represents depletion due to vapor deposition. The rate of change of the air velocity is described by the continuity equation for an incompressible fluid
dFa ∂Fa 1 ∂(AuFa) ) + dt ∂t A ∂x
d(uFa) dx
)
∑m˘
ij
(5)
(6)
i
where m ˘ ij is the depletion or formation rate of the jth species in the ith bin. Brief details on the processes that are included in the model are given in the following sections. Homogeneous Binary Nucleation. In general, particle formation in the human respiratory tract is not probable for typical environmental concentrations of precursor gaseous species. However, particle formation by nucleation may occur in specific exposure scenarios when ambient concentrations of vapor species, such as sulfuric acid and ammonia, are much higher than the background values. Such conditions may occur in occupational environments of specific industries (i.e. petrochemical, power generation, glass, metal) and, perhaps, in urban areas with adverse dispersion conditions and high emissions of combustion products. Particle formation in the human respiratory tract is modeled by either the classical nucleation theory for a two-
component (binary) homogeneous aerosol or by a new steady-state model in which energy fluctuations trigger new particle formation (30). In the classical theory, the nucleation rate is expressed as
(
J ) Joexp -
F* kBT
)
(7)
where Jo is a kinetic prefactor [particles/m3s], F*/kBT is the normalized free energy of formation of the critical droplet, and kB is Boltzmann’s coefficient. The kinetic prefactor for binary nucleation may be written in a form similar to the one used in unary nucleation (31)
Jo ) FvApBZ
(8)
where Fv is the combined density of the condensable vapors (Fv ) F1v + F2v), Ap is the surface area of the droplet, B is the average growth rate, and Z is Zeldovich’s nonequilibrium factor. The second approach for calculating the nucleation flux is based on a solution of the appropriate Fokker-Planck equation, where the effect of energy fluctuations is included (30). This model predicts much higher nucleation rates than the classical theory and is used as a conservative estimate of nucleation. Condensation. Fine particles, either inhaled or formed within the human respiratory tract, can grow due to vapor condensation or absorption. Steady-state solutions for mass and heat diffusion toward a single droplet in an infinite fluid were obtained from first-order phenomenological transport equations for a binary-mixture (32). These equations, derived for the continuum regime, were modified to apply for particles in the transition regime (33), corresponding to Knudsen numbers Kn ) λ/rp∼O (1), where λ is the molecular mean free path and rp is the particle radius. For the near-saturation physiological conditions found in the lungs, accurate calculation of the growth rate requires excessive computation time. To speed the calculation, the growth process is modeled using the well-known Mason eq 34 with a zeroth-order approximation for the mass balance equation. Evaluation of the vapor pressure at the particle surface is obtained via a linearized Clausius-Clapeyron equation, using the first term of its series expansion and assuming that TpT∞ ≈ T∞2, where the subscripts p and ∞ refer to the particle surface and the bulk fluid, respectively. The error in using this approximation is small when (Tp - T∞)/T∞ , 1, as is the case in the human lung. Using the above approximations, a modified expression of Mason’s equation is obtained for the mass flux toward the particle (33, 35)
dm
)
dt
∑
4πrp(S∞, j - Sr,j) NM,j
j
+
βM
(9)
NT,j βT
where
NM,j )
R gT ∞
, NT,j )
D∞,jMw,jpsj (T∞)
(
Lj LjMw,j -1 kT∞ RgT∞
)
(10)
βM and βT are correction factors for mass and heat transfer toward a particle of size corresponding to the transitional regime, D∞,j is the binary diffusion coefficient of the jth species, Mw,j and Lj are the vapor molecular weight and latent heat of condensation at T∞, respectively, and S∞,j and Sr,j are the saturation ratios far from the particle and on its surface, respectively. In this formulation only insoluble (yet wettable)
aerosol particles are considered. The initial mole fraction of the liquid film on the particles is calculated utilizing the heterogeneous binary nucleation theory (36). Coagulation. Coagulation results in changes in the particle size distribution and a net loss of one particle per coagulation, although the total mass is preserved. The general equation controlling coagulation is the Smoluchowski equation
dnk dt
) 1/2
∑
i+j)k
∞
K(mi, mj)ninj - nk
∑K(m , m )n i
k
(11)
i
i)1
where ni is the number of particles in the ith size bin (of a radius rp,i), mi is the mass of all the particles in the ith bin (37), and K(mi, mj) is the kernel for coagulation of mass mi and mj particles. The first term on the LHS of eq 11 represents summation over those particles for which mi + mj ) mk, whereas the second term represents removal of particles from the kth bin due to size increase resulting from coagulation. The coagulation kernel depends on the mechanisms causing particles to approach each other and collide. Here, coagulation processes due to Brownian diffusion, gravitational settling, and turbulence are considered (26). Brownian coagulation (with a kernel KB) is modeled following the work of Wagner and Kerker (38), based on matching particle fluxes in the continuum and the free molecular regimes. As a result of the laryngeal jet that is induced by the glottis aperture, turbulent coagulation of particles that are smaller than the characteristic length of the turbulent eddies is calculated in the first few generations of the tracheobronchial tree, following the work of Saffman and Turner (39). Turbulent coagulation comprises the mechanisms of inertial agglomeration and relative particle motion, both arising from microscale velocity gradients in the isotropically assumed turbulent field. The turbulent coagulation kernel is therefore a sum of the inertial turbulent kernel KI and the turbulent agglomeration kernel KS (40). Gravitational agglomeration (with a kernel KG) occurs due to a relative settling velocity between particles of different size. The total coagulation rate is obtained by combining the kernels of the different coagulation mechanisms, with the turbulent and gravitational contributions added quadratically (39)
Ktotal ) KB + (K2G + K2S + K2I )1/2
(12)
Particle Deposition. Inhaled particles and gases deposit in the human airways by different deposition mechanisms, depending to varying extent on particle size, local flow field, and lung architecture. Calculation of particle deposition is required to evaluate the size resolved loss term in the aerosol general dynamic eq 1. Thermophoresis, Brownian diffusion, inertial impaction, and gravitational settling are the particle deposition mechanisms addressed in this study. Gases deposit as a result of vapor condensation on particles and airway surfaces. An option of including heterogeneous chemical reactions of gaseous deposits on the airway surfaces is also provided. This module will be tested in a future work. The main features of the deposition mechanisms included in the model are discussed below. Thermophoresis. The body core temperature is independent of the ambient temperature, the extreme values of which in certain climate zones may differ significantly from 37 °C. During breathing, air passes proximate to biological heat exchangers in the extrathoracic and upper airways, causing conditioning of its temperature and humidity. Radial temperature gradient may induce thermophoretic force on small aerosol particles, since gas molecules impact the particle with a momentum proportional to the local (differential) temperature. In general, thermophoretic force tends to move fine aerosol particles from high to low-temperature VOL. 35, NO. 18, 2001 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
3729
regions, i.e., enhancing deposition when the inhaled air is at temperature higher than body core temperature. The thermophoretic particle flux to the airway walls is expressed as
Jc ) ngυth
(13)
where ng is the bulk particle concentration and υth is the thermophoretic velocity. A detailed overview of expressions for υth in different flow regimes appears in ref 27. The temperature gradient is estimated here in terms of the width of the thermal boundary layer adjacent to the airway walls
δth ) da/Nu
(14)
where da is the airway diameter, Nu ) dah/kg is the Nusselt number, h is the heat transfer coefficient, and kg is the air thermal conductivity. Typically, the flow in the lower airways of the lung is characterized as undeveloped laminar flow (29), i.e., it is not fully developed along a large portion of the airway and, hence, entrance effects are important. These effects are incorporated here by defining a position-dependent Nusselt number whose asymptotic value is 3.66, in accordance with the solution for the fully developed flow. Laminar Diffusion. The diffusional deposition velocity in laminar flow, υl, is calculated in the model in terms of the thermophoretic deposition velocity, υth (27)
υl )
υth exp(υthδp/Dp) - 1
(15)
where Dp, the particle diffusion coefficient, is calculated from Stokes-Einstein relation (41)
(16)
Cc is Cunnigham’s correction factor, and µ is the viscosity of air. The particle boundary layer thickness, δp, is estimated from an analogy between particle and thermal transport equations (40)
( )
Dp 1.73 δ 2 th,p ν
1/3
(17)
where δth,p is the particle thermal boundary layer. The width of the latter is estimated using the expression for calculating the width of thermal boundary layer adjacent to the airway walls, eq 14. However, because particles larger than ∼0.5 µm are characterized by a Schmidt number >1, i.e., the particle boundary layer is contained within the particle thermal boundary layer, the asymptotic value of the pertinent Nusselt number is 4. In cases where (quasi-) isothermal conditions pertain and therefore the thermophoretic velocity is negligible, the laminar deposition velocity (in the limit υth f 0) is Dp/δp. Turbulent Diffusion. In the first few airway generations the increased particle diffusivity induced by the prevailing turbulent flow enhances particle deposition on the airway walls. Diffusional deposition from turbulent flows is calculated by a simplified eddy deposition model, assuming that particles are transported by turbulent fluctuations (coherent structures) from the edge of the hydrodynamic boundary layer to a layer extending to a stopping distance away from the airway walls. The particles complete their motion by diffusing across the laminar viscous sublayer of the hydrodynamic boundary layer (40). Like diffusion, thermophoresis is considered only within the laminar viscous sublayer. 3730
9
Fv,j )
vv,jMw,j v [pj - psj (T)] RT
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 35, NO. 18, 2001
(18)
where vv,j is the vapor deposition velocity of species j, pvj and psj (T) are the partial vapor pressure of species j and its saturation pressure at temperature T, respectively, R is the universal gas constant, and Mw,j the molecular weight of species j. The vapor deposition velocity vv,j can be approximated by (44)
vv,j ) Dj/δc
kBTCc(Kn) Dp ) 6πrpµ
δp )
Gravitational Settling and Inertial Impaction. The gravitational settling velocity, the particle terminal velocity, is expressed by the product υg ) gτp, where g is the gravitational accelaration, τp ) 2Fpr2pCc/9µ is the particle relaxation time (37), and Fp is the particle density. Inertial impaction of particles on the airway walls in bifurcation zones is accounted for by the expression developed by Cheng and Wang (42) for particle collection efficiency in smooth bends. Here, the bends extend to a length determined in terms of the input branching angle. Interception deposition efficiency, which is used to account for the finiteness of real particles as opposed to referring to point particles in the Eulerian approach, is neglected in the simulations because of the smallness of the governing parameter 2rp/da , 1. Interception is considered to be an efficient deposition mechanism only for elongated and fibrous particles (43). Vapor Deposition. The vapor deposition flux can be defined as a product of the deposition velocity and the difference in vapor concentration between the airway lumen and walls (41). Taking the pressure gradient that drive the condensation flux to be the pressure difference across the concentration boundary layer, the deposition flux of species j, Fv,j, can be expressed as (27)
(19)
where Dj is the molecular diffusion coefficient of species j and δc the thickness of the concentration boundary layer. Like the thermal boundary layer, the concentration boundary layer is given in terms of a dimensionless transport parameter
δc ) da/Sh
(20)
where Sh, the Sherwood number, is the equivalent of the Nusselt number in mass transport processes. Since mass and heat transport are comparable in air, entrance effects are incorporated in the calculation of the concentration boundary layer thickness by defining a position-dependent Sherwood number whose asymptotic value is 3.66. Reactions on Surfaces of the Human Airways. Calculations of chemical reaction rates account for two processes: (i) diffusion of vapor from the airway lumen to the wall and (ii) chemical reaction between different species. The “overall” (effective) chemical reaction rate can be expressed as
Rchem )
(
1 1 + vv,j kchem
)
-1
pvj Mw,j RT
(21)
where kchem is the net reaction rate on the airway surfaces. Kinetically controlled reactions of different gaseous species are included in a submodule of the model, for predicting the fate of highly reactive inhaled agents. A detailed description of the model numerics can be found in the Supporting Information.
The Architecture of the Human Lungs The architecture of the human lungs incorporated in the model is based on Weibel’s Model A (45). This morphological model of the tracheobronchial airways forms the basis of
most lung ventilation and inhalation dosimetry models (14). The tracheobronchial tree and the pulmonary region are described as a regular symmetrically dichotomous branching system of straight, smooth, circular tubes. Every bifurcation results in two identical daughter airways, usually of smaller diameter than the parent airway. The tree-structure of the lung model consists of 24 airway generations, the trachea being the generation i ) 0 and the alveoli and alveolar sacs, the fundamental gas exchange units, being lumped in the last generation i ) 23. As airways successively bifurcate, there are 2i airways in each generation for a total of 224 - 1 airways. Since the total cross-sectional area usually increases with penetration distance, the average air velocity in distal airways decreases. Airway gravity angles, which in single pathway lung models are uniform within any generation, vary in accordance with the specific body posture considered (upright position, sleeping position, etc.) In particular, the generation-specific gravity angles implemented in the model for a standing (upright) position follow Yeh and Schum’s (16) typical pathway lung model. Deterministic symmetrical branching causes all airways in any given generation to have identical dimensions and location along the thoracic pathway (a single pathway lung model). Clearly, this picture is based on a simplification of the actual morphology of the human tracheobronchial tree in which irregular branchings are usually encountered, and where respiratory disease may cause changes to the human airway geometry, such as inducing variation of the thickness of the mucus layer lining the ciliated airways (14). To account for the large intersubject morphological variability while maintaining a simple geometrical anatomic representation of the human lungs, we implemented within the Weibel’s lung model bifurcation zones that are characterized by branching angles from Yeh and Schum’s (16) typical pathway lung model. The lack of comprehensive morphological data on the human lung, and in particular of morphometric studies on lungs of children, the elderly, and individuals with chronic obstructive pulmonary disease suggests that present model predictions are more suitable for healthy adults rather than for individual at the high risk end. Inded, Weibel’s lung model is generally accepted as a representative model of the mean morphology of a large normal population, rather than as a description of the lung of a specific individual. Hence, model predictions should be considered representative (average) of a large population. Clearly, by using morphometric data of a specific individual, obtained by MRI, CT, or other medical imaging techniques, model results can be tuned to predict particle deposition in the individual’s lung. Future plans include replacing the single pathway lung module with a stochastic lung model for which anatomical and physiological inputs are obtained from pertinent distributions (46). This will allow the calculation of particle deposition as a distribution characterized by a mean, a standard deviation, and a skewness for any given set of physicochemical properties of the inhaled aerosol. Although the present model is general enough to permit the use of time-dependent air flows, to simplify the computation unidirectional air flow with flat velocity profile (plug flow) is specified in the airways. To make the results more realistic, however, boundary layer effects are incorporated by introducing appropriate dimensionless transport parameters for estimating the corresponding thermal and species boundary layers. Different breathing rates are considered, corresponding to sedentary and extensive exercise conditions.
Results Comparison of Model Predictions with Experimental Data. For evaluation of the purpose of model, theoretical predictions were compared to independently published experimental data. The experimental data were obtained by state-
FIGURE 1. Comparison of tracheobronchial deposition data from Emmett et al. (49) with theoretical predictions. Model parameters include branching angles of 45°, tidal volume of 1000 cm3, inert aerosol initially composed of particles having a log-normal size distribution with a 1 µm count median diameter, and a 1.7 geometric standard deviation. of-the-art techniques and are considered very reliable data. Over the years these data were used for evaluation of numerous dosimetry models (47, 48). In fact, empirical correlations of particle deposition extracted form these data are used in the International Commission on Radiation Protection ICRP66 recommended dosimetry model (15), the predictions of which are also compared with results of the present fully mechanistic model. Evaluation of the model against more recent experimental data is under progress. Emmett et al. (49) measured total and regional deposition of monodisperse aerosols in the respiratory tract of 12 healthy subjects that breathed through the mouth in an upright position. Radioactively labeled polystyrene particles with aerodynamic diameters ranging between 3.5 and 10.0 µm and density of 1060 kg/m3 were used. The total deposition rate increased with increasing particle size. In the simulations, downward air flow corresponding to a 1000 cm3 tidal volume (in agreement with the experimental conditions) and 45° branching angles were assumed. The aerosol considered comprised inert particles with 1 µm count median diameter (CMD) and 1.7 geometric standard deviation (GSD). Particle agglomeration and deposition were the only mechanisms considered. Model predictions are in qualitative and quantitative agreement with the experimental deposition profile in the tracheobronchial region (Figure 1), even though the experimental data are quite scattered. The increase in deposition with increasing particle size for the size range considered in Figure 1 is documented in the literature and is in agreement with predictions of the ICRP66 (15). Chan and Lippmann (50) used monodispersed Fe2O3 particles with mass median aerodynamic diameter larger than 2 µm and density of 2.56 g/cm3 for in vivo deposition determination in 26 healthy nonsmokers and for in vitro measurements in hollow casts. The tidal volume was approximately 1000 cm3. To compare model predictions with experimental data it was assumed that the tidal volume was inspired in 1 s, corresponding to an air flow rate of 1000 cm3/s. This is a relatively high air flow, corresponding to nonsedentary breathing conditions under which respired particles undergo enhanced inertial deposition. The deposition fraction in the 18th generation (the first alveolated airway in Weibel’s model) was chosen to surrogate deposition in the alveolar region. Since diffusional deposition calculations in the model tend to diverge for minuscule volumes and vanishing air velocities such as those found in the last airway VOL. 35, NO. 18, 2001 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
3731
FIGURE 2. Model predictions of tracheobronchial deposition compared to predictions of the ICRP66 empirical model (15) and to the experimental data of Lippmann (51) and Chan and Lippmann (50). Theoretical predictions were obtained for air flow of 1000 cm3/s and bifurcation angles of 45° (typical deposition conditions), 65° (maximum deposition conditions), and 25° (minimum deposition conditions).
FIGURE 4. Predicted thracheobronchial deposition vs experimental data. Model parameters include air flow of 1000 cm3 s-1, log-normally distributed inert aerosol with an initial 1 µm count median diameter, and a 1.7 geometric standard deviation, and bifurcation angles of 45° (typical deposition conditions), 65° (maximum deposition conditions), and 25° (minimum deposition conditions).
FIGURE 3. Model predictions of alveolar deposition compared to predictions of the ICRP66 empirical model (15) and to experimental data of Lippmann (51) and Chan and Lippmann (50). Theoretical predictions were obtained for air flow of 1000 cm3/s and bifurcation angles are 45° (typical deposition conditions), 65° (maximum deposition conditions), and 25° (minimum deposition conditions).
FIGURE 5. Predicted alveolar deposition vs experimental data. Model parameters include: air flow of 1000 cm3 s-1, log-normally distributed inert aerosol with an initial 1 µm count median diameter and a 1.7 geometric standard deviation, and bifurcation angles of 45° (typical deposition conditions), 65° (maximum deposition conditions), and 25° (minimum deposition conditions).
generations, the accuracy of model predictions decreases in this region. Model predictions in Figures 2 and 3 were obtained for three typical bifurcation angles θ, the only parameter not rigidly fixed in Weibel’s lung model. Typical deposition conditions correspond to θ ) 45°, maximum (most favorable) deposition conditions correspond to θ ) 65°, and minimum (least favorable) deposition conditions to θ ) 25°. Figure 2 shows that particle deposition data in the tracheobronchial region seems bound by the theoretical results for the most and least favorable deposition conditions. At present, varying the bifurcation angles affects only inertial deposition, which takes place mainly in the proximal airways. Therefore, alveolar deposition is less affected by these morphological variations (Figure 3). Overall, the natural variability of the lung morphology, which is transformed into model uncertainty with respect to the bifurcation angles, significantly affects the particle deposition profile along the tracheobronchial airways. Indeed, deposition of large particles in the upper airways can double under adverse conditions (i.e. larger surfaces for impaction in the carinal ridges). In contrast,
deposition of small particles, which are less affected by inertial impaction, is less affected by the variation in branching angels. Figures 4 and 5 show model predictions of deposition in the tracheobronchial and the alveolar regions, respectively, vs data from different studies (49-55). The experimental data were obtained by tracking radioactive labeled poorly soluble particles of diameter larger than 0.1 µm. In vivo measurements were based on the amount of radioactivity retained in the lung as a function of time, where the “fast-cleared” and “slow-cleared” deposition fractions correspond to deposition in the ciliated tracheobronchial and alveolar regions, respectively (14). Although experimental data are quite scattered, owing both to nonconsistency in methodology (56, 57) and to intersubject variability (14), the theoretical predictions are in fair agreement with the general deposition pattern in the tracheobronchial region (Figure 4). Since these data are considered to be the benchmark for comparisons between predicted and measured lung deposition of thoracic and respirable particles, the model seems to be adequate in describing the upper and lower deposition bounds. With
3732
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 35, NO. 18, 2001
new data recently emerging on deposition of ultrafine particles in the human respiratory tract, the model prediction power for particles in the size range 0.01-0.5 µm will be investigated. Evaluation of Model Predictions Against Results of the ICRP66. The present mechanistic model was evaluated against the International Commission on Radiological Protection ICRP66 benchmark model (15), the recommended model for predicting deposition of particulate matter in the human respiratory tract. The ICRP66 model is an empirical model, which is based on experimental deposition data of a limited particle-size range (1-10 µm aerodynamic diameter) in adult Caucasian males, combined with semiempirical expressions for calculating regional deposition in the human airways. The respiratory tract is modeled by two extrathoracic and three thoracic compartments. Deposition results from gravitational settling, inertial impaction, and Brownian diffusion. Figures 2 and 3 show predictions of both the new and the ICRP66 models against Lippmann’s (51) and Chan and Lippmann’s (50) experimental data on tracheobronchial and alveolar deposition. Agreement with the experimental data is satisfactory for both models. The close agreement of the present mechanistic model with the experimental data (Spearman’s rank correlation rSP ) 0.31), despite the numerous simplifications used and the large number of physiological inputs, may suggest that the underlying nature of the physicochemical processes responsible for PM deposition in the human airways is captured by the model. In contrast, agreement of the ICRP66 model predictions with the experimental results (Spearman’s rank correlation rSP ) 0.6) is expected, since the model was tuned to reproduce these data. The nature of empirical models is that they are based on, and therefore fit, experimental data. As such, it is rarely possible to extrapolate empirical models to conditions different from those for which they were developed. In contrast, mechanistic models provide the possibility for spatially detailed calculations, parametric study of the system, and extrapolation to different conditions and species. Obviously, there are many more requirements for input parameters of the lung’s structure and physiology in mechanistic models. As discussed above, the limitation of the experimental data with respect to both its uniform and passive chemistry, truncated polydispersity, and huge variability among subjects and among studies does not provide adequate opportunity to demonstrate the superiority of mechanistic vs empirical models for deterministic (single prediction per particle size) modeling. Sensitivity Analysis. The variability in lung anatomy, breathing physiology, environmental conditions, and aerosol properties are usually transformed into model uncertainty in terms of incomplete knowledge of model parameters and inputs (58). A complete study of the propagation of data uncertainty through the model is beyond the scope of this work, yet the promise of such an analysis is demonstrated here by a simple sensitivity analysis of selected parameters involved in the calculation. The carinal ridges appear to be “hot spots” for particle deposition. Figure 6 depicts the size distribution of persistent nongrowing particles at the end of the 16th airway generation for two uniform bifurcation angles along the lung pathway. In agreement with previous work (14), higher particle deposition occurs for wider bifurcation angles. Thermophoretic deposition is studied assuming a constant airway wall temperature of 36 °C and a hypothetical temperature profile of the inhaled air, linearly decreasing from 55 °C. For such a high temperature difference thermophoresis has a pronounced effect on the size distribution evolution, enhancing deposition along the airways. When thermophoresis is considered, the mass distribution of the
FIGURE 6. Particle mass distribution at the end of the 16th airway generation (the ultimate conducting airways in Weibel’s model A for two uniform bifurcation angles. Particle initial size is log-normally distributed with 1 µm count median diameter and a 1.7 geometric standard deviation.
FIGURE 7. Effect of thermophoresis on the particle mass distribution at the end of the 16th airway generation (the ultimate conducting airways in Weibel’s model A). The airway wall temperature is 36 °C, whereas the hypothetical temperature profile of the inhaled air decreases linearly from 55 °C such that the temperature difference in the 16th airway generation is 3 °C. persistent aerosol is smaller than that calculated under isothermal conditions (Figure 7). Clearly, smaller temperature differences will have a lesser effect on both the deposition profile and the persistent aerosol size distribution. Accounting for the thermophoretic effect under realistic conditions is under progress.
Discussion Differences in airway morphology and breathing pattern among subjects significantly affect the deposition of inhaled particles. Because of the complex nature of the lungs, mathematical modeling of aerosol transport and fate in the human respiratory tract is based on many idealized assumptions, which constitute a set of approximations to the actual system. Here, one-dimensional airflow is assumed, with boundary layer effects incorporated in the calculation by introduction of appropriate dimensionless transport parameters. Experimental deposition data show considerable scatter, which results from large intersubject and intertechnique variabilities. To address this issue model results can provide upper and lower estimates of particle deposition. Specifically, VOL. 35, NO. 18, 2001 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
3733
the envelope bounding deposition due to a single morphological change was demonstrated. This approach can be extended to study the n-dimensional bounding surface due to independent or synergistic changes in model parameters. Unlike clinical methods, mathematical models provide an efficient option for sensitivity analysis of lung deposition, where processes and mechanisms that govern the inhalation route of the exposure-dose-response scenario can be examined independently. Accordingly, the model accounts for and keeps track of the contribution of different processes to the overall deposition profile. In particular it directly considers the wide polydispersity, multimodality, and heterogeneous composition of common aerosols. As such, the model may be useful for regulatory policy development in the fields of industrial hygiene and environmental and occupational health, in predicting the dose resulting from exposure via inhalation to airborne pollutants such as fugitive dust, fine diesel PM, agrotoxic sprays, lubricating fluid mists, and bioactive airborne agents and in analyzing clinical data and protocols in the field of aerosol therapy.
(19) (20) (21) (22) (23)
Acknowledgments
(33)
This work was supported by the U.S. DOE Cooperative Agreement DE-FC01-95EW55084 (CRESP), the U.S. EPAfunded Center for Exposure and Risk Modeling (CERM), and U.S. EPA Cooperative Agreement CR827033.
(34)
Supporting Information Available A detailed description of the model numerics. This material is available free of charge via the Internet at http:// pubs.acs.org.
Literature Cited (1) Georgopoulos, P. G.; Walia, A.; Roy, A.; Lioy, P. J. Environ. Sci. Technol. 1997, 31, 17-27. (2) Roy, A.; Weisel, C. P.; Lioy, P. J.; Georgopoulos, P. G. Risk Analysis 1996, 16, 147-160. (3) Roy, A.; Georgopoulos, P. G. Technical Report CCL/EDMAS-02; Environmental and Occupational Health Sciences Institute: Piscataway, NJ, 1997. (4) Isukapalli, S. S.; Georgopoulos, P. G. Submitted for publication, 2000. (5) Pope, C. A.; Schwartz, J.; Ransom, M. R. Arc. Environ. Health 1992, 47, 211-217. (6) Schwartz, J.; Dockery, D. W. Am. Rev. Respir. Dis. 1992, 145, 600-604. (7) Schwartz, J. Fundam. Appl. Toxicol. 1994, 64, 26-35. (8) Pope, C. A.; Thun, M. J.; Namboodiri, M. M.; Dockery, D. W.; Evans, J. S.; Speizer, F. E.; Heath, D. W. Am. J. Respir. Cr. C. Med. 1995, 151, 669-674. (9) Schlesinger, R. B. Inhal. Toxicol. 1995, 7, 99-109. (10) NRDC. Natural Resources Defence Council Report; Washington, DC, 1996. (11) McClellan, R. O.; Miller, F. J. Chem. Ind. Institute Toxicol. 1997, 17, 1-23. (12) Schwartz, J.; Dockery, D. W.; Neas, L. M. AWMA 1996, 46, 2-12. (13) Wilson, W. E.; Suh, H. H. AWMA 1997, 47, 1238-1249. (14) U.S. Environmental Protection Agency. Air Quality Criteria for Particulate Matter; Technical Report EPA/600/P-95/001aF; Research Triangle Park: NC, 1996. (15) International Commission on Radiological Protection. Human Respiratory Tract Model for Radiological Protection; Pergamon Press: Oxford, U.K., 1994. (16) Yeh, H. C.; Schum, G. M. Bull. Math. Biol. 1980, 42, 461-480. (17) Koblinger, L.; Hofmann, W. J. Aerosol Sci. 1990, 21, 661-674. (18) Martonen, T. B.; Zhang, Z. Inhal. Toxicol. 1993, 5, 165-187.
3734
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 35, NO. 18, 2001
(24) (25) (26) (27) (28) (29) (30) (31) (32)
(35) (36) (37) (38) (39) (40) (41) (42) (43) (44) (45) (46) (47) (48) (49) (50) (51) (52) (53) (54) (55) (56) (57) (58)
Taulbee, D. B.; Yu, C. P. J. Appl. Physiol. 1975, 38, 77-85. Yu, C. P. J. Aerosol Sci. 1983, 5, 599-609. Nixon, W.; Egan, M. J. J. Aerosol Sci. 1987, 18, 563-579. Darquenne, C.; Paiva, M. J. Appl. Physiol. 1994, 77, 2889-2898. Rudolf, G.; Gebhart, J.; Heyder, J.; Schiller, C. F.; Stahlhofen, W. J. Aerosol Sci. 1986, 17, 350-355. Broday, D. M.; Georgopoulos, P. G. Aerosol Sci. Technol. 2001, 34, 144-159. Batchelor, G. K.; Shen, C. J. Col. Int. Sci. 1985, 107, 21-37. Williams, M. M. R.; Loyalka, S. K. Aerosol Science, Theory and Practice; Pergamon Press: Oxford, U.K., 1991. Hontanon, E.; Lazaridis, M.; Drossinos, Y. J. Aerosol Sci. 1996, 27, 19-39. Asgharian, B.; Anjilvel, S. Aerosol Sci. Technol. 1995, 22, 261270. Pedley, T. J.; Kamm, R. D. In The Lung: Scientific Foundations; Crystal, R. G., West, J. B., Eds.; Raven Press: NY, 1991; Vol. 1, pp 995-1010. Lazaridis, M.; Drossinos, Y. J. Phys. A: Mathematical General 1997, 30, 3847-3865. Lazaridis, M.; Kulmala, M.; Gorbunov, B. J. Aerosol Sci. 1992, 23, 457-466. Fuchs, N. A. Evaporation and Droplet Growth in Gaseous Media; Pergamon Press: NY, 1959. Jokiniemi, J.; Lazaridis, M.; Lehtinen, K. E. J.; Kauppinen, E. J. Aerosol Sci. 1994, 25, 429-446. Mason, B. J. The Physics of Clouds; Clarendon Press: Oxford, U.K., 1971. Lazaridis, M.; Koutrakis, P. J. Aerosol Sci. 1997, 28, 107-119. Lazaridis, M.; Ford, I. J. J. Chem. Phys. 1993, 99, 5426-5429. Friedlander, S. K. Smoke Dust and Haze; Wiley & Sons: NY, 1977. Wagner, P. E.; Kerker, M. J. Chem. Phys. 1977, 66, 638-646. Saffman, P. G.; Turner, J. S. J. Fluid Mech. 1956, 1, 16-30. Im, K. H.; Chung, P. M. AIChE J. 1983, 29, 498-505. Hinds, W. C. Aerosol Technology; Wiley & Sons: NY, 1982. Cheng, Y. S.; Wang, C. J. Atmos. Environ. 1981, 15, 301-306. Broday, D.; Fichman, M.; Shapiro, M.; Gutfinger, C. Phys. Fluids 1998, 10, 86-100. Im, K. H.; Ahluwalia, R. K.; Chung, P. M. Aerosol Sci. Technol. 1985, 4, 125-140. Weibel, E. R. Morphometry of the Human Lung; Academic Press: NY, 1963. Koblinger, L.; Hofmann, W. Phys. Med. Biol. 1985, 30, 541-556. Hofmann, W.; Koblinger, L. J. Aerosol Sci. 1990, 21, 675-688. Edwards, D. A. J. Aerosol Sci. 1995, 26, 293-317. Emmett, P. C.; Aitken, R. J.; Hannan, W. J. J. Aerosol Sci. 1982, 13, 549-560. Chan, T. L.; Lippmann, M. Am. Ind. Hyg. Assoc. J. 1980, 41, 399-409. Lippmann, M. In Handbook of Physiology; Lee D. H. K., Falk H. L., Murphy, S. D., Geiger, S. R., Eds.; American Physiological Society: Philadelphia, PA, 1977; pp 213-232. Foord, N.; Black, A.; Walsh, M. J. Aerosol Sci. 1978, 9, 343-357. Stahlhofen, W.; Gebhart, J.; Heyder, J. Am. Ind. Hyg. Assoc. J. 1980, 41, 385-398. Stahlhofen, W.; Gebhart, J.; Heyder, J. Am. Ind. Hyg. Assoc. J. 1981, 42, 348-352. Stahlhofen, W.; Gebhart, J.; Heyder, J.; Scheuch, G. J. Aerosol Sci. 1983, 14, 186-188. Stahlhofen, W.; Gebhart, J.; Heyder, J.; Scheuch, G. J. Aerosol Sci. 1986, 17, 333-336. Smaldone, G. C.; Perry, R. J.; Bennett, W. D.; Messina, M. S.; Zwang, J.; Ilowite, J. J. Aerosol Med. 1988, 1, 11-20. Isukapalli, S. S.; Roy, A.; Georgopoulos, P. G. Risk Analysis 1998, 18, 351-363.
Received for review August 2, 2000. Revised manuscript received May 30, 2001. Accepted May 30, 2001. ES001545W
