INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

CONVECTION AND DIFFUSION ARE the two major mechanisms of mass transport for gas molecules and submicrometer aerosols in the pulmonary acinus. For gas transport, diffusion dominates at distances comparable to acinar size and over times comparable to breathing frequencies, and, therefore, theories based on diffusion are probably adequate. By contrast, the particle diffusivity of submicrometer-sized aerosols is very small, and, therefore, acinar convection, even though it is in a quasi-Stokes viscous flow regime (26), is correspondingly more important and may dominate aerosol transport. However, current theories describe aerosol transport as a dispersion (diffusion-like) process (e.g., Refs. 7, 10, 11, 19, 34). These theories are based on the following two key assumptions: 1) acinar flow is basically kinematically reversible (i.e., during expiration each fluid particle retraces the path taken during inspiration) (9, 45), and 2) all processes (including the coupling of Brownian diffusivity with the convective flow field and any kinematic irreversibility that may be present) that contribute to irreversible aerosol bolus spreading can be characterized as axial mixing with an effective longitudinal diffusivity (D_eff). The first assumption is based on classical fluid mechanics (36), and the second assumption is substantially equivalent to Taylor dispersion (35). As most aerosol studies are currently interpreted in the framework of these dispersion theories, experimental data are often reduced and analyzed through the use of some D_eff (e.g., Ref. 30), and many of the recent theoretical research efforts are focused on refining D_eff for better fit to experimental data, through which new insights into acinar transport mechanisms are sought (e.g., Ref. 6).

Our laboratory's recent findings (37-41), as well as those of others (8, 18), however, have questioned the validity of the basic assumptions that formed the basis of the dispersive theories. We have demonstrated that, because of the peculiar geometry of the alveolated duct and its time-dependent motion associated with tidal breathing, under certain conditions, alveolar flow can be chaotic (16, 40). As a consequence, acinar flow can be kinematically irreversible, even though it is a low-Reynolds number viscous flow, and lung expansion and contraction are approximately self-similar and reversible (1, 13, 14, 24, 46). These findings are supported by a discovery in fluid mechanics that chaotic mixing can occur even in a viscous flow (2, 25). We have also developed a theoretical analysis (5) to show that chaotic mixing is radically different from diffusive mixing. In chaotic acinar flow, a tracer bolus undergoes cyclic stretch-and-fold deformation, resulting in the induction of finer and finer scales in tracer profile with repeated breaths. This process soon reaches a critical moment at which the lateral distance between adjacent tracer striations becomes comparable to the diffusion distance. A burst of mixing occurs at this moment, and mixing is quickly completed.

The objectives of the studies reported here are, through numerical simulations of bolus experiments, to provide key data showing that the behavior of particles in a rhythmically expanding, multiply alveolated duct flow, with a saddle point and associated vortexes in each air pocket, does not satisfy the fundamental assumptions of any dispersion theory and that the shape of a tracer bolus evolves to a stretch-and-fold fractallike pattern, similar to those found in flow visualization experiments in rat lungs (Tsuda A, Butler JP, and Rogers RA, unpublished observations). The results suggest that 1) kinematic irreversibility is the origin of aerosol transport, 2) axial transport cannot be characterizable by an effective diffusivity, and 3) fractal trajectories can occur in most of the alveoli in the acinar tree. The alternative mechanism of aerosol transport that we propose here may, in fact, be the dominant mechanism determining deposition of submicrometer particles deep in the lung.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

