INTRODUCTION

Top Abstract Introduction Background Methods Discussion References

IMAGING TECHNIQUES are becoming increasingly important in physiology because of the need for noninvasive determination of regional structure and function in intact subjects. Xenon-enhanced computed tomography (Xe-CT) is one such method for the measurement of regional pulmonary ventilation, determined from the washin and washout rates of the radiodense nonradioactive gas Xe, as measured in serial computed tomography (CT) scans. Whereas the feasibility of this technique was first demonstrated many years ago (7, 8), advances in CT technology now make practical the quantitative determination of regional ventilation distribution in the lung on a length scale as small as several millimeters. Combined with the unique capability of CT to describe detailed anatomy (26) and pulmonary perfusion (11), this single-imaging modality can thus create a nearly complete noninvasive structural and functional characterization of the lung. However, along with these advances in noninvasive measurement techniques must also come methods for the estimation of variability and confidence intervals (CI) so that the statistical significance of the information obtained may be evaluated. This problem is particularly important for CT-based measurements, since obtaining repeated measurements is frequently impractical because of the time and expense involved.

Monte Carlo (MC) techniques have been successfully applied in physiology for determining confidence limits and the sensitivity of model parameters to noise and for performing modeling and simulations that include realistic biological variability (1, 6). In this paper, MC techniques are used to estimate confidence limits for regional ventilation measurements made by using the Xe-CT method, and these estimates are compared with values determined by repeated measurements. The MC approach is evaluated in terms of its unbiasedness and coverage of the CI. Finally, MC simulation is used to compare three different experimental protocols and parameter estimation techniques for Xe-CT studies.

	BACKGROUND AND MODEL

Top Abstract Introduction Background Methods Discussion References

Xe-CT ventilation measurement. Because stable Xe gas is denser than air, its concentration may be measured by CT and is linearly related to its CT enhancement. Repeat CT scans taken at constant lung volume (usually at end expiration) during the washin and/or washout of 30-70% Xe yield exponential local density curves for a chosen lung region of interest (ROI). The time constant () of this curve, determined from a nonlinear curve-fitting procedure, is equal to the inverse of the local ventilation per unit volume (specific ventilation) if single-compartment behavior is assumed. By selecting different ROIs, spatial patterns of ventilation may be analyzed. Ventilation to the entire lung may be obtained by placing a ROI over the trachea. A given slice may be subdivided and analyzed according to any desired regional scheme, and, by combining the results from multiple slice locations from apex to base, the total spatial distribution of ventilation throughout the lung can be characterized. Thus this technique provides a direct noninvasive measurement of regional ventilation, with precise anatomic correlation provided by the CT image.

The following is a description of a simple washin (wi) protocol, although the general principles apply to washout (wo) and combined washin/washout (wi/wo) protocols as well. The subject is positioned in the CT scanner, and the level of the desired imaging run is determined from a scout image. CT images at the same location are obtained after each breath at end expiration. Typical settings for the GE9800 CT scanner are 80 kV, 120 mA, 10-mm slice thickness, 2-s scan, and 512 × 512 pixel image. It is not necessary that mechanical ventilation be used, but the same end-expiratory volume must be achieved for each image. Two or three baseline images while subjects are breathing room air or air-O₂ mixture are first acquired, and then the inspired gas is changed to the Xe-O₂ mixture. Typically, 50% Xe-40% O₂ in air is used for animal studies to maximize the change in density; lower Xe concentrations may also be used. Sequential end-expiratory images are acquired until the lung is equilibrated with the Xe-O₂ mixture, usually 15-30 breaths later. The process is then repeated at other lung locations or after various interventions. Variations on this basic protocol include imaging some subset of breaths, continued imaging as the tracer gas is washed out of the lung (wi/wo protocol), or imaging only the washout phase after the lung is first equilibrated with tracer gas (wo protocol). Issues related to the choice of a particular protocol are considered below (see DISCUSSION).

The images are analyzed by selecting a desired ROI, which may range from the entire lung field visible in the cross-sectional CT slice to an arbitrary region as small as 0.5 cm² drawn on the image. The mean density in this ROI is measured in each image and plotted as a function of time or image number. The time constant  is then determined by fitting this curve to a single-compartment exponential model by using a nonlinear least squares curve-fitting procedure. Example ideal wi, wo, and wi/wo curves along with their respective model equations and parameters are depicted in Fig. 1. Depending on the  relative to the number of images obtained, full equilibration may or may not be achieved (particularly in the wi/wo protocol). Note that in all three models D₀ is used to denote the baseline density value (without Xe), whereas D_f refers to the density after full equilibration with Xe.

