| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1U2R2M, Unité de Recherche en Résonance Magnétique Médicale, CNRS, Univ Paris-Sud, Le Kremlin-Bicêtre; 2Centre de Recherche Claude Delorme (Air Liquide Research Center), Les Loges-en-Josas; 3Laboratoire Jacques Louis Lions, CNRS UMR 7598, Université Paris VI, Paris; 4Biomécanique Cellulaire et Respiratoire, Institut National de la Santé et de la Recherche Médicale UMR651, Université Paris XII, Créteil, France
Submitted 22 December 2005 ; accepted in final form 4 February 2007
 |
ABSTRACT |
|---|
 TOP ABSTRACT MATERIALS AND METHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
tracheobronchial tree; patient-based geometry; airway velocity profiles
In the early seventies, Olson et al. (35) measured axial velocities with a hot wire anemometer on two-dimensional (2D) planar cross sections located from the trachea down to the subsegmental bronchi in a post mortem airway cast extracted from a normal subject. This study demonstrated the importance of convective patterns of flow during inspiration. Ten years later, Chang and El Masry (3) and Isabey and Chang (21) explored, with the same measurement technique, not only axial velocities but also secondary motions in an idealized multiple-generation airway model based on Horsfield et al. (18) morphological study. They observed and quantified the complexity of flow, including the development of eddies influencing gas mixing. The swirls highly depend on the specificity of geometry such as the asymmetrical and three-dimensional (3D) character of the airway branching system. Such patterns have never been validated in vivo even within large airways and only preliminary results were recently obtained in the trachea of a volunteer during inhalation (9).
Lately, computational fluid dynamics (CFD) and medical imaging were combined to study airflow in more realistic morphologies. In a first step, medical imaging techniques, such as computed tomography (CT) lung scans, provided a wide range of possible morphologies and CFD studies were performed in realistic airway models (34, 39). In a second step, small airways were added to the image-based proximal airway model (26, 42). In these studies, the absence of a common reference anatomy precluded inter-study comparisons.
Since the emergence of rapid prototyping, realistic airway casts have been generated from medical imaging data. Clinkenbeard et al. (5) constructed a resin airway cast composed of five to eight generations with 78 end branches. However, most of the branches after the fourth generation did not have any lumen. This airway cast was used to calculate flow rates in airway branches using Pitot tubes to measure dynamic pressures and pressure drops between the tracheal inlet and the opened end branches (40).
Only a few experimental velocimetry techniques can be used to probe fluid dynamics on such complex geometries. Optical velocity mapping techniques like particle image velocimetry and laser Doppler velocimetry need a transparent model and a limited surface irregularity to perform accurate measurements (10). Hot wire anemometry introduces an invasive probe that is susceptible to perturbate the flow and that also needs a precise placement within the flow. Magnetic resonance velocimetry (30) has proven to be an efficient technique to map velocity in liquids containing proton 1H, and it was used both in vitro in various geometries and in vivo mainly for cardiovascular applications (17). Although, barely any MR signal can be measured from within the airways mainly because of the very low proton density. Hyperpolarized (HP) gas imaging (22), introduced in the mid nineties (1, 12), circumvents the intrinsic limitations of proton-based MR imaging in airways by using an HP gas that acts as a high MR signal source. This technique relies on drastically increased polarization (up to 105 times) to overcome the three orders of magnitude loss in atom density compared with water. Mapping gas velocity within the airways was only recently shown to be feasible (8, 9, 29). Indeed, thanks to the combination of HP 3He, phase-contrast (PC) velocity mapping, and rapid imaging, full 2D maps of the airway velocity vector can be obtained within a few seconds.
Here, we directly compare, for the same realistic patient-based geometry, velocity patterns obtained experimentally with HP 3He MR velocimetry to those obtained by CFD simulations. On the one hand, a realistic airway cast, issued from thoracic CT scans on a healthy patient was built using rapid prototyping and velocity maps were obtained by HP 3He MR velocimetry at different cross sections (CS) down to the fourth generation in that physical airway model, conveying a steady flow. On the other hand, CFD simulations were carried out in the corresponding digital airway geometry using the same dynamic conditions as for the experiments. Finally, numerical and experimental velocity maps were qualitatively and quantitatively compared.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Details on the segmentation procedure were given in previous papers (16, 44). Images were segmented using a method based on quantitative 3D CT bronchography (14). This fully automated 3D reconstruction approach allows the segmentation of the bronchial tree down to the fifth to seventh generations and the extraction of a central axis that provided navigation and quantification functionalities.
To obtain the airway geometrical description composed of triangular meshed lateral surface from the trachea down to generations 5–7, a marching cube algorithm (25), used to render iso-surfaces in volumetric data, was implemented in AMIRA (Mercury Computer System, Chelmsford, MA). After a filtering and a reduction of the number of triangular cells, the mean segment size was 
1.2 mm, which was adapted to match mesh quality requirements for CFD.
A tetrahedral mesh was then obtained using standard mesh generation algorithm in Tgrid (Fluent, Lebanon, NH). The total number of generated nodes was 203,196. The mean tetrahedral volume was 0.166 mm3 and mean segment sizes were 1.19 mm. The total number of cells was 106 (1,024,713) covering an airway lumen of 169.9 cm3. For comparison, the tetrahedral volume was twice lower (56%) than the voxel volume on the CT images. The model contained 58 end branches, and mesh density, which was uniform throughout the cast, was chosen to have a sufficient number of cells in the end branches. The specific quality criterion was skewness, i.e., the difference between the perfectly equilateral tetrahedron and the considered tetrahedron volumes divided by the volume of the perfectly equilateral tetrahedron. In that case, the perfect equilateral cell has the same circumradius as the considered cell. According to this definition, a skewness of 0 corresponds to perfectly equilateral cells, whereas a skewness of 1 reveals completely degenerate cells. A second simulation on a finer mesh (26% additional nodes and 35% additional cells) was also performed to investigate the local convergence of the CFD algorithm. For both meshes, skewness ranged from 0 to 0.8 with a mean value of 0.5. Maximum skewness was chosen to fulfill the commercial solver manufacturer (FIDAP, Fluent) recommendation to limit discretization errors and convergence time.
The triangular surface mesh was then converted into volume by giving a 1.5 mm thickness in the external direction of the bronchi by extruding the initial mesh with 3DsMax (Autodesk, San Rafael, CA) and by manually correcting the error in the generated volume. From this volumetric description of the surface, the cast was materialized with a scale of 1:1 using rapid prototyping (stereolithography). This method consisted of fixing a photosensitive resin with a laser, layer by layer (thickness: 0.1 mm), to generate a 3D object of the inner airway geometry.
Gas production and administration. MR measurements were performed using HP 3He as a high MR gaseous tracer, which was mixed in a vector gas (N2) and injected through the airway cast while the imaging sequence was running. This technique was previously described and validated in straight and U-shaped pipes (9). 3He was hyperpolarized on site (4) for 1 h to get two doses (140 ml polarized at 10%). Nondepolarizing vector gas, 650 ml of N2, was added to 70 ml of HP 3He before the mixture was transferred to a Tedlar bag for transportation to the magnet located 15 m away.
An administration system (Fig. 1) was designed to automate the injection of the HP gas mixture from the bag into the cast, continuously measure the input flow rate, and trigger the imaging sequence. This system was MR compatible (no magnetic materials, no radiofrequency (RF) noise production, low sensitivity to power RF emissions), and adapted to HP gases (nondepolarizing materials to maintain the prepared nuclear polarization).
|
Flow rate was measured with a nonmagnetic airflow sensor (AWM700, Honeywell International, Morristown, NJ). This sensor was made of plastic materials, contained in a shielded box, and powered with a nonmagnetic battery. The measurement principle consisted in the creation of a pressure drop between the two sides of a honeycomb. One percent of the total flow went through the bypass line, where the flow rate was measured via a thermal method.
Error sources in flow measurement were carefully analyzed. This flowmeter provided an approximate linear response for airflow up to 420 ml/s. It was calibrated with a 3-liter syringe (Hans Rudolph, Kansas City, MO) with air at room temperature and pressure. This calibration showed to be reproducible within 1%. It is known that the resulting pressure drop throughout this type of flowmeter depends on the velocity profile and disturbances to flow associated with flowmeter geometry and response may vary by a few percent (15). Response of the sensor was proportional to its voltage supply provided by the MR compatible battery, which decreased by no more than 3.7% while used over an entire day. Furthermore, while the gas physical properties transiently changed (i.e., density and viscosity) secondary to the injection of HP gas over a limited period of time during the experiment, the sensor, calibrated on air, underestimated flow rate by <2.8%. Considering the hyperpolarized gas production setup, able to produce two doses in 1 h, and its limited ability to keep the polarization for a long time, experiments on different CS were done separately. Consequently, flow rate was not maintained constant during gas production, but was set for a few minutes before each experiment with the error sources described above. Henceforth, the overall precision of the input flow rate between experiments on different CS could be estimated to be 7.5%, which does not differ from SEs performed with this type of flowmeter and administration protocol (13, 15, 31).
Valve switching and sequence triggering were controlled by an in-house program using HPVEE (Hewlett Packard, Loveland, CO) and flow rate was continuously acquired (DT VPI, Data Translation, Marlboro, MA) and stored at the end of each experiment for comparison with the imaging data.
All connecting tubes were taken from respirator hardware (standard flexible 22-mm diameter plastic tubes, Flextube #1573, Intersurgical, Wokingham, UK). An adaptation piece was placed at the cast entrance. It was composed of first 10-cm-long, 1-mm-diameter thin capillaries placed to fill entirely the 22-mm-diameter tube, and then, a respiratory filter directly connected to the cast. This stabilization segment aimed at reducing the curvature influence on the velocity profile at the entrance of the cast leading to a rather flat velocity profile, as checked by the velocity profile measurement performed at this location.
The airway cast laid horizontally in the MRI system (Fig. 2). Its position was adjusted so that the explored CS was centered in the magnet bore to reduce Maxwell-term artifact in the velocity measurement (33) and orientation was set normal to the static magnetic field B0. Positioning errors were estimated to be smaller than 1 mm and 2° for the angular orientation.
|
exp = 1.164 kg/m and viscosity 
exp = 17.9 µPa·s. The Reynolds number (Re) was 770 ± 58 (mean ± SD) at the tracheal inlet where the SD depended on the uncertainty on the input flow measurement.
Imaging technique.
Experiments were performed on a whole body 1.5 T MRI unit (GE Signa, Milwaukee, WI; inside tunnel diameter 55 cm) with a maximum gradient strength of 22 mT/m and a slew rate of 120 T·m–1·s–1. A transmit/receive switch and radio-frequency preamplifier as well as an 8-cm diameter surface coil (quality factor = 300) were specially designed and manufactured for 3He imaging (resonance frequency at 1.5 T = 48.65 MHz).
The sequence was based on the combination of radial MR imaging and phase-contrast velocimetry. A 2D radial-based PC spoiled gradient-echo sequence was implemented (see Fig. 3). It relies on the acquisition of a set of diameters in k-space (i.e., Fourier transform of projection of the imaged object) and on the reconstruction of the image from a sufficient set of projections separated by uniform angles. Each projection was realized with first a spatially selective excitation (1-cm-thick slice, flip angle 10–45°) followed by a bipolar gradient [field of speed (FOS) = 6 m/s, corresponding to a 2
phase shift], a readout gradient with a 1.56 mm spatial resolution [field of view (FOV) = 10 cm, 64-point full-echo readout with a 16-kHz bandwidth, echo time (TE) = 6 ms] and finally a spoiler gradient. This sequence pattern lasted TR = 12 ms and was repeated with different phase-contrast encoding gradients to get the three velocity components with a four-point measurement (37) before changing the orientation angle of the projection using a dichotomous algorithm. For each CS, a total of 1,024 projections (256 per velocity encoding steps) were continuously acquired during HP gas flow, leading to a total scan time (Tacq) = 12.3 s. The order used for projections allowed to reconstruct a posteriori, with complex filtered back projection, dynamic images of the flowing gas. The technique allowed playing on different aspects of the reconstruction to favor either temporal or spatial resolution. A trade-off of 32 projections per reconstructed images was used. It led to an effective total scan time per image of 1.5 s for the four velocity steps. From the four velocity-encoded images, a velocity map of each component was calculated by taking the phase difference between the corresponding encoded image and the reference image. The corresponding velocity-error map was calculated from the signal intensity image. Noise SD 
was calculated in a region void of signal and the reference image was normalized by this noise. The velocity uncertainty v for each voxel containing the signal S was then deduced using the following relation describing the statistical error on the measured phase (i.e., velocity) resulting from the measurement process (11):
 | (1) |
|
|
Airflow was simulated with FIDAP using a Linux 64-bit workstation, with a single AMD Opteron 2.5 GHz processor and 8 Go RAM. A segregated-solver algorithm with relaxation and upwinding was used. The tetrahedral elements were linear for both velocity vector and pressure, and hence pressure stabilization was used (19). Galerkin formula for intern laminar calculation and standard finite element stabilization methods were used (streamline upwind Petrov-Galerkin and pressure stabilization Petrov-Galerkin). The overall criterion for convergence was met when the solution-vector relative-error norm was <10–5. To quantify local convergence, comparison was made between simulations with the smaller cell size. Velocities were interpolated on a structured grid on all CS and a Bland-Altman analysis (2) was performed. Bias and SD of the difference (CFD) on axial velocity were evaluated and divided by the mean CS velocity to get relative statistical parameters.
Experimental and numerical comparisons. Experimental cross sections were registered with a 0.4-mm precision using contours from the numerical data, and mean signal-to-noise ratio (SNR) was calculated (Fig. 4 and Table 1). Numerical results were interpolated on a structured grid at the experimental cross-section locations (bilinear interpolation). Flow rates were assessed both for experiment and simulation by integrating the axial velocity over the CS. Flow rate error was estimated by propagating the measured velocity errors (usually 1–2% on the MR-imaged flow) and by adding the error obtained on the input flow measured by the flowmeter. Mean secondary motion amplitude and mean axial velocity on each section were calculated. Position and magnitude of selected flow features, such as CS axial-velocity maximum and vortex center positions were manually and independently selected on each velocity map and point-by-point compared.
|
diff) of the axial velocity difference was calculated to get statistical parameters that describe how well CFD and MRI velocities matched. |
RESULTS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
100 iterations in a calculation time of 7.5 h. At every CS, no bias on velocity was observed between the simulations performed with different cell sizes. CFD ranged from 0.7 cm/s (CS2,5,7) to 3.3 cm/s (CS9), with a mean value of 1.2 cm/s (3.7% of the mean velocity) for all CS, indicating that the cell size was small enough on each considered CS (Fig. 5).
|
Overview of experimental and numerical velocity comparison.
SNR was given for each CS. The mean SNR over the CS was 38, corresponding to a typical error on velocity of 3.5 cm/s (from Eq. 1). The measured and simulated mean axial and secondary velocities are given in Table 1, as well as diff for axial velocity. Mean axial velocity directly reflects flow rate and consequently matched for all sections except CS9. Mean secondary velocity tended to be overestimated by MRI, thus the ratio of mean secondary velocity to mean velocity was only calculated for CFD. This ratio varied from 15 (CS7) to 45% (CS6). It tended to increase before flow dividers (CS2,5,9) and after curvatures (CS3,6). From the Bland-Altman analysis that provided parameters to determine how well CFD and MRI corresponded for a given CS, no significant bias was observed, except for CS9. Mean bias was between 0.1 cm/s (CS7) and 11.8 cm/s (CS9), with a mean value of 1.4 cm/s. The SD diff varied from 6 cm/s (CS6) to 32 cm/s (CS9), with a mean value of 14 cm/s. Figure 6 plots MRI vs. CFD axial velocity for CS7. Most of the data ranged within the 95% confidence interval, confirming the good agreement between numerical and experimental results.
|
|
|
|
|
As found by both CFD and MRI, velocities within CS2 were smaller than those within CS1 due to its larger surface area. The small backflow region that appeared on the simulated data was not observed on the experimental data. In both cases, the velocity profiles indicated the onset of flow division before the fluid reaches the geometrical bifurcation. Furthermore, in-plane velocity exhibited two main directions corresponding to the downstream bronchus orientations. Maximum in-plane velocity was near 0.2 m/s both for CFD and MRI.
Flow in the left main bronchus. For CS3, the axial velocity profile was highly asymmetrical with a maximum velocity toward the inner edge of the first bifurcation and a minimum toward the outer edge. Despite higher numerical velocities, CFD and MRI agreed on this asymmetry. Secondary motions exhibited similar patterns with comparable velocity values. Two eddies [originally described by Olson, quoted by Pedley (36)] developed in the airway with circulation oriented from the minimum (i.e., the inner edge of the curvature, outer edge of preceding bifurcation) toward the maximum axial velocity (outer edge) and skirting along the edge. Due to the asymmetry of the section and bifurcation, eddies were asymmetrical and of unequal strength, the posterior eddy having, in our model, much higher intensities and spatial occupancy.
Flow in the right main bronchus. For CS4, the axial velocity profile was highly asymmetrical with a maximum displaced toward the inner edge of the first (carinal) bronchial bifurcation and a minimum toward the outer edge. Normal to the bifurcation plane, an M-shaped profile was seen along the anteroposterior (AP) diameter (Fig. 8). The peak axial velocity was higher for CFD than for MRI. Two eddies developed with the circulation oriented from the outer to the inner edge of bifurcation and then skirted along the edges, as expected for a symmetrical bifurcation. Similar to CS3, the posterior eddy had higher intensities and extended more than the anterior eddy. Vortex intensities and location of the eddy centers agreed between CFD and MRI.
CS5, located one diameter downstream from CS4 exhibited a profile similar to that of CS4 with a maximum axial velocity located toward the inner edge of the first bifurcation. The axial profile indicated the onset of flow division before the fluid reached the downstream bifurcation (Fig. 8). Lower velocities were observed in the center of the section and a second peak toward the outer edge of the first bifurcation.
Flow in the upper right bronchus. Velocity intensities were similar for CFD and MRI, lower than the values encountered upstream on CS4–5. Simulated axial velocity exhibited an asymmetrical shape with a peak (0.35 m/s in both cases) directed toward the outer edge of the second bronchial bifurcation, i.e., the outer edge of the first bifurcation, thus a continuation of CS5. This result was unexpected, as previous studies observed a higher influence of the immediately upstream bifurcation with maximum intensity toward the inner edge (3). This was probably due to the interaction between bifurcations, because this section was fed by a nonnegligible part of secondary motions visible at CS5.
Flow in the intermediate bronchus. The intermediate bronchus branching angle was 15°, nearly providing a straight path from the right main bronchus. Consequently, the CS7 velocity pattern was a continuation of that of CS4–5. The velocity profile was highly asymmetrical, with a maximum velocity near the outer edge of the upstream bifurcation and a minimum near the inner edge. Two eddies developed with the circulation oriented from the minimum to the maximum axial velocity. The eddies were asymmetrical, the posterior one being much stronger. A two-peak axial velocity profile right-left (RL) and a pronounced M-shape profile in the plane normal to the bifurcation plane (AP) was observed (Fig. 8). Axial as well as in-plane velocities agreed fairly well (Fig. 9A). The location of maximal velocities matched between CFD and MRI.
CS8, located one diameter downstream from CS7, exhibited the same type of asymmetrical velocity pattern in relation with unequal eddies. Flow development was illustrated by a more homogeneous axial velocity profile. We emphasize two findings: secondary motion intensity increased up to 21% (15% for CS7), and axial velocity peak was shifted toward the posterior region, the latter likely being the consequence of a predominant posterior eddy visible at CS7. This shift was precisely confirmed by the experiments (Table 2).
Flow in the lower right bronchus. CS9 was located in a fourth generation branch, slightly upstream from a fifth-order bifurcation. Due to the small branch length (<2 cm) and a small cross-section area, the experimental flow and velocity measurement suffered from partial volume effects. Nevertheless, velocity maps still gave similar results: two axial velocity peaks could be observed with higher values than encountered upstream (0.95 m/s) and the velocity profile started to divide, illustrating the influence of the fifth order bifurcation.
|
DISCUSSION |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
In the presented geometrical model, left-right flow distributions were unbalanced in favor of the right lung as expected in a normal lung and previously observed in asymmetrical idealistic models (55–45%; Refs. 3, 18, 24). The flow rate measured at the fourth generation reached 23% of the inlet flow and axial velocities were high. Approximately the same proportion (16%) was obtained at the same location in a post mortem airway cast composed of 31 generations (18). The tendency of the gas to follow the straightforward flow path (an effect related to gas inertia) may explain this previously described feature.
As shown in Figs. 7 and 8, experimental peak velocities exhibited lower intensities and blunter velocity profiles than those found in the numerical simulations. The inlet velocity profiles used for simulations were not simultaneously imaged with the selected CS. Therefore, one might expect a slight variation of the experimental inlet conditions when the cast was repositioned for each CS. However, the two successive upstream bifurcations effectively reduce the influence of the inlet flow conditions; furthermore, as observed in Figs. 7–9, velocity profiles matched better on CS7 and 8 than on upstream CS. Except for CS9, numerical mean axial velocities fit experimental values but mean secondary motions slightly differed (Table 1). The latter discrepancy could be due to an error on the imaged slice position, which can induce a small angular deviation with respect to the airway axis, thus projecting the axial velocity onto the transverse measurement plane. Thus, despite slight differences in the tracheal region of the model, understandable on the basis of slight discrepancies in entry conditions, we find similar flow features such as velocity extrema and vortex centers in similar regions (Figs. 7 and 9). Local comparison on Table 2 showed that extremum magnitudes were similar, whereas their locations were almost superimposed. Point-by-point quantitative comparison presented in Fig. 9 confirmed the good agreement between numerical and experimental velocities.
Numerical and experimental investigations revealed flow characteristics that may have important impact on airway resistance, particle deposition, and gas mixing. In the main bronchi, the axial velocity peak was along the inner side of the bifurcation and secondary motions form Dean motions (3, 24) with magnitudes higher than those found in previous studies. In the lower right subsegmental bronchus, due to the small branching angle, the upstream bifurcation had less influence in the velocity distribution than in the first (carinal) bifurcation. This effect has already been simulated in a similar configuration (6). CS7 and 8 presented similar velocity profiles but with a shift in the velocity peak location and a development of the dorsal clockwise secondary motion from CS7 to CS8. The smaller cross-sectional area and shape variation of CS8 resulted in slightly higher mean axial velocity component and secondary motion than those on CS7, whereas a simple flow development would be expected in an idealized bifurcation (36). Secondary motions appear in the straight tubes with noncircular section (rectangular and triangular; Ref. 38) and are associated with section changes (divergence; Ref. 45). The present observation of higher secondary motions suggests that gas mixing, characterized by vortex intensities, increases in realistic geometry not only due to the branching geometry and asymmetry, but also with size reduction along branches, noncircular sections, and probably local branch curvature.
Experimental limitations. MR gas velocimetry used to measure 3D flows in complex geometries was validated in a previous in vitro work by our group and applied in vivo (9). As most experimental velocity measurement techniques, MR velocimetry suffered from higher uncertainty in the low velocity range so experimental errors increased near the airway walls. This is mainly due to the combination of partial volume (41) and time-of-flight (32) effects. Consequently, mean secondary motion quantification given in Table 1 was done without considering data at the edges because they artificially increased discrepancy between CFD and MRI. Except near the edge, comparison of local velocities (Fig. 9) showed that in-plane CFD and MRI velocities matched fairly well.
Furthermore, the velocity error due to the effect of averaging over the nonnegligible size of the measurement volume (the voxel was 10 mm thick) was numerically evaluated. It was found to be smaller than the experimental error bars for CS1–8. Therefore, this smoothing effect could be neglected. This approximation, which is related to this measurement technique, was previously validated in other reference geometries (e.g., curved tube; Ref. 9). For CS9, the branch length was on the order of the slice thickness, resulting in nonnegligible averaging effects that might contribute to the discrepancy between MRI and CFD, especially for flow rate. Another source of error in the quantitative comparison concerns the position of the experimental CS. A spatial 1-mm, 2° angular precision and 0.4-mm registration could easily introduce an additional term in the calculated error diff.
Gas properties. The steady-flow condition that was studied herein, has been previously considered to be a fair approximation of inspiratory peak-flow conditions (27). Helmholtz number was null during the experiments, a reasonable approximation to in vivo resting conditions. For the CFD simulation, the fluid was assumed incompressible (Mach number << 1, as during normal breathing), and Newtonian, a common assumption for airway gas. In addition, because the Re was <1,000 in the airway model, airway flow was simulated in laminar conditions.
Geometry. The airway model used for both experiments and simulations was highly realistic as it was generated from images of a human subject. As in most previous in vitro studies, branches were considered rigid and surfaces dry, an assumption that has been discussed in an earlier paper (3). Considering that the CT data were acquired during apnea, no major tracheal or bronchial deformation can be expected [contrary to substantial deformation or even airway collapses, which are currently observed during forced expiration (43) and deep breathing]. As a consequence, the model's geometry reflects resting conditions, but, in addition to previous studies based on more idealistic geometries, our model includes several realistic features: 1) 3D curvatures and unequal lengths of airway branches, 2) change in CS shapes and areas along the airway tree, and 3) asymmetrical airway junctions and variety of branching angles. Due to the spatial resolution of CT imaging data, five to seven generations were included in the airway model. Note that this cast built by stereolithography accurately reproduced the highly complex airway geometry to within 0.1 mm of the numerical geometry. It was noteworthy that no similarity condition was necessary to reproduce breathing conditions, as the studied airway geometry was life size. This was not the case in previous studies where enlarged airway models were used to permit flow measurement (3).
Inlet boundary conditions. The airway cast was repositioned to align each CS on a sagittal plane because the MR pulse sequence presently used did not allow oblique plane imaging and to minimize Maxwell term artifacts (33). The inlet velocity profile was thus measured at the exit of a sharply curved tube necessitated by the spatial limitation of the magnet bore. However, thanks to the smoothing segment, a flat velocity pattern was observed confirming its efficiency (Fig. 7).
The comparison between CFD and MRI was done here by considering the same inlet geometrical features. However, upper airways were not included in the model, and thus important geometrical features at the inlet, such as glottis and larynx, were lacking. These features could have an impact on the flow patterns in the trachea, including the generation of significant 3D gas motions, as shown by both in vitro (3, 28, 35) and in vivo (9) studies. Additional studies including upper airways at the entrance of the present model could be advantageously performed to address this limitation (7). Despite this, the velocity profiles after the first bronchial bifurcation might only be moderately affected by the inlet profile (28).
Outlet boundary condition. Our five to seven-generation airway model accounts for most of airway resistance (36), as maximum airway resistance is known to occur up to the fourth generation at resting conditions. CFD outlet boundary conditions were identical to those experienced in MRI experiments where all outlets were opened (zero relative pressure). Physiological conditions might have to include the effect of smaller bronchi through the addition of dissipative pressure drop and gas distribution. To the best of our knowledge, no in vivo measurement of resistance nor pressure at that level has been done due to the lack of experimental technique. Airway resistance might correctly be modeled by Poiseuille flow beyond the tenth generation (36). Even if a Poiseuille resistance seems relevant to predict resistance or flow distribution in small airway (18), in vivo, asymmetric action of the respiratory muscles or heterogeneous physiological or pathological airway obstruction might lead to an infinite number of possibilities for the outlet boundary conditions. However, in a flow distribution study in a six-generation model (20), it was shown that branch obstruction has a limited spatial effect, again because local geometry predominates in rapidly changing branching systems, meaning that exit conditions would not deeply modify the normalized shape of the velocity profile.
Reynolds number.
The relative SD over the input flow rate was estimated to be 7.5%, which is probably overestimated given the sharp correspondence of only 3% between CFD and MRI. Considering for an adult at rest a tidal volume of 500 ml, a breathing frequency of 12 cycles/min, an effective inspiration lasting one-third of the cycle duration, and a tracheal mean diameter 2 cm, inspiration flow rate is estimated 300 ml/s and Re 1,200. Here, the considered Re and flow rate were slightly lower than this estimation but remains characteristic of the convective flow regimen encountered in this case. Furthermore, results might be extended to higher Re, as in various in vitro studies, it was shown that both flow distribution (20) and flow patterns (3) depend more on the geometrical airway structure than on Re.
In conclusion, the work presented in this paper concerns the application of an original gas velocity mapping technique. We demonstrated that velocity data can be obtained on the same realistic airway geometry by coupling experimental flow measurements and numerical flow computation. The HP 3He MRI phase-contrast velocimetry was applied for the first time in a realistic human airway model and agreed well with the numerical simulation. The experimental technique thus provides a powerful tool to map gas velocity in complex and even opaque structures, including human bodies.
The data confirmed the complexity of flow and the strong geometry dependence of convective flow patterns within large airways. Expected global flow patterns were observed as a result of the global branching structure, but local flow patterns strongly depended on local geometry as previously demonstrated on the basis of the comparison of data with literature results (44). Because of its small branching angle, the right lower subsegmental bronchus exhibited a flow that is mainly influenced by the geometry of the first bifurcation (including the carina). Secondary motions were higher than in previous studies due to the high noncircularity of CS, large area variations as distance from entry increases, and most likely 3D curvatures of real airways. These results suggest that gas mixing is enhanced in realistic airway morphologies compared with idealistic models.
The present experimental results elegantly validate a full computational process made of thoracic CT acquisitions, specific airway segmentation, and CFD on a patient-based geometry. The final item is a morphofunctional simulator of airway flow serving as a new numerical tool to study respiratory airflow and potentially physiological and pathophysiological flow-related mechanism.
The current flow model simulates resting conditions and the tool might be used to probe the effect of different gas properties. As flow patterns are mainly dependent on the patient-specific airway geometry, it opens a route for personalized and spatially localized diagnosis of respiratory diseases. The functional consequences of geometrical modifications induced by diseases (e.g., airway tumor) or by surgical interventions could also be assessed. Patient-based simulations may also provide a useful tool to optimize respiratory assistance and help the development of inhaled drug therapies.
|
GRANTS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
ACKNOWLEDGMENTS |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
The first two authors contributed equally to this work.
|
FOOTNOTES |
|---|
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
|
REFERENCES |
|---|
|
TOPABSTRACTMATERIALS AND METHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
| ||||||||||
, a bipolar gradient encoding velocity was applied just before the radial imaging readout gradients Gread. Finally, the tipped magnetization was spoiled. This sequence lasted for TR = 12 ms. Echo time (TE) and signal observation time (Tobs) are also indicated. The sequence was repeated 4 times for the 4 velocity encoding steps before varying the readout gradient orientation 
. The total imaging time for each CS was Tacq = 12.3 s.