In a previous investigation (40), we used the single-alveolus model to explore the basic physics operating in a viscous flow subjected to cyclic alveolar wall motion. The alveolus model used in the present investigation is also axisymmetric but comprises a central circular channel around which are placed nine tori, equispaced in the axial direction (Fig. 1A). Details of a typical cell are given in Fig. 2. The duct and alveolar walls move in a perfectly kinematically reversible, simple sinusoidal manner with a specific volume excursion C of 25% [C = (V_max  V_min)/V_min, where V_max and V_min are the maximum and minimum volumes of the model, respectively] and a cycle period T of 3 s. These correspond roughly to typical tidal ventilation and respiratory period in human. Any length scale L of the model changes as L(t) = [1 + K sin(nt)], where is the mean L value, t is time, n = 2/T, and K = (  1)/( + 1), where . A smooth-walled, closed-end duct is fitted to the distal end of the model. The closed-end duct is used to approximate the airway distal of the site of interest up to the terminal alveoli, and it is also used to control the Reynolds number of the flow in the model. The walls of this duct also move in the same sinusoidal manner as the main alveolus model. The mean length of this duct, _D, controls the bulk velocity U of the fluid entering the model; i.e., U = /R = 3nK_D cos(nt), where is the volume flow rate and R_D is the duct radius. The root mean square (RMS) Reynolds number Re_RMS = U_{RMS D}/, where U_RMS is the RMS U and equals U_max/, where U_max is maximum U, _D is the mean R_D, and  is the kinematic viscosity. Re_RMS ranges from 0.006 to 0.728 in the alveolated region of the model. Using a closed-end duct has the added advantage of avoiding the difficulty of defining appropriate boundary conditions at the downstream boundary of the alveolated section. Typically, L_D/R_D is >100, and, hence, it can be assumed that the end of the closed-end duct was sufficiently far from the downstream boundary of the alveolated section so as not to affect the flow across this boundary adversely. As a further refinement to the original model, the sharp corner at the intersection between the alveolus and the duct is replaced by a more natural circular section. For comparative purposes, two other cases are also simulated. One is the flow in an isolated alveolus model (Fig. 1B), similar to that used in Tsuda et al. (40), and the other is the flow in a nonalveolated, rhythmically expanding straight tube (Fig. 1C).

View larger version (12K): [in this window] [in a new window]   Fig. 1.   Schematics of moving walled models. A: 9-cell alveolated duct model. B: isolated cell alveolated duct model. C: nonalveolated smooth-wall model.

View larger version (14K): [in this window] [in a new window]   Fig. 2.   Expanded view of a typical alveolar cell. RD (= 250 µm), duct radius; RA (= 200 µm), alveolar radius; LA (= 346 µm), alveolar opening length; LC (= 413 µm), alveolar cell length; D, ductal volume flow rate; A, alveolar volume flow rate;  (= 60°), alveolar half-opening angle; CL, duct center line. Note that the alveolar corners were rounded with the corner radius = 0.02 RD.

The flow field is defined by the full, incompressible Navier-Stokes equations, which are solved numerically on a multiblock, body-fitted moving grid using the finite volume code CFX-4 (CFDS, AEA Technology, Harwell, UK). This general-purpose, pressure-correction code offers a variety of discretization schemes and solution techniques. In these calculations, central differencing is used to model the convection terms, and the implicit backward Euler method is used to advance the solution in time. A combination of the SIMPLEC pressure-correction method (44) and the Rhie-Chow (28) algorithm to eliminate pressure oscillations on the collocated gird (12) is used in the formulation of the discrete equations. Stone's method (33) is used to solve the discrete velocity equations, and the method of preconditioned conjugate gradients (see, for example, Ref. 22) is used to solve the discrete pressure correction equation. At the inlet, a constant-pressure boundary is defined. The value of the pressure on this boundary is set arbitrarily to zero, and, as for incompressible flow, only the pressure gradient is of importance. The no-slip condition is enforced on all solid surfaces, which, in the case of moving walls, means that the fluid matches the wall velocity at the fluid-wall interface. Tests were carried out to ensure that the solutions are grid independent and converged. For example, increasing the number of grid cells by 30% above that which was eventually used produces an increase of only 0.16% in the predicted maximum velocity. The final 73-block grid had a total of 38,349 active cells. With the use of a measure of error due to using backward Euler time stepping, suggested by Roache (29), a time step of T/240 is found to give sufficiently accurate results in that further reductions in the size of the time step used does not produce any significant reduction in the error. These simulations are computationally intensive, with one breathing cycle taking ~350 h on a Sun Ultra 10.