View larger version (22K): [in this window] [in a new window]   Fig. 1.   Example washin (wi), washout (wo), and washin/washout (wi/wo) curves and their respective model equations and parameters. See text for further explanations and definitions.

Stochastic model. In reality, of course, the curves are not perfectly smooth but contain noise (Fig. 2). There are several sources of this variability. Measurement error, in terms of the accuracy of the CT imaging system, is probably only a very small component of the total variability (25). Most of the variability comes from the fact that CT measures density, and one cannot distinguish between changes in density due to the tracer gas (Xe) vs. the underlying substrate (lung tissue). Lung density will change with small changes in lung volume, registration of the identical ROI from image to image, and changes in blood volume within the ROI. The action of the heart, which both moves the adjacent lung tissue and changes the instantaneous pulmonary blood volume, also contributes to this artifact. This variability should add a random fluctuation to the actual density value, and thus it represents a stochastic component that can be explicitly added to the deterministic model and then utilized for the statistical analysis. Thus we chose to add a stochastic term _t to the baseline density D₀ at every point to simulate the noise component of CT-derived Xe tracer curves. Analysis of the residuals of actual fitted data by using normal probability plots (Fig. 3) and the Shapiro-Wilks test for normality (Stata, 4.0, College Station, TX) showed that the noise was normally distributed, so _t was modeled as a normally distributed variable with zero mean and variance ² [_t ~ N(0,²)]. For the washin model, this is represented as

View larger version (13K): [in this window] [in a new window]   Fig. 2.   Example wi/wo computed tomography (CT) density data () and fitted curve (solid line) from left lung of an anesthetized dog. Residuals for each point are plotted above. Curve fitting resulted in a time constant () of 4.75 and normalized summed squared residuals (SSR) of 5.5. HU, Hounsfield units.

View larger version (17K): [in this window] [in a new window]   Fig. 3.   Normal probability plot of representative wi/wo curve fit residuals (xi), demonstrating that they are normally distributed. P = 0.338 for Shapiro-Wilks normality test applied to this data set. Normal F[(xi  µ)/], normal cumulative probability of xi with mean (µ) and SD (); P(i), empirically observed cumulative probability of xi in ascending order; n, no. of data points.

	ANALYTIC METHODS AND RESULTS

Top Abstract Introduction Background Methods Discussion References

Parameter estimation. Parameters D₀, D_f, and  are estimated from each set of data by using a nonlinear least squares curve-fitting procedure (Igor Pro, Wavemetrics, Lake Oswego, OR). The starting points of the washin and/or washout segments (t₀, t₁) are initially estimated by visual inspection of the curves and then determined by the curve-fitting procedure. Although data points occur once per breath, t₀ and t₁ are permitted to take on fractional values, reflecting that the arrival of the front of Xe to the lung periphery may occur in midcycle because of the variable amount of dead space relative to tidal volume in the system. Wi/wo fits are constrained to fit the same  for both wi and wo segments. The units of density used are Hounsfield units (HU), offset by 1,000 so that air has a density of 0 and water 1,000 (instead of 1,000 and 0 as with standard HU). Because images are obtained once per breath, the time constant  is conveniently given in units of breaths, which may be converted to time units by dividing by the respiratory frequency.

The variance ² of the stochastic noise parameter _t was determined from the normalized summed squared residuals (SSR) of a given fitted data set where D(t_i) is the density measurement at t_i, (t_i) is the corresponding value of the fitted curve, and n is the number of data points. Simulated wi/wo curves generated by adding this "matched" noise to a curve generated from the parameters determined from experimental curves are similar to a set of curves obtained by repetition of the entire imaging study the same number of times (Fig. 4). Note that in Fig. 4A there is some baseline variation between the experimental curves, which is not modeled by this process. Ignoring this baseline variation between curves does not affect the estimation or CIs for .

View larger version (20K): [in this window] [in a new window]   Fig. 4.   Comparison of 10 repeated experimental (A) and simulated (B) wi/wo lung density curves for region of interest (ROI) 1 of Table 1. Note that baseline shift in A is real and represents changes in baseline lung density (D0) over the 90 min required to obtain 10 repeated studies. Monte Carlo model used mean D0 and density after full equilibration with Xe (Df) from A for all simulated curves.

MC CI determination. CIs for a given set of Xe-CT data are determined as follows. A set of simulated curves is generated by adding normally distributed random noise to a wi, wo, or wi/wo curve generated from the actual parameters determined from the fitted curve. The variance of this applied noise is set equal to the normalized SSR of the fitted curve. Each of these simulated curves is fitted, the time constants _i are determined, and the 95% CI is obtained directly from the distribution of _i. CI for the other fitted parameters D₀ and D_f may also be obtained from the same procedure.

