Abstract
Simon, Brett A., Catherine Marcucci, Mansheung Fung, and Subhash R. Lele. Parameter estimation and confidence intervals for XeCT ventilation studies: a Monte Carlo approach.J. Appl. Physiol. 84(2): 709–716, 1998.—Xenonenhanced computed tomography (XeCT) is a technique for the noninvasive measurement of regional pulmonary ventilation from the washin and/or washout time constants of radiodense stable xenon gas, determined from serial computed tomography scans. Although the measurement itself is straightforward, there is a need for methods for the estimation of variability and confidence intervals so that the statistical significance of the information obtained may be evaluated, particularly since obtaining repeated measurements is often not practical. We present a Monte Carlo (MC) approach to determine the 95% confidence interval (CI) for any given measurement. This MC method was characterized in terms of its unbiasedness and coverage of the CI. In addition, 10 identical XeCT ventilation runs were performed in an anesthetized dog, and the time constant was determined for several regions of varying size in each run. The 95% CI, estimated from these repeated measurements as the mean ± 2 × SE, compared favorably with the CI obtained by the MC approach. Finally, a simulation was performed to compare the performance of three imaging protocols in estimating model parameters.
 xenonenhanced computed tomography
 lung imaging
imaging techniques are becoming increasingly important in physiology because of the need for noninvasive determination of regional structure and function in intact subjects. Xenonenhanced computed tomography (XeCT) 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 singleimaging 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 CTbased 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 XeCT 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 XeCT studies.
BACKGROUND AND MODEL
XeCT 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 curvefitting procedure, is equal to the inverse of the local ventilation per unit volume (specific ventilation) if singlecompartment 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, 10mm slice thickness, 2s scan, and 512 × 512 pixel image. It is not necessary that mechanical ventilation be used, but the same endexpiratory volume must be achieved for each image. Two or three baseline images while subjects are breathing room air or airO_{2}mixture are first acquired, and then the inspired gas is changed to the XeO_{2} mixture. Typically, 50% Xe40% O_{2} in air is used for animal studies to maximize the change in density; lower Xe concentrations may also be used. Sequential endexpiratory images are acquired until the lung is equilibrated with the XeO_{2} 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 crosssectional CT slice to an arbitrary region as small as 0.5 cm^{2} 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 singlecompartment exponential model by using a nonlinear least squares curvefitting 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_{0} is used to denote the baseline density value (without Xe), whereas D_{f} refers to the density after full equilibration with Xe.
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_{0} at every point to simulate the noise component of CTderived Xe tracer curves. Analysis of the residuals of actual fitted data by using normal probability plots (Fig.3) and the ShapiroWilks 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 ς^{2}[ε_{t} ∼N(0,ς^{2})]. For the washin model, this is represented as
ANALYTIC METHODS AND RESULTS
Parameter estimation. Parameters D_{0}, D_{f}, and τ are estimated from each set of data by using a nonlinear least squares curvefitting procedure (Igor Pro, Wavemetrics, Lake Oswego, OR). The starting points of the washin and/or washout segments (t_{0}, t_{1}) are initially estimated by visual inspection of the curves and then determined by the curvefitting procedure. Although data points occur once per breath, t_{0} and t_{1} 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 ς^{2} of the stochastic noise parameter ε_{t}was determined from the normalized summed squared residuals (SSR) of a given fitted data set
MC CI determination. CIs for a given set of XeCT 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_{0} 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 XeCT 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 (Table1). 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 longrun 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, τ,D̂_{0}, andD̂_{f} for 500 repetitions of simulated curves with added noise. The analysis was repeated for noise (variance ς^{2}) 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_{0} = 300, D_{f} = 400, and 40 points per curve. These density values are typical for wellexpanded 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.
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_{0} 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_{0} 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_{0}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_{0} 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_{0}, D_{f}} and noise level ς^{2}, a simulated curve was generated and curve fitted to provided parameter and noise (SSR) estimates (θ̂,ς̂^{2}). 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. Table3 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_{0} 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 XeCT 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 (120) and noise levels (ς^{2}4–100). The results are presented as onehalf CI width normalized to the parameter value, a measure of relative error (Fig.6). The onehalf 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.
As seen in Fig. 6, both D_{0} 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_{0}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
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 imagereconstruction 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 Xray 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 XeCT measurement of regional pulmonary ventilation. XeCT 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 repeatedmeasurements experiment.
Three different XeCT 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 firstorder (monoexponential) kinetics, and, as such, a deterministic model can be described for each protocol (Fig. 1). The assumption of firstorder behavior for regional lung washout implies that the ROI is behaving as a single wellmixed 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 singlecompartment 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 XeCT 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 XeCT 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 XeCT 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 XeCT 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 secondexponential 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 endexpiratory 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 beattobeat 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 electronbeam 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 40image wi/wo protocol in a normal animal averages around 8–12 HU^{2}, ranging from 3–5 HU^{2} (for very smooth curves taken from the lung apex) to 40–80 HU^{2} (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 nongascontaining portions of a ROI. This phenomenon is the basis for the XeCT 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 gasdensity 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 XeCT 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.
Acknowledgments
The authors thank Praxair Pharmaceuticals, Tarrytown, NY, for providing the Xe gas; Diversified Diagnostic Products, Houston, TX, for providing the Enhancer 3000 Xedelivery system; and Lifecare International, Westminster, CO, for providing the computercontrolled PLV102 ventilator.
Footnotes

Address for reprint requests: B. A. Simon, Dept. of Anesthesia, Tower 711, Johns Hopkins Hospital, Baltimore, MD 212878711 (Email:bsimon{at}welchlink.welch.jhu.edu).

This work was supported by the Johns Hopkins Department of Anesthesiology and Critical Care Medicine and by the American Heart Association Research Grant MDBG0495. B. A. Simon was supported by a Johns Hopkins School of Medicine Clinician Scientist award.
 Copyright © 1998 the American Physiological Society