Particle trajectories are calculated in all three models using a special-purpose tracking routine. This routine reads a full cycle of flow field and grid data produced by CFX-4 and uses this data repeatedly, cycle after cycle, and a predictor and corrector method to track individual particles (fluid elements) over as many cycles as required. The time step used in the particle tracking routine is set independently of that for the flow-field solution. Before the particle track is advanced to the next time, the time step is recalculated, using local flow conditions, to ensure that the particle does not step out of the solution domain. At each particle-track time, a grid and flow field are created from the CFX-4 data using bicubic interpolation in space (27) and linear interpolation in time. For numerical efficiency, the routine solves the particle tracks in a stationary computational space, and, at each time step, the particle position is mapped back to physical space before it is written to an output file. The maximum error in the predicted particle position is estimated to be 0.2% of the distance traveled by the particle, and, in most cases, the error is considerably smaller than this.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Flow patterns. Solving the velocity field of the carrier gas on a moving grid over the physiologically relevant range of flow parameters (Re_RMS < 1) in a rhythmically expanding and contracting, multiply alveolated duct, we often detected the presence of slowly rotating recirculation in each alveolus (see Fig. 3, A and B). The size of the recirculation flow depends on the ratio between the alveolar flow (_A) (i.e., flow produced by the volume change of the alveolus) and the volumetric ductal flow (_D) (_A/_D). Similar to our previous study (40), we found that the smaller the _A/_D, the larger the alveolar recirculation. When _A/_D is larger than ~0.1, however, the alveolar flow is largely radial without recirculation (Fig. 3C). Because, under normal breathing conditions or during moderate exercise, the value of _A/_D is usually <0.05 in the majority of alveoli along the acinar tree (from the respiratory bronchioles to the last few generations) (37), we expect that most of the alveoli are likely to possess recirculation in their flow field. Importantly, the presence of alveolar recirculation in the cyclically expanding alveoli is topologically associated with the existence of a stagnation saddle point in each alveolar flow field (40), implying that chaotic mixing could originate in most of the alveoli (see DISCUSSION).

View larger version (19K): [in this window] [in a new window]   Fig. 3.   Typical streamline patterns at the peak inspiration with various flow conditions. A: Reynolds number root mean square (ReRMS) = 0.714, ratio of A to D (A/D) = 0.005. B: ReRMS = 0.077, A/D = 0.047. C: ReRMS = 0.030, A/D = 0.119. Insets: saddle points (arrowhead) near the proximal corner of alveoli.

Interalveolar kinematic mixing. To demonstrate the effects of a series of saddle points on fluid flow irreversibility, the motion of massless particles was tracked over one ventilation cycle in a nine-cell alveolar model for three different ranges of _A/_D (0.0050 < _A/_D < 0.0053, 0.040 < _A/_D < 0.069, and 0.081 < _A/_D < 0.577 shown in Fig. 4, A, B, and C, respectively). The particles were initially placed on radial lines across the duct midway between alveoli at three different initial axial locations (shown as three different colors, brown, green, and pink, in Fig. 4). As soon as inspiration begins, particles near the center line convect distally, proportionally to the bulk mean velocities (see t/T = 0.25 in Fig. 4). The lines of particles quickly approach each other and comigrate along the central channel, particularly in the cases of smaller _A/_D (Fig. 4, A or B). The particles that were located initially near the alveolar opening enter the alveolus (Fig. 4). The depth of particle penetration into the alveolus during inspiration depends on the size of alveolar recirculation. When the alveolar recirculation is large (Fig. 4A), the particles penetrate deep into the alveoli, and the particles from the three different lines appear to be mixed during inspiration (see t/T = 0.25 and 0.5 in Fig. 4A).

View larger version (23K): [in this window] [in a new window]   Fig. 4.   Interalveolar mixing. Massless (fluid) particles were tracked, starting from 3 different initial axial locations (shown as 3 different colors) with 3 different flow conditions. A: 0.0050 < A/D < 0.0053, 0.728 > ReRMS > 0.690. B: 0.040 < A/D < 0.069, 0.092 > ReRMS > 0.053. C: 0.081 < A/D < 0.577, 0.045 > ReRMS > 0.006. Five different time points in each simulation were shown: time (t)/cycle period (T) = 0 (initial), t/T = 0.25 (peak inspiration), t/T = 0.5 (end inspiration), t/T = 0.75 (peak expiration), and t/T = 1.0 (end expiration). Note that, in A and B at end inspiration, many particles are distal of the alveolated portion of the model and hence are not visible in the views given.