Evaluation of MC CI method. MC CI estimates for all three fitted parameters were determined for several data sets at several noise levels and by using 10, 100, 500, and 1,000 repetitions. There was minimal change in the CI estimated above 100 repetitions, and thus all results reported here used 100 repetitions unless noted otherwise. This approach was evaluated by comparison of MC CI estimates with CI determined from repeated experimental measurements and by the evaluation of the "unbiasedness" and "coverage" of the parameter estimates and CI (5).

Ten identical, repeated Xe-CT wi/wo ventilation studies were performed in an anesthetized, paralyzed, and mechanically ventilated dog, and the  values were determined for several ROIs in the left lung for each run (Fig. 4). The 95% CI for these repeater measurements, estimated as the mean ± 2 × SE, was compared with the 95% CI estimates for the same data when using the MC approach (Table 1). ROI chosen included the entire left lung field on the CT slice (ROI 1), the arbitrarily divided top and bottom halves of the left lung field (ROI 2 and 3), and arbitrary circular regions from within the top and bottom halves (ROI 4 and 5). ROI volumes ranged from 21.5 ml (ROI 1) to 0.67 ml (ROI 4). The mean SSR values from each of the repeater series were used to set the noise levels for the respective MC simulations. The results indicate that the MC technique provides a reasonable and conservative estimate of the CI determined from repeated measurements. This held true even for the very noisy ROI 4, which was <1 ml in volume and located adjacent to the ascending aorta.

A statistical estimator is considered unbiased if its long-run average approaches the true value of the parameter as the number of experiments increases. The unbiasedness of the parameter estimates was determined by calculating the mean ± SD of the estimators, , ₀, and _f for 500 repetitions of simulated curves with added noise. The analysis was repeated for noise (variance ²) levels of 4, 25, and 100 and for nominal  values of 1, 2, 5, 10, and 20 for each of the three models (Table 2). Other nominal parameters used were D₀ = 300, D_f = 400, and 40 points per curve. These density values are typical for well-expanded lung, and the  values span the range obtained in normal animal studies. The distribution of 411 normalized SSR (noise) values obtained in four dogs during analysis of the vertical distribution of ventilation in the prone and supine positions (14) is presented in Fig. 5, demonstrating that the typical range of SSR values in healthy animals is within this modeled range.

View larger version (21K): [in this window] [in a new window]   Fig. 5.   Frequency histogram of normalized SSR from ROI obtained in 32 separate wi/wo studies in 4 dogs in both prone and supine positions (n = 411), demonstrating range of noise values found in healthy animals. Curves with SSR >150 typically have readily identifiable artifacts, such as partial volume effects or extreme cardiac motion.

Table 2 demonstrates that the parameter estimates are unbiased, even at high levels of applied noise and large values of . Results are presented for wi/wo and wo models only; the wi results were nearly identical to the wo values, with the important difference that the SD values for D₀ and D_f were transposed. In general, the SD increases with the SD of the applied noise (), as expected. As  increases at a given noise level, the SD of the parameter estimates also increases, although not uniformly. As a result, the relative error (SD/mean) falls slightly as increases, reaching a minimum at  = 5 and then increasing beyond the initial level by  = 20. This reflects the amount of information available at fixed sampling frequency for fitting the exponential at the extremes of  and is discussed in more detail below (see Simulation study: model comparison). Note also that the  SD values for the wi/wo model are consistently smaller than the values for the wo (or wi) model. For the wi/wo and wi models, the SD of D₀ remains relatively constant as  increases, but the SD of D_f increases more than four- (wi/wo) or ninefold (wi). Again, this difference reflects that D₀ is determined from the initial baseline values independent of , whereas the estimation of D_f becomes more difficult if a plateau is not obtained. The results for the wo model are similar to those of the wi model, with the roles of D₀ and D_f reversed.

To determine whether the MC CI procedure provides adequate coverage, it must be demonstrated that the "true" parameter value actually lies within the calculated CI a sufficient percentage of the time. To do this, the following simulation was performed for the wi/wo model. Given a set of nominal parameters  = {, D₀, D_f} and noise level ², a simulated curve was generated and curve fitted to provided parameter and noise (SSR) estimates (, ²). The MC 95% CI was then determined as above for this fitted set of parameters , and this CI was compared with the original parameters . This procedure was repeated 500 times, and the percentage of CI that enclosed the original parameters was calculated. Table 3 gives the CI coverage for  at two nominal values of  and three noise levels. Note that the coverage remains the same as the noise increases: the CI width expands with greater noise levels to contain the parameter. The coverage of the CI for the D₀ and D_f parameters was 100% under all conditions. These results for the 95% CI estimates suggest that the MC CI procedure adequately covers the CI.