Motion of interface between inhaled tracer fluid and the host alveolar residual fluid. Regarding massless tracer particles as marked fluid elements, the interface between one fluid and another can be approximated by chords connecting initially adjacent particles (Fig. 5A). In our studies, this piecewise linear approximation to the interface can represent the front of incoming tidal gas or a tracer bolus facing the host alveolar residual gas. As we have demonstrated in Fig. 4, A and B, the kinematically irreversible deformation of the tracer interface seems to be due to its association with alveolar recirculation flow over the respiratory cycle. Thus characteristics (e.g., size and strength) of the alveolar recirculation may be major determinants in these processes. Our studies described here focus on the situation that presumably occurs deep in the acinus, where the alveolar flow exhibits a medium- to small-sized recirculation (Fig. 3B). Larger alveolar recirculation flow (Fig. 3A) is likely to occur near the entrance of acinus (e.g., respiratory bronchioles), and that case is discussed elsewhere (38).

View larger version (27K): [in this window] [in a new window]   Fig. 5.   Fluid interface between incoming air and alveolar residual gas. A: initial position of fluid interface shown by straight-line chords connecting initially adjacent particles. For definition of symbols Q, R1, R2, R3, s1, s2, s3, Alv-1, Alv-2, Alv-3, see DISCUSSION. B: the interface shows "fingerlike" protrusions at the end of the first cycle. Inset 1: each protrusion consists of complex stretched-and-folded patterns. Inset 2: similar and finer stretched-and-folded patterns appear as the scale becomes finer. Flow conditions in the 9-cell model: 0.040 < A/D < 0.069, 0.092 > ReRMS > 0.053.

View larger version (15K): [in this window] [in a new window]   Fig. 6.   Axial variance (2) of particle distribution plotted vs. cycle number (N) as a family in the number of the particles (P). , P = 512; , P = 1,024; , P = 2,048; , P = 4,098; , P = 8,192; dashed and dotted lines show linear growth (2 = 0.0271N) and exponential growth [2 = 0.0474 (e0.457N  1)], respectively. Flow conditions in the 9-cell model: 0.040 < A/D < 0.069, 0.092 > ReRMS > 0.053.

View larger version (14K): [in this window] [in a new window]   Fig. 7.   Tracer stretching (L/Lo, where L is length and Lo is initial length) at end expiration plotted vs. N as a family in the P. , P = 512; , P = 1,024; , P = 2,048; , P = 4,098; , P = 8,192. Flow conditions in the 9-cell model: 0.040 < A/D < 0.069, 0.092 > ReRMS > 0.053.

View larger version (33K): [in this window] [in a new window]   Fig. 8.   A: distributions of massless (fluid) particles at N = 1 (green), N = 2 (blue), and N = 3 (red). Approximately 16,000 particles are shown for each N. B: box-counting analysis. A linear relationship on log[coefficient of variation (CV)] vs. log[edge length (E2)] with a slope of 0.2 shows that the pattern of particle distribution is fractal with a fractal dimension D 1.2. Flow conditions in the 9-cell model: 0.040 < A/D < 0.069, 0.092 > ReRMS > 0.053.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

The principal findings of this study are that 1) chaotic fluid motion occurring in a rhythmically expanding and contracting, multiply alveolated duct induces substantial kinematic irreversibility in the acinus, even under low-Reynolds number flow conditions, and 2) mixing due to this kinematic irreversibility is fundamentally different from the mixing described in dispersive processes.¹ A tracer (bolus) subjected to the chaotic flow field is deformed both axially and laterally. The ² of particle distribution increases exponentially, rather than linearly, with increasing cycle time; thus axial bolus spreading does not obey the basic rules described in classical diffusive transport theories (35). The cycle-by-cycle evolution of lateral particle distribution is even more complex. The tracer forms fingerlike protrusions, even after one cycle. The tracer length exponentially increases as N increased, forming characteristic stretch-and-fold fractal-like patterns.