Simulation study: model comparison. A quantitative comparison of the three Xe-CT protocols was undertaken to determine whether there was an optimal protocol for ventilation measurement. It was hypothesized that, given the same nominal parameters, noise level, and number of images, the wi/wo model would yield better  estimates because most of the information in estimating an exponential is contained in the steep initial portion of the curve, and the wi/wo model has two such segments compared with only one for the wi and wo models. A simulation was performed by comparing the MC CI width for each model over a range of  values (1-20) and noise levels (² 4-100). The results are presented as one-half CI width normalized to the parameter value, a measure of relative error (Fig. 6). The one-half CI width is comparable to 2 × SE for determining differences between means. The wi/wo protocol is superior to either the wi or wo approach for determination of , with the wi/wo relative errors 65-70% smaller at the lower  values, falling to 50% at  = 20 and the highest noise level. Results for the wi and wo protocols were similar. Note that at very small  values the relative error increases, reflecting that there are few data points on the steep portion of the curve under these conditions. At the largest values of , the relative error also increases, reflecting that the curves do not fully plateau when  is large relative to the number of images, as discussed above. In addition, the relative errors for  increase roughly with  for all three models.

View larger version (22K): [in this window] [in a new window]   Fig. 6.   Results of model comparison simulation. Relative error, expressed as half confidence interval (CI)/mean, for parameter estimates of , D0, and Df for wi (), wo (), and wi/wo (×) protocols as functions of simulation noise level (2) and nominal .

As seen in Fig. 6, both D₀ and D_f have small relative errors for all three protocols at  values of 10 or less. At  = 20, however, the relative error for D₀ increases fivefold for the wo model and that for D_f increases comparably for the wi model, again reflecting the difficulty in estimating the plateau. The wi/wo model has a slightly larger CI for D_f than the wi model for  10, since there is less time defining the plateau, but does better at  = 20 and with higher noise levels, probably because of the better overall fit from the more precise estimate of .

	DISCUSSION

Top Abstract Introduction Background Methods Discussion References

The use of imaging techniques for physiological measurement permits the noninvasive determination of regional function with a high degree of anatomic correlation. The measurements themselves typically involve complex image-reconstruction algorithms, and correlations between measurements made within an image or from a series of images may be impossible to determine. In addition, these measurements are frequently time consuming to obtain, expensive, and may involve significant radiation or X-ray exposure. The use of MC methods provides an alternative approach for the determination of CI and, furthermore, allows the exploration of various analytic approaches that may be difficult or laborious to approach experimentally. In effect, the MC method substitutes a statistical experiment for multiple physical experiments.

In this paper, we apply MC techniques to the Xe-CT measurement of regional pulmonary ventilation. Xe-CT uses a series of density measurements made from a sequence of CT images taken at the same location during the washin and/or washout of radiodense Xe gas, which are then fitted to an explicit exponential model. The method presented assumes that there is a component of real variability in these measurements that is likely both technical and biological in nature. This variability is estimated for each data set from the curve fit and then used to scale the stochastic noise used in the MC CI estimate. Evaluation of the approach shows that the parameter estimates are unbiased, that there is appropriate coverage of the CI, and that the CI estimates are reasonable and conservative compared with those obtained from a repeated-measurements experiment.

Three different Xe-CT protocols are presented and evaluated by MC simulation. All three are based on the assumption that the washin and/or washout of the tracer Xe gas occurs with first-order (monoexponential) kinetics, and, as such, a deterministic model can be described for each protocol (Fig. 1). The assumption of first-order behavior for regional lung washout implies that the ROI is behaving as a single well-mixed compartment. This model has been used extensively for the measurement of ventilation with radioactive tracers (2, 21, 22, 24). In the presence of heterogeneity, however, the single-compartment model may be inappropriate and result in an underestimation of the actual effective ventilation of a region (18). For radioactive tracers, multicompartment behavior is evident as a departure from linearity of the semilog washout curve and the presence of elevated counts at the conclusion of the washout, and a number of correction schemes have been proposed (17, 18, 23). For Xe-CT studies, one would expect to see a pattern to the residuals from the curve fits in the presence of significant heterogeneity. No such pattern was evident in any of the curves examined from normal lung Xe-CT studies, and heterogeneity may be less of a problem when small ROIs are examined. However, as this technique is applied to diseased or injured lungs, it will be important to watch for this type of behavior and consider more complex multicompartment models.

Because the Xe-CT method has potential for use in patients (9, 10, 19), possible hazards to humans must be minimized. Xe is an anesthetic gas [~30% more potent than nitrous oxide (3)], and it is recommended to limit its concentration to <50% (9). Thus it is desirable to determine an imaging protocol that maximizes the quality of the data obtained while minimizing the exposure to radiation and inhaled Xe. Wi, wo, and combined wi/wo protocols have all been used with different imaging techniques. Wi techniques have certain advantages over wo, including minimal exposure to the tracer gas, decreased consumption of expensive tracer gas, decreased artifact due to tracer uptake by the blood (related to the mean concentration and time of exposure), and ease of use, but may be adversely affected by incomplete tracer equilibration. On the other hand, wo methods, in which the lung is fully equilibrated with tracer before the imaging of the tracer washout by fresh gas, have other advantages. The wo protocol can provide an index of lung gas volume from comparison of the equilibration image with the baseline and is also more sensitive for detecting regions with long s. Combined wi/wo protocols have been used with radioactive tracers in the lung (16, 17) and Xe-CT cerebral blood flow measurements (12, 15) to take advantage of shorter exposure times to the tracer gas while theoretically making up for the lack of full equilibration with the information obtained from the second-exponential phase. The MC simulation quantitatively confirms this advantage, demonstrating that the wi/wo approach provides a narrower CI for  for the same nominal parameters and noise levels than the wi or wo models. Other important protocol issues, such as the effect of different Xe concentrations and further refinement of the imaging protocol to minimize the number of images obtained, could also be explored by using this simulation technique.

The stochastic model chosen superimposes noise on the baseline density parameter of the deterministic model to create fluctuations in the simulated density measurement. Noise could have been added to the dynamic term of the model either alternatively or in addition to the baseline term with equivalent results. There are multiple real sources for this variability. Motion of the lung, which inflates and deflates between images, results in small differences in end-expiratory volume and changes in the exact tissue elements imaged in each ROI. This effect will be greater as slice thickness and ROI size decrease, because small changes in the amount of denser tissue elements (airway walls, vessels, connective tissue, etc.) alter the average ROI density according to their fractional volume. Cardiac motion, in which the beating heart may affect lung tissue density by motion of the adjacent parenchyma or beat-to-beat changes in pulmonary blood volume, also contributes to the background variability. The local magnitude of the cardiac effect depends on the proximity of the ROI to the ventricle. Gating of the image acquisition to the elecrocardiogram, which can be done with the ultrafast electron-beam CT scanners, greatly decreases this artifact (20). Finally, there is a small interscan variation in the CT number of a fixed target of 1-2 HU or less, attributable to measurement error within the scanner (13). Given these various sources of noise, the actual SSR for a typical 40-image wi/wo protocol in a normal animal averages around 8-12 HU², ranging from 3-5 HU² (for very smooth curves taken from the lung apex) to 40-80 HU² (for very noisy curves such as might be found in a ROI adjacent to heart). Occasionally, SSR >100 may be obtained in very small regions such as was seen in ROI 4 of Table 1, particularly if they are located such that there is a significant motion artifact. However, ultimately, it is this noise level, rather than an absolute limit imposed by ROI size, which determines the CI of the measurement.

Finally, Xe is a moderately soluble gas in blood or tissue, and uptake of Xe will alter the density of the non-gas-containing portions of a ROI. This phenomenon is the basis for the Xe-CT measurement of cerebral blood flow (4). Xe uptake by blood and lung tissue would result in changes in background density that would lag behind the gas-density changes and would be greater with longer exposure times and Xe concentrations. Examination of density changes of muscle or large blood vessels in the CT field of view suggests that the peak magnitude of this effect is 5-8 HU, when 60% Xe is used. This effect is not random in nature and thus is not included in the noise added by the MC method, but density changes due to tissue uptake could be modeled explicitly and their effect on the calculation of  determined through further simulation.

Conclusions. The use of imaging techniques such as Xe-CT for the noninvasive measurement of regional pulmonary ventilation requires the development of practical methods for the assessment of measurement CI and reliability. We have described the use of MC techniques for this purpose. The parameters and CI determined have the requisite properties of unbiasedness and coverage and they provide a conservative estimate of the CI as determined from repeated measurements. In addition, MC simulations were used to evaluate the relative merits of different imaging protocols over a wide range of realistic experimental conditions. The results of these simulations indicate that, of the three protocols examined, the wi/wo protocol provides the best estimate of the regional ventilatory , with the narrowest CI for a given set of imaging conditions. These approaches thus allow the study of methodological issues that would be difficult or time consuming to investigate experimentally.