Chaotic mixing in the pulmonary acinus. Viscous flow has been considered kinematically reversible if the boundary motion is reversible (36, 45). In the mid-1980s, however, there was a breakthrough discovery in fluid mechanics, namely, that even Stokes flow can be kinematically irreversible if the structure of the flow is chaotic (2, 25). In the last several years, applying this new concept to respiratory fluid mechanics, our laboratory has been studying the role of chaotic flow phenomena in the experimentally observed, yet theoretically unexplained, convective mixing occurring in the lung periphery (5, 16, 37, 40). In our laboratory's previous numerical study (16, 40), we reported that chaotic flow and chaotic mixing can occur in the alveolated duct because of its peculiar geometry and time-dependent motion associated with tidal breathing. We found that acinar flow was often slowly rotating in the alveolar air pocket, and the velocity field near the alveolar opening was complex with a stagnation saddle point typical of chaotic flow structure. Performing Lagrangian fluid particle tracking, we further demonstrated that, in such a flow structure, the motion of fluid, x(t), could be highly complex, irreversible, and unpredictable even though it was governed by simple deterministic equations [x(t) =  v(x,t)dt, x(0) = x_o, where v(x,t) denotes Eulerian velocity field].

Fingerlike protrusion. In this study, a line of massless particles was introduced in the alveolated duct to represent a fluid-fluid interface (Figs. 4 and 5A). The behavior of this interface, such as its reversibility and irreversibility, changes in its shape and size and contains crucial information for understanding the mechanism of mixing between inhaled particles and alveolar residual gas. By tracking the motion of the tracer particles, we have followed the motion of this interface over several cycles. Because the boundary of this line, shown as the point Q in Fig. 5A, is stationary on the wall because of the no-slip condition, the line expands when the alveolar walls expand during inspiration, and the line also tends to contract when the walls contract during expiration. However, the reversibility of this process depended on the nature of the flow fields sampled during expansion and contraction. The segments of the line near the channel center line are enormously stretched axially (see Fig. 4) and sample mostly reversible Poiseuille-like ductal flow fields (Fig. 3). Consequently, the interfacial line near the center line also shows approximately reversible behavior (Fig. 5B). By contrast, the segments of the lines near the walls (e.g., segments s1, s2, and s3 in Fig. 5A) sample a series of irreversible alveolar flow fields (e.g., Alv-1, Alv-2, and Alv-3, respectively, in Fig. 5A) and, consequently, do not return to their original positions (Fig. 5B). Each of these line segments, separated by points that sample approximately reversible flow fields between alveoli (e.g., R1, R2, and R3) basically form one finger after a cycle (Fig. 5B). The number of "fingers," therefore, roughly matches the number of alveoli distal to the initial position. Although, in the present computational model, the number of alveoli that the tracer samples is limited to six (due to computational constraints), in a real acinus, the inhaled bolus is expected to encounter a larger number of alveoli (~200) along the acinar longitudinal pathway. This implies that the acinar airways are likely to be filled with many longitudinal fingers after a few breathing cycles. This prediction has been confirmed experimentally in our laboratory's recent flow visualization studies performed in rat acini (39).

Axial phenomena-a departure from the conventional dispersion theory. In current theories, aerosol transport in the pulmonary acinus is described as a dispersion (diffusion-like) process, and the mixing phenomena are couched in the language of a D_eff. The most important feature of this approach is that, in any process that is dispersive in the sense that it can be characterized by an effective diffusivity, the variance of a bolus asymptotically increases linearly in time (or N). Phrased differently, the variance is an additive function over time. The reduction of experimental data through the use of some D_eff has thus been the framework by which many aerosol studies have been conducted (e.g., Ref. 30). For instance, bolus spreading in the tracheobronchial tree is commonly characterized by the difference between the inhaled and exhaled variances (proportional to H², where H is the bolus width at half height), given by , where H_exp is expiratory H and H_insp is inspiratory H (43). This numerical maneuver is clearly based on the underlying assumption of additivity of variances.

Lateral phenomena (fractal patterns). The results of our study show that complex flow phenomena can occur in the lateral direction. Because of alveolar flow irreversibility, the tracer (i.e., a series of line segments, s1, s2, etc., in Fig. 5A) deforms at every cycle, forming fine characteristic stretch-and-fold striation patterns (see Fig. 5B, insets), on top of the global fingerlike protrusion (Fig. 5B). Consequently, an enormously large lateral diffusion surface (which is represented as a stretched tracer line in our study) evolves over every cycle, suggesting a substantial enhancement of lateral particle transport and subsequent deposition on the acinar walls. Interestingly, we found that exact estimation of tracer length L was not possible because L was strongly dependent on the number of particles forming the tracer (Fig. 7). At sufficiently fine scales (i.e., when the number of test particles used is sufficiently large), the apparent L of the tracer increases exponentially with increasing N.

Physiological origins of mixing. Our laboratory has recently proposed two possible origins of "stretch-and-fold" kinematics in the pulmonary acinus: first, that induced by a small departure (asynchrony) from kinematically reversible motion of alveolar walls (16, 41), and, second, that due to the presence of saddle points associated with alveolar recirculation flows in the acinar flow field (40). With respect to the first mechanism, Miki et al. (24) reported the presence of a small but consistent geometric hysteresis (i.e., temporal asynchrony) in lung expansion during normal tidal ventilation in live rabbits. By matching this degree of geometric hysteresis, we have generated physiologically realistic asynchrony in physical models (41) and computational models (16), which demonstrated that geometrical hysteresis, even if small, can produce stretch-and-fold patterns and, consequently, induce substantial acinar flow irreversibility. In a related work, Smaldone et al. (31), using gravitational sedimentation of aerosol particles to estimate mean linear intercepts, showed significant geometric hysteresis in excised lungs with large-volume excursions from minimal volume and interpreted their data in terms of respiratory unit recruitment and derecruitment. Those studies, however, being restricted to single-volume histories, did not address the issue of mixing during cyclic ventilation and so are not strictly comparable to the present work. Finally, flow/volume hysteresis (different flow magnitudes at isovolume points on inspiration and expiration) can lead to differences in the velocity profiles in the central airways; this, in turn, can cause a difference in aerosol deposition between inspiration and expiration (4).

Significance and conclusions. It is well established that exposure to environmental aerosol pollutants is associated with health risks, ranging from mild to life threatening (47). The exposure and pathophysiological consequences are clearly linked by two independent and separate causal pathways: the exposure-dose relationship (i.e., given an aerosol concentration in the ambient air, what is the actual dose delivered to the lung) and the dose-response relationship (i.e., for a given deposited dose or burden, what is the biological consequence). Despite the large literature on exposure assessment methodologies as well as on the pathophysiological consequences of short- and long-term exposures, there is much less known about the mechanisms contributing to the first of these links, and much deposition and mixing data are inconsistent with previous theories of mixing deep in the lung. It is important, therefore, to identify new potential mechanisms that may dominate aerosol transport, even when boundary motion is approximately reversible. We argue in this paper that chaotic mixing is such a candidate. We have shown, through realistic numerical simulation of the low-Reynolds number alveolated duct flow (and by comparison with experimental results in rat lungs; unpublished observations), that the peculiar geometry of the alveolated duct structure within the pulmonary acinus and its cyclic motion during breathing can give rise to 1) a chaotic type of mixing associated with the presence of saddle points, 2) slow recirculatory flow within the alveoli, and 3) stretching and folding of stream surfaces. These, in turn, can significantly increase mixing, especially laterally, and will also contribute to an increasing ² (which increases faster than linearly with breath number, an observation that is inconsistent with any dispersal mechanism that can be characterized by an effective axial diffusivity). We suggest that chaotic mixing may be the dominant mechanism of aerosol transport and deposition deep in the lung.