Abstract
In internal radionuclide therapy, a growing interest in voxel-level estimates of tissue-absorbed dose has been driven by the desire to report radiobiologic quantities that account for the biologic consequences of both spatial and temporal nonuniformities in these dose estimates. This report presents an overview of 3-dimensional SPECT methods and requirements for internal dosimetry at both regional and voxel levels. Combined SPECT/CT image-based methods are emphasized, because the CT-derived anatomic information allows one to address multiple technical factors that affect SPECT quantification while facilitating the patient-specific voxel-level dosimetry calculation itself. SPECT imaging and reconstruction techniques for quantification in radionuclide therapy are not necessarily the same as those designed to optimize diagnostic imaging quality. The current overview is intended as an introduction to an upcoming series of MIRD pamphlets with detailed radionuclide-specific recommendations intended to provide best-practice SPECT quantification–based guidance for radionuclide dosimetry.
In radionuclide therapy, nonuniform activity distributions arise over a range of temporal and spatial dimensions from subcellular and cellular to organ, tumor, and whole-body levels (1). These activity nonuniformities lead to corresponding nonuniformities in tissue-absorbed dose. Traditionally, however, the dosimetric quantity that is reported is the mean absorbed dose to normal-organ or tumor volumes, and estimation of this quantity is based on the assumption of a uniform activity distribution within the source volumes. Furthermore, the calculation of the mean absorbed dose over the target volume is implicitly made assuming that this quantity is predictive of biologic effects. Although this may be true for stochastic biologic effects (e.g., cancer induction), there is mounting evidence that this is not the case for deterministic effects (e.g., tumor control and normal-organ toxicity). The implications of nonuniform dose and dose-rate distributions can be significant, as shown by mathematic modeling studies. For example, O’Donoghue demonstrated that a nonuniform tumor-absorbed dose distribution, characterized by “underdosing” certain parts of the tumor while “overkilling” other parts, might contribute to treatment failure by delivering sublethal doses to some clonogenic cells within the tumor (2). The spatial resolution (full width at half maximum [FWHM]) of current in vivo nuclear medicine imaging modalities (SPECT and PET) is on the order of 5–25 mm and thus can be used to determine the 3-dimensional (3D) spatial distribution of activity in volumes that are large relative to these dimensions. Serial quantitative images from these modalities can be coupled with methods for nonuniform dosimetry to obtain voxel-level absorbed dose rates. Corresponding dose–volume histograms can be coupled with computational models to yield radiobiologic quantities such as the equivalent uniform dose (EUD) and the biologic effective dose (BED), which, respectively, account for the biologic consequences of spatial and temporal nonuniformities in absorbed dose (2,3).
Many assumptions are implicit in internal dose calculations, and errors on the order of 30%–100% have been suggested (4). In the case of diagnostic radiopharmaceuticals, such errors are acceptable because of the large risk–benefit ratio and the fact that dosimetry is applied to patient populations to estimate stochastic risk and not to individual patients to plan treatment. However, in therapeutic applications, where toxicity and efficacy are of concern and where there is much less tolerance for inaccuracies, patient-specific dosimetry is justified. In particular, SPECT/CT-based 3D dosimetry may be essential in cases in which the average absorbed dose to an organ does not provide the information needed to predict or anticipate potential biologic effects. A regional or voxel-based calculation is thus required. Using a pretherapy tracer study, such calculations may be performed as a crucial step in optimized, treatment planning–based delivery of radiopharmaceutical therapy.
MIRD pamphlet no. 16 (5), a document on quantitative imaging techniques for dosimetry of internal emitters, focused on 2-dimensional planar imaging methods for determining the mean absorbed dose to target regions within the patient’s body. In contrast, the present MIRD pamphlet presents an overview of 3D SPECT requirements for dosimetry at both the regional level and the voxel level. This is the first in a series of reports on this topic. It provides general recommendations without specific examples. Subsequent reports in this series will provide specific recommendations and examples for individual commonly used radionuclides.
IMAGING-BASED DOSE ESTIMATION
The absorbed dose is defined as the mean energy imparted to target tissue per unit tissue mass. According to the MIRD schema, the mean absorbed dose, D(rT, TD), to target tissue, rT, over a dose-integration period, TD, from a radioactive material distributed uniformly within a source tissue, rS, is given by Equation 1 (6):
To determine the time-integrated activities in source regions or source voxels, serial measurements of activity must be made after administration of the radioactive material. Historically, the largest contributor to the error in absorbed dose estimation has been considered to be the inaccuracy in estimating time-integrated activity. The accuracy of dose calculations has improved with the use of the patient’s own images in Monte Carlo radiation transport calculations instead of generic anatomic models, such as reference man (Fig. 1) (10–12). As a result, accurate estimation of the time-integrated activity has become even more pressing.
Most dosimetric calculations in internal radionuclide therapy have relied on conjugate-view (anterior–posterior) planar imaging and geometric mean attenuation compensation for activity quantification, as described in MIRD pamphlet no. 16 (5). The planar method, however, cannot resolve the source depth or reliably correct for counts emanating from activity in tissue overlying or underlying structures of interest. The improved quantitative accuracy of 3D tomographic imaging modalities such as SPECT and PET over 2-dimensional planar imaging is well established (4,5,13). Presently, SPECT-based dosimetry is more widely applicable than PET, as most therapeutic radionuclides are not positron emitters and emit photons suitable only for SPECT (Table 1). In addition, advances in image reconstruction and quantification methods for SPECT have substantially improved both its accuracy and its precision. For 90Y, which does not have a γ-ray emission suitable for SPECT, a γ-emitting surrogate such as 111In is typically used. Recently, however, quantitative SPECT of 90Y bremsstrahlung has markedly improved and has enabled direct in vivo quantification of 90Y at therapeutic activity levels (14).
CURRENT CAPABILITIES FOR SPECT QUANTIFICATION AND VOXEL-LEVEL DOSIMETRY
Quantitative SPECT of therapeutic radionuclides can be more difficult than quantitative SPECT of diagnostic radionuclides such as 99mTc because of the higher-energy or multiple emissions typically associated with the former (Table 1). For example, the Monte Carlo–generated energy spectra in Figure 2 demonstrate the high fraction of scatter and penetration events in the SPECT acquisition window for typical therapy radionuclides. Nevertheless, with iterative reconstruction, confounding physical factors such as attenuation, scatter, and collimator–detector response (CDR) can be compensated for by modeling their effects in the projection and backprojection steps of the reconstruction algorithm (15,16). In addition, the recent clinical availability of hybrid SPECT/CT has provided the opportunity to further improve the accuracy of SPECT-based activity measurements for dosimetry. In particular, the CT anatomic image set facilitates compensation for nonuniform attenuation and for partial-volume effects. The CT image sets can also be used to define the anatomy and the density map for patient-specific dosimetry. The sequential SPECT and CT acquisitions during a single imaging session with a hybrid system eliminate much of the error and complexity associated with coregistration of SPECT and CT images acquired on different systems. Figure 3 shows patient images used in recent SPECT/CT-based internal dosimetry studies. The figure compares SPECT images reconstructed using the ordered-subsets expectation maximization (OS-EM) iterative algorithm without and with correction for image-degrading physical factors. The SPECT/CT acquisition and reconstruction methods corresponding to these images have been previously reported (11,14,17).
Table 2 summarizes the results of recent SPECT studies in which physical phantom measurements were used to evaluate the reliability of lesion or organ activity quantification of therapeutic radionuclides, as well as 99mTc. These studies used OS-EM reconstruction with attenuation correction, scatter correction, and in some cases CDR compensation and partial-volume correction (PVC). From Table 2, one can expect quantification accuracies of approximately 10% for most organs and lesions, but larger errors are to be expected for volumes with dimensions that are small relative to the SPECT image resolution. In vivo validation is generally more difficult than phantom validation. In a recent study, Zeintl et al. reported quantitative accuracy within 17% in 16 patients undergoing 99mTc-diphosphonate bone imaging by comparing SPECT/CT-based estimates of bladder activity with well counter measurements of urine activity immediately after imaging (18). Willowson et al. reported clinical quantification accuracies within 7% in lung ventilation–perfusion studies by comparing the SPECT/CT-based activity estimate with the known injected activity of 99mTc-macroaggregated albumin (19).
SPECT/CT-based dose calculation obviates model-based approximations and allows for nonuniform dosimetry down to the voxel level. The simplest approach for voxel-level dosimetry is to assume that the emitted energy is locally absorbed, that is, completely absorbed in the voxel in which the radiation was emitted. Such an assumption is valid only for particles whose maximum range in tissue is equal to or less than the voxel dimensions (typically >4 mm for SPECT). Local energy deposition has been assumed in previous voxel-level dosimetry calculations for 90Y β-particles (14) and for the β-component of the dose for 131I (12). Current computational approaches to voxel-level dosimetry that do not require an assumption of complete intravoxel energy deposition include the voxel S value method, the dose point–kernel method, and the Monte Carlo radiation transport methods (9). The voxel S value and dose point–kernel approaches are considered to be a reasonable compromise between simplified body or organ model–based calculations and more computer-intensive and time-consuming methods based on Monte Carlo radiation transport. Functional (e.g., from SPECT) and anatomic (e.g., from CT) imaging coupled with direct Monte Carlo radiation transport is generally considered to be the most accurate and most patient-specific of all currently available dose estimation methods (10).
QUANTITATIVE SPECT TECHNIQUES FOR NONUNIFORM DOSIMETRY
This section provides an overview of current best practices for quantitative SPECT to determine regional or voxel-level activities of therapeutic radionuclides or their surrogates.
Acquisition
In nearly all therapy imaging studies, parallel-hole collimation is used with selection of a low-, medium-, or high-energy collimator depending on the energy of the imaged photon and any significant higher-energy photons, as well as the desired balance of spatial resolution and sensitivity. For the radionuclides listed in Table 1, the suggested collimators are low-energy for 117mSn, 153Sm, and 186Re; medium-energy for 111In, 177Lu, and 67Cu; high-energy for 131I; and medium- or high-energy for 67Ga, 90Y, 166Ho, and 188Re. For some radionuclides that are imaged using a relatively low-energy photopeak, such as 166Ho and 188Re, a medium- or high-energy collimator is recommended because of septal penetration by higher-energy γ-emissions or by the considerable amount of bremsstrahlung caused by the β-emissions. The 9.5-mm (3/8-in)-thick NaI(Tl) crystal commonly found in current SPECT systems is best suited for photon energies up to 200 keV, but systems with thicker crystals should be used for imaging higher-energy photons such as the 364-keV γ-ray of 131I.
The projection-image matrix size (typically 642 or 1282) should be selected by considering the appropriate balance of spatial resolution and image noise. The ideal pixel size is considered to be smaller than half the spatial resolution (FWHM) of the SPECT system, measured at the center of rotation, for the isotope being imaged. The pixel size (with no zoom) can be calculated by dividing the digital field of view (FOV) of the camera by the number of pixels per row. For a SPECT camera with a FOV of 40 cm and an expected resolution of 20 mm in FWHM, a 642 matrix (pixel size, 6.25 mm) would satisfy the sampling requirement; if the expected resolution were 10 mm, however, a 1282 matrix (pixel size, 3.125 mm) would be required. SPECT data can be acquired using step-and-shoot or continuous gantry rotation. The latter is the more efficient and is the method of choice for a large number of projections. Body contouring is preferred over a circular orbit, as it minimizes the object-to-detector distance, thereby reducing the resolution-degrading effects of the distance-dependent collimator blurring. To minimize undersampling, the number of angular views over 360° should be at least equal to the projection-image matrix size (i.e., 128 views for a 1282 matrix). To minimize patient-motion artifacts, the total imaging time should be less than 30 min. When target contours can be obtained from a registered anatomic image, such as CT, shorter SPECT acquisition times (similar to those used for planar imaging) can give target activity estimates for which statistical noise is not the primary limiting factor in reliability (20). The acquisition energy window should be at least twice as wide as the energy resolution (FWHM) of the detector to avoid excessive count losses. Adjacent narrow windows below and above the main energy window are used if triple-energy-window scatter correction (21) is to be implemented (Fig. 2). For radionuclides that have multiple photon emissions, the energies of the emissions, as well as counting statistics, should be considered when one is selecting the photopeak for imaging. For example, for 111In, previous studies report similar accuracy using either one or both photopeaks, but with reduced noise with the latter (22,23). For bremsstrahlung imaging, the energy window setting is important even though there are no photopeaks in the energy spectrum (Fig. 2). Windows that include energies below 100 keV may prove challenging to quantify because of the large amount of down-scatter and the significant contribution of lead x-rays from the collimator. In a recent 90Y bremsstrahlung SPECT study, for example, a 105- to 195-keV window was used (14). Ideally, for each radionuclide, the energy window setting should be optimized (considering both quantitative accuracy and noise) using phantom experiments with well-calibrated activity or using Monte Carlo phantom simulations where the truth is known (24,25).
When SPECT/CT is performed, in most cases an unenhanced low-dose CT acquisition is sufficient for attenuation correction and is recommended to minimize the patient’s x-ray exposure. A low-dose CT acquisition is also sufficient for obtaining a reliable density map for Monte Carlo–based dosimetry. Low-dose CT typically results in a CT effective dose of less than 4 mSv (26). With integrated systems, the mechanical alignment of the 2 acquisition units is an important part of routine quality control. In addition, respiratory motion and patient movement between the sequential SPECT and CT scans must be minimized, and postacquisition evaluation of the degree of misregistration is needed. If the SPECT/CT misregistration is significant, the images need to be realigned either manually on the basis of visual comparison or using image registration software.
Dead-Time Correction
In radiation imaging systems, events occurring in close temporal proximity to a preceding event will be lost or malpositioned because of the short but finite time interval required to process each recorded event (system dead time) (27). In SPECT after therapy, where injected activities are often over 4 GBq, dead-time losses can be substantial and dead-time–related count losses must be corrected. Even with dead-time correction, however, such high counting rates may result in prohibitive image distortion. Dead-time correction is particularly important for radionuclides with multiple photon emissions such as 131I, as photons not included in the acquisition window also contribute to dead time. In a study on the dosimetric impact of dead-time correction after a 4-GBq therapeutic injection of 131I, correction for count losses led to an 11% increase in whole-body time-integrated activity (28). In another study, after a 3- to 4-GBq therapeutic injection of 131I, dead-time correction increased the measured activity at 1 d after therapy by up to 8% (29).
Before SPECT reconstruction of posttherapy data, the measured projection counts should be corrected for camera dead time. The simplest method for dead-time correction is based on monitoring counts corresponding to a reference source placed at the edge of the camera FOV (30,31). Here, the dead-time correction factor is determined as the ratio of counts recorded in the region of interest for the reference source without and with the patient present. The dead-time–corrected counts (true counts) for the patient study are then estimated by multiplying the measured counts by the correction factor. This method is problematic to implement, however, because of interference between photons emitted from the patient and those emitted by the reference source. Therefore, analytic dead-time correction methods based on mathematic models that characterize the system as paralyzable (every event can lead to dead time) or nonparalyzable (events during dead time are ignored) are preferable. These methods require preliminary calibration to experimentally determine the system dead time, which should be performed for the radionuclide of interest using the same window setting and scattering conditions as those encountered in patient studies (30–32). Once the dead time of the system and the measured counts are known, the dead-time–corrected counts can be determined from an analytic equation if the system is nonparalyzable or can be estimated by graphical or numeric methods if the system is paralyzable (27,30–32). The dead-time model that best approximates the behavior of a particular SPECT system should be determined empirically from additional measurements of observed counting rate over a range of different activity levels. In a graph of observed counting rate versus activity, the counting rate increases asymptotically toward a maximum in the nonparalyzable case, whereas it rises to a maximum and then decreases in the paralyzable case (27). In recent high-counting-rate quantitative studies involving SPECT/CT systems, a combined paralyzable–nonparalyzable model was used for 99mTc dead-time correction (19) and a paralyzable model was used for 131I dead-time correction (33).
Image Reconstruction
Although filtered backprojection is still widely used in clinical practice for SPECT reconstruction, iterative reconstruction should be used for quantitative studies. Iterative reconstruction allows for optimal correction for image-degrading physical effects and for improvement in noise properties (15,16). In iterative reconstruction, the set of algebraic equations that model the SPECT acquisition process is solved by successive refinement. Iterative reconstruction consists of a criterion that is optimized and an iterative algorithm to find the activity distribution that optimizes that criterion. Because the projection data are corrupted by statistical noise, statistical criteria are preferred. The most widely used criterion is the maximum-likelihood (ML) criterion, which seeks the activity distribution that maximizes the Poisson likelihood of the projection data. The first algorithm proposed to seek the ML criterion in emission tomography reconstruction was the expectation maximization (EM) algorithm (34). One limitation of the ML-EM algorithm is its slow convergence. An alternative is the OS-EM algorithm (35), now available with most commercial SPECT systems. The OS-EM algorithm updates the estimate multiple times per iteration using a different subset of the projections in each update. Although OS-EM is practically useful, it is not guaranteed to converge and can give suboptimal results for noisy data. The acceleration achieved with OS-EM, compared with ML-EM, is usually on the order of the number of subsets. A reasonable compromise between speed and reconstruction quality is to have at least 4 projections per subset (15), though for noisy data (e.g., fewer than 50,000 total counts per slice), more projections per subset should be used.
Iterations
Important considerations when using iterative reconstruction are the number of iterations and the use of postreconstruction filters to suppress noise. The number of iterations should be determined by considering the bias–variance trade-off, the effects compensated for during the reconstruction, the dosimetric quantity of interest, and practical considerations such as the computation time. The bias, which is the difference between the mean reconstructed image and the true activity distribution, is reduced as iterations proceed. However, with more iterations, the variance (noise level) of the image increases. When CDR compensation is included in the reconstruction to provide resolution recovery, edge artifacts are observed, and these also become more prominent with additional iterations (36,37). These artifacts are particularly problematic in 3D dosimetry applications because of the resulting distortion of the activity distribution. To determine the number of iterations to use in patient imaging, experimental studies with a physical phantom or Monte Carlo simulations can be performed to evaluate bias and noise as a function of number of iterations for the radionuclide and SPECT system of interest. In such studies, the true value of the quantity to be estimated (e.g., the total source-region activity or the activity distribution) is, of course, known. For example, Figure 4 shows the 131I source-region activity recovery as a function of the number of 3D OS-EM iterations for a physical phantom experiment with multiple hot spheres in a warm background performed on a SPECT/CT system. (Recovery is defined as the ratio of activity estimated from the image to the true activity in the object). As shown in Figure 4, high reconstruction accuracy (recovery > 90%) is achieved for larger objects after about 30 iterations, but for smaller objects, significantly more iterations are needed to improve the reconstruction accuracy. In general, more iterations are required for convergence when compensation for image-degrading factors is performed, as compared with iterative reconstruction with no compensation. When the quantity of interest is the mean absorbed dose to a target region, a larger number of iterations can be used because averaging over the target region will tend to reduce the effects of noise and edge artifacts (both of which increase with the iteration number). However, when the quantity of interest is 3D dose metrics such as the dose–volume histogram, or quantities estimated from it, a smaller number of iterations may be optimal.
Filtering
Postreconstruction low-pass filtering is often applied to images to reduce noise and edge artifacts, improving the perceived image quality. However, this filtering can degrade image resolution and thus increase partial-volume effects. Consequently, postreconstruction filtering is not desirable for quantifying total target activity. On the other hand, when the images are used to estimate activity distributions to calculate 3D dose metrics such as the dose–volume histogram, postfiltering is often desirable to suppress noise effects (with CDR compensation the need for filtering may be less, but some level is desirable). As an example, Figure 5 shows cumulative dose–volume histograms of the liver for OS-EM reconstructions with and without postreconstruction filtering (with projection data generated by phantom simulations). Also shown are the cumulative dose–volume histograms for the true activity distribution (phantom). Noise and partial-volume effects can explain the deviation of the histogram obtained with OS-EM from the true histogram: noise results in voxels having estimated activities (and thus dose rates) higher or lower than the true values, and partial-volume effects reduce the estimated activities (and thus dose rates) at the edges of objects with higher activities than their surroundings. Postreconstruction low-pass filtering improved the histogram, despite the fact that it tends to increase partial-volume effects. This is because the filtering tends to reduce the effects of noise on the histogram in the interior of objects with near-constant activities.
Compensation for Image-Degrading Effects: Attenuation, Scatter, and Detector Response
Iterative reconstruction–based compensation for attenuation is performed by including attenuation factors in the system–transition matrix of the system model (38). Similarly, compensation for the CDR is performed by incorporating the depth-dependent CDR in the system matrix (39). Because the distribution of scattered photons is a complex function of the emitter and the absorber, incorporating the scatter estimate directly in the system matrix leads to a considerably less sparse matrix than if only attenuation factors and CDR are included. To avoid the long computations associated with a nonsparse system matrix, a constant precalculated scatter estimate could be included in the statistical model of the iterative reconstruction by adding it to the forward-projection estimate (40). This approach of including the scatter estimate in the iterative algorithm itself yields better noise properties than if the scatter estimate is subtracted from the projections before image reconstruction (41). Methods for obtaining the attenuation map, scatter estimate, and CDR are described below.
Estimation of the Attenuation Map
Because most anatomic regions in the body are heterogeneous in tissue composition, determination of an accurate and patient-specific attenuation map for nonuniform attenuation compensation is critical. The attenuation map is a voxel-by-voxel representation of the linear attenuation coefficients (μ) at the SPECT photon energy and must be well registered to the SPECT image. Such maps can be generated by transmission scanning with a radionuclide source, but CT-based attenuation maps are preferred as they have lower noise, better spatial resolution, and better contrast (38) and are generally faster and easier to acquire.
For generation of the CT-based attenuation map, CT images that are reconstructed in a larger image matrix (typically 512 × 512) with a smaller transaxial voxel size and slice thickness than the SPECT image (typically 64 × 64 or 128 × 128), must first be down-sampled to the SPECT format and then spatially aligned with the SPECT image set. The coregistration and down-sampling of SPECT and CT are facilitated by hybrid imaging, in which the patient is not repositioned between the SPECT and CT scans and the software included with the scanner generates appropriately registered and down-sampled attenuation maps. The registration can be problematic when CT images acquired on a separate system are used to generate the SPECT attenuation map. After coregistration and down-sampling, the CT voxel values, expressed in Hounsfield units, must be transformed to the μ values of the corresponding tissue. Although transformation methods based on segmentation can separate the CT image into regions that represent different tissue types (with known values of μ), the commonly used approach is based on voxel-by-voxel scaling and thus does not require image segmentation (42). Because linear attenuation coefficients are energy-dependent, the CT values at the x-ray energy (typically a polychromatic beam having an effective energy of 50–80 keV) must be scaled to the energy of the SPECT photons. For low-Z materials such as air, water, and soft tissue, a single scaling factor can be used to convert from x-ray energies to SPECT photon energies. For spongiosa and cortical bone, however, simple linear scaling is not adequate. Bilinear scaling has been proposed, noting that CT numbers in the range of −1,000 to 0 Hounsfield units correspond to regions that contain mixtures of lung (largely air) and soft tissue whereas CT numbers with greater than 0 Hounsfield units correspond to regions that contain mixtures of soft tissue and bone (42,43). To determine the appropriate scaling factors to convert Hounsfield units to μ values, a system-specific CT calibration measurement should be performed with materials of known composition. For example, Figure 6 shows the relationship between CT number and μ for 16 materials in a calibration phantom imaged on an integrated SPECT/CT system. For hybrid SPECT/CT systems, such calibrations are typically performed by the manufacturer.
Estimation of the Scatter Component
Photons that undergo Compton scatter in the patient, the camera, or the surrounding material degrade the quality of the reconstructed image and the accuracy of quantitative analysis. Scatter is particularly problematic for imaging of bremsstrahlung photons and for radionuclides that have γ-ray emissions with energies greater than the photopeak energy of interest, as these photons can be scattered and result in counts within the acquisition window, a process often referred to as down-scatter (Fig. 2). Several methods have been proposed for SPECT scatter correction (40); these can be categorized as energy distribution–based or spatial distribution–based. However, not all these methods are practical in the clinical environment and only a few, such as the dual-energy or triple-energy-window (TEW)–based estimation methods and the effective scatter source estimation (ESSE) method, have been implemented commercially. For therapy radionuclides, the TEW scatter correction is often most appropriate because of its ease of implementation, especially when down-scatter is an issue. The more sophisticated ESSE or similar model- or Monte Carlo–based corrections are recommended if multiwindow acquisition is not available or if window-based correction is not suitable, as in the case of bremsstrahlung imaging.
In energy distribution–based methods, the number of scattered photons in the photopeak energy window is estimated on the basis of counts in one or more scatter windows. In the TEW method, the scatter is estimated as the area of the trapezoid beneath the line joining the 2 adjacent narrow scatter windows (for radionuclides with a single γ-ray emission, the upper scatter window can be omitted). For each projection pixel i, the number of scatter counts in the photopeak window, Ci,scat, is calculated as:
Because of the use of narrow windows, the TEW estimate is noisy, and low-pass filtering should be used to reduce noise. Excessive filtering distorts the scatter estimate, whereas minimal filtering accentuates noise (44). Generally, wider scatter windows can reduce the noise but result in poorer accuracy of the trapezoidal approximation of scatter. Ideally, the scatter window settings should be optimized considering both accuracy and noise, preferably by Monte Carlo simulation studies since they allow the separation of scatter events (Fig. 2) and thereby allow comparison of the true and TEW-estimated scatter components. For example, in a 111In Monte Carlo study investigating multiple window settings, good agreement between the TEW-estimated scatter-to-total count ratio and the true scatter-to-total count ratio was demonstrated using 6% lower and upper scatter windows adjacent to the 171- and 245-keV photopeaks (25). In a 131I study, good quantitative accuracy was demonstrated using a 20% photopeak window at 364 keV and adjacent 6% upper and lower scatter windows (45). These studies did not evaluate counting rate effects; the narrow windows make the TEW correction sensitive to changes in the energy spectrum because of, for example, pulse pile-up during high-counting-rate posttherapy imaging (33).
ESSE (46) is an example of a spatial-domain scatter estimation method in which modeling of scatter in the projection data is based on the estimate of activity distribution. The advantage of model-based methods over energy-window–based methods is that they do not require additional acquisition windows or optimization of parameters such as the window width and filtering. In addition, window-based methods are not suited for bremsstrahlung SPECT because of the continuous energy distribution of primary (unscattered) bremsstrahlung photons (Fig. 2). ESSE works by calculating an effective scatter source. The attenuated projection of the effective scatter source then gives the scatter component of the SPECT projections. The major assumption in ESSE is that the photon propagates through water from the point of emission to the last scatter point. Consequently, inaccuracies will be largest for regions where the attenuation distribution is nonuniform and cases in which down-scatter is important. However, for many practical isotopes of interest (such as 99mTc, 123I, 111In, and 131I), ESSE has been shown to provide good estimates of scatter. A shortcoming of ESSE and of any spatially based scatter estimation method is the difficulty of compensating for scatter from activity outside the FOV, resulting in reduced accuracy at the edge of the FOV.
Other, more sophisticated, approaches to scatter modeling, which are expected to provide better quantitative accuracy, include Monte Carlo simulation–based methods (45,47,48) and analytic calculations using the Klein–Nishina formula for Compton scattering (49). These methods are computationally intensive and have been limited to use in research settings but are expected soon to enter the clinical arena as computing power continues to improve.
Modeling the Camera Response
The CDR of the system, which refers to the image generated by a point source, has 4 components: the intrinsic response of the detector and the 3 components of the collimator response (geometric, septal penetration, and scatter). For scintillation cameras, the intrinsic response is well modeled as a Gaussian function corresponding to a spatial resolution of 4 mm or less in FWHM for modern systems. The geometric response of the collimator can be calculated theoretically and is often approximated as a Gaussian function with a FWHM that is a linear function of source-to-collimator distance. In the absence of collimator scatter and penetration, which is a reasonable approximation for low-energy photon emitters, the CDR can be obtained by convolving the intrinsic response of the detector with the theoretic geometric response of the collimator (39). Such a gaussian CDR model is often used in commercial iterative reconstruction software, regardless of the energy of the photons being imaged. A gaussian CDR model, however, does not account for the septal penetration and scatter components of the collimator response function. The Monte Carlo–simulated images of Figure 7 demonstrate the significance of these components for imaging of therapy radionuclides that emit medium- to high-energy photons or that produce bremsstrahlung. Unlike the geometric component, the scatter and penetration components of the CDR function cannot be determined theoretically but can be determined by Monte Carlo simulations (Fig. 7). When full CDR compensation is desired, the total (i.e., including intrinsic, geometric, penetration, and scatter components) distance-dependent CDR can be determined experimentally by measuring the in-air point-source response at various distances from the collimator face (50,51) or by Monte Carlo simulations of this geometry (52). The total CDR should be determined for the radionuclide of interest, as both septal penetration and scatter are energy-dependent (in contrast to the geometric component). For experimental measurement, a pointlike source can be prepared by saturating a small (∼5-mm diameter) disk of filter paper or an ion-exchange-resin bead with a radioactive solution and sealing it with tape. The source can then be placed on a low-scatter medium (such as a Styrofoam [The Dow Chemical Co.] cup) positioned at the required distance from the collimator. The planar acquisitions to determine CDR should be performed using a fine image matrix (512 × 512) and for at least 5 source–collimator distances in a clinically relevant range (2–30 cm) (50,51). For medium- and high-energy photon emitters, a gaussian function in combination with decreasing exponential terms can be used to fit the measured CDR (50,51). The exponential terms are used to account for the wider tails of the response function due to collimator scatter and penetration.
Definition of Target
Because 3D dosimetry provides an absorbed dose estimate for each image voxel, definition of the target (organ, tumor) for estimation of the target mass is not always critical. However, after the 3D (i.e., voxelwise) dose calculation, the quantity of interest could be dose–volume histogram statistics or EUD over the target volume, which do require definition of the target volume. In previous SPECT-based 3D dosimetry studies, both high-resolution anatomic imaging (11,12) and SPECT thresholding (53) have been used to define the target volume of interest.
Although PET is emerging as an important modality for the delineation of target volumes (54), SPECT-based image segmentation is difficult because of the relatively coarser spatial resolution and higher noise content of SPECT images. Consequently, the preferred method is to define the target (manually or using automatic threshold-based contouring) on a high-resolution anatomic image such as CT. The contour is then copied to the SPECT image after the low-resolution SPECT image has been interpolated to the matrix size of the reference CT image and has been spatially aligned with it. This process is facilitated by hybrid SPECT/CT; however, CT acquisitions on such systems are typically not of diagnostic quality (performed in low-dose mode and without contrast) and can have poor soft-tissue contrast, which can affect the accuracy of target contouring. In addition, volumes defined on anatomic imaging modalities may not actually be the same as the volume within which the radionuclide localizes because of differences in physiology and anatomy within the volume of interest. When anatomic imaging–based target delineation is not feasible, for SPECT-based target delineation either a fixed threshold or more sophisticated adaptive thresholding methods based on the gray-level histogram of pixels within the volume of interest can be used (55). When fixed thresholds are used, the threshold required for accurate target delineation is a function of object shape and size relative to the system spatial resolution and target-to-background image contrast (56). The optimum threshold for patient studies should therefore be based on system-specific phantom calibration experiments (to determine threshold-volume curves) and a priori estimates of the volume and contrast of the anatomic structure. In 2 such optimized studies using threshold-volume calibration curves for spheres (acting as surrogates for lesions), the average difference between the lesion volumes estimated from thresholding the emission image and those defined on the CT image were less than 15% (57,58). The threshold that gives the best estimate of the object volume is not necessarily the same as the threshold that gives the best estimate of the object activity.
Camera Calibration Factor for Absolute Quantification
In nuclear medicine diagnostic imaging, only relative quantification is typically needed. However, absolute SPECT quantification is a requirement for imaging-based internal radionuclide dosimetry. After reconstruction, the system sensitivity or calibration factor must be known to convert the reconstructed SPECT voxel values to activity or activity concentration. The most reliable method of determining the calibration factor is to perform an experimental measurement with a known amount of activity or activity concentration of the isotope of interest. This calibration measurement can be done in 2 ways, the simplest of which is a planar acquisition of a pointlike source (or a small volume of activity in a syringe or Petri dish) to determine the in-air camera sensitivity (19,49,52,59). Such a calibration, though practical to implement, is appropriate only if the patient SPECT data are reconstructed with high accuracy, including perfect correction for scatter and attenuation. Because this standard is difficult to achieve, the more reliable calibration measurement is a SPECT acquisition with a source geometry that better approximates the scatter and attenuation properties in patient imaging, such as a tank of uniform activity or hot spheres in uniform background activity, thereby reducing the effects of imperfect compensation (18,50). Acquisition, reconstruction (including correction for attenuation, scatter, and CDR), and target definition in the calibration study should be performed in the same manner as in the patient study. The calibration factor (in units of counts/s/MBq) is then determined by dividing the total reconstructed counts within the target volume of interest by the product of the known activity and the SPECT acquisition time.
Because camera sensitivity may vary with time, the calibration experiment needs to be repeated periodically. However, repeating phantom experiments may not be practical in a clinical environment, especially with longer-lived radionuclides. In this case, an approach combining SPECT data from a single phantom measurement with periodic in-air point-source planar measurements can be used. The calibration factor from the initial phantom study is updated (i.e., scaled) on the basis of the ratio of the measured point-source sensitivity performed at the time of patient imaging to that measured at the same time as the measurement of the phantom-based calibration factor.
PVC
The partial-volume effect affects activity determination mostly in small structures (i.e., having dimensions less than 3 × FWHM of the spatial resolution of the reconstructed images), but it can also affect activity at the edge of larger objects because of count spill-in and spill-out (60). In iterative reconstruction, partial-volume effects can be reduced by including the distance-dependent CDR function in the model of the image formation process, but this does not result in complete resolution recovery. Therefore, PVC should be used for improved activity quantification even when CDR compensation is used. Experimental and simulation studies of 131I and 111In SPECT have shown the importance of PVC after 3D OS-EM reconstruction for accurate quantification of activity in small structures or structures close to regions with high activity (61,62).
Most PVC methods either are empiric methods based on physical phantom measurements or are anatomy-based postreconstruction methods that use coregistered high-resolution CT or MR images. In the empiric methods, physical phantom measurements with simple geometric shapes such as spheres of varying sizes are used to determine PVC factors or recovery coefficients (RCs) for similarly sized anatomic structures (Fig. 4). Such phantom measurements must be performed on the same system as used in clinical imaging, and data must be acquired and reconstructed in the same manner as the patient data. This empiric recovery correction is easy to implement and is most applicable to anatomic structures that can be approximated by simple geometric shapes such as spheres. For example, a recent experimental study demonstrated that the RCs for ellipsoid lesions with major-to-minor axis ratios of less than 2 are accurately approximated (to within 10%) by the RCs for volume-equivalent spheres (63). Instead of physical phantom measurements, Monte Carlo simulations with accurate characterization of the imaging system can also be used to determine RCs. In a recent study, RCs determined as a function of simulated sphere size, position, and reconstruction parameters were used when quantifying activity in anatomic structures of patients (18).
Volume-based PVC using RCs can be insufficient because RCs depend on a variety of factors other than size, including shape, object-to-background contrast, and position in the image (because of spatially varying resolution). In addition, although RCs can be used to improve activity quantification at the regional level, they cannot provide corrections on a voxel-by-voxel basis. More robust approaches to PVC than the empiric methods use coregistered anatomic information. Anatomy-based PVC is highly susceptible to misregistration and segmentation errors; hence, these methods are best suited for hybrid SPECT/CT. Anatomy-based PVC can be performed for the object as a whole (with the assumption of uniform uptake within that volume) and on a voxel-by-voxel basis to represent the true activity distribution. One common method at the object level is the geometric transfer matrix method, which was originally developed for brain PET and has been shown to significantly improve SPECT quantitative accuracy (64). The perturbation-based geometric transfer matrix method, which is a refinement to account for possible nonlinear effects inherent in iterative reconstruction, has been used for SPECT PVC of organ activities in imaging-based dosimetry applications (52,62). Voxel-by-voxel PVC, which has been applied mostly to PET, is significantly more complex than object-level corrections (65). For SPECT, there have been some attempts at voxel-by-voxel PVC (66,67), but there is currently no widely accepted, well-validated method.
Time-Integrated Activity
For imaging-based dosimetry, the time-integrated activity in source regions or voxels must be determined from serial quantitative imaging. The temporal distribution and number of imaging samples are important for defining the values of the measured transit time distribution. The number of samples required is related to the number of compartments, or separate exponential terms, in the source region time–activity function, whereas their optimal timing is related to the clearance rate for each exponential term. Rapidly clearing activity, for example, will necessitate frequent and early sampling, whereas slowly clearing activity requires fewer and more widely dispersed measurements. The area under the time–activity curve represents the time-integrated activity and can be determined by numeric integration or by analytic methods in which the measured data are fit to a summation of exponential terms or another mathematic function that can be analytically integrated. An alternative approach for determining time-integrated activity is based on compartmental modeling. Details of pharmacokinetic modeling, including sampling requirements, are beyond the scope of this report, and the reader is referred to other publications (5,68,69). Determination of time-integrated activity at the voxel level introduces additional challenges due to image noise, noise propagation, and the need for accurate coregistration of the serial images (70).
Organ- or tumor-level time-integrated activities can be determined with or without registration of serial SPECT images. If the image series can be accurately registered to a reference scan, a single volume of interest defined on the reference scan can be applied to all the images in the time series. This approach, however, is not recommended for defining volumes that change significantly among the imaging sessions; malignant lymphomas, for example, can shrink dramatically within days of treatment. In this case or if the serial images cannot be accurately registered, organ and tumor volumes should be redefined on each of the serial scans to generate the appropriate time–activity curves. This step avoids the need to register the serial images, but the variability in target volume definition on each scan introduces an additional source of error. A study evaluating the effects of misregistration and inaccuracy of delineation of organ volumes on SPECT activity quantification concluded that, generally, larger errors were associated with misregistration (71).
Voxel-level determination of time-integrated activity requires registration of the SPECT images to achieve spatial concordance among the images in the time series. The registered SPECT image sets can be integrated, voxel by voxel, over time to obtain the 3D time-integrated activity image. However, SPECT–SPECT image registration is challenging because of the coarse spatial resolution, lack of anatomic information, noise, and temporal changes in the activity distribution. Manual rigid alignment, achieved by rotation and translation of 2 SPECT image sets based on visual comparison, is generally insufficient for voxel-level accuracy. Consequently, automated and semiautomated techniques that rely on external fiducial markers attached to the body surface, as well as intrinsic methods that rely only on the patient-generated image content, have been used to register sequential SPECT images for 3D dosimetry (72–74). Registration using external markers overcomes some of the resolution and noise limitations of the SPECT data and does not require complex optimization algorithms, assuming that internal structures maintain the same respective locations relative to the external markers. For intrinsic methods, the registration is based on aligning anatomic landmarks or segmented structures or on optimizing voxel similarity metrics of the multiple datasets. Voxel-based registration using similarity metrics operates directly on the image gray-scale values in the entire volume of the multiple image sets, in contrast to other registration methods, which are based on identifying equivalent features in the image sets. Normalized mutual information has been used as a similarity metric for both rigid and nonrigid registration of serial SPECT/CT scans in which the registration of the SPECT images was accomplished by first performing CT–CT registration and applying the resulting transformation parameters to the corresponding SPECT images (74). When accurately registered SPECT/CT data are available, as given by hybrid imaging devices, CT-based registration of serial SPECT images is preferred because it avoids problems associated with the coarse SPECT spatial resolution and noise, as well as temporal changes in the activity distribution. Even with hybrid imaging, measures to minimize spatial mismatch due to mechanical misalignment and patient motion must be considered, and registration of the SPECT/CT images should be visually verified.
Hybrid Planar/SPECT Methods
Three-dimensional pharmacokinetic input is required for 3D dosimetry calculations. To date, most SPECT-based dosimetry calculations have used multiple planar studies to obtain pharmacokinetic data and a single SPECT scan to scale the pharmacokinetic curve for the region of interest (12,17,75,76). This approach enables 3D dosimetry analysis with only a single SPECT scan. The 3D pharmacokinetic input is derived from the single SPECT scan by assuming that the spatial activity distribution measured at the single time point remains fixed and that the kinetics for the region of interest in the planar image applies to the 3D activity distribution within the corresponding volume of interest. However, serial PET using 124I in thyroid cancer patients has shown that the activity distribution of radionuclides within tumors is not fixed over time and that adoption of this approach will lead to absorbed dose distributions that do not reflect the activity distribution obtained by a 3D representation of the kinetics (i.e., by multiple PET or SPECT scans) (70). Here, the kinetics in 3 high-activity foci within a 5-cm3 tumor were compared with overall tumor kinetics, and a 2-fold greater uptake rate was found in 2 foci than in the entire tumor. Depending on the time point chosen to measure the 3D activity distribution, a hybrid method that applies a single time–activity curve for the entire tumor would have led to a 56% underestimate of the absorbed dose to one of the foci and a 10% overestimate to another focus. Similar concerns might be anticipated for normal-organ dosimetry in some instances (e.g., for kidneys in radiopeptide therapy). In general, the error introduced by assuming a fixed relative spatial activity distribution is exponentially magnified if the resulting absorbed dose (or BED) distribution is used to calculate the EUD for tumors or to assess the probability of normal-tissue complications. Analysis of quantitative SPECT acquisition protocols has shown that the time required for SPECT can be reduced substantially so that it is competitive with the time required for planar imaging if the objective is activity quantification rather than clinical diagnosis (20). In light of these considerations, acquisition of serial SPECT studies is recommended to obtain the pharmacokinetic input for fully 3D dosimetry calculations. If this is not feasible, a hybrid planar/SPECT approach may be adopted as long as its limitations are understood.
CONVERSION OF QUANTITATIVE SPECT IMAGES TO ABSORBED-DOSE OR BED IMAGES
After quantitative SPECT the main steps for 3D dosimetry are as follows:
Register serial SPECT or SPECT/CT images to each other as outlined previously.
Calculate a dose-rate image for each of the SPECT time points to obtain a dose-rate map, using one of several methods for dose-rate estimation at the voxel level: assuming complete local (i.e., intravoxel) absorption of short-range electrons, voxel S values, dose-point kernel convolution, or Monte Carlo radiation transport. For Monte Carlo–based calculation, the corresponding registered CT image at each time point of step 1 can be used to account for density and composition variations in the dose-rate calculations.
Integrate the dose-rate images over time, either analytically, by fitting a mathematic expression, or numerically, or by a combination of the 2 approaches to obtain absorbed dose images. Alternatively, the BED distribution may be calculated directly from the time dependence of the dose rate distribution.
Present the 3D absorbed dose (or BED) distribution as parametric images, isodose contour plots, or dose (or BED)–volume histograms for tumors and normal organs.
The results of step 4 can be used as input for radiobiologic models to calculate EUD or normal-tissue complication probabilities over segmented structures of interest.
The calculation of BED (3) and EUD (2) is beyond the scope of the current report, but the reader is referred to the references for details. Normal-tissue complication probability modeling using either 3D dose or BED distributions has been applied to external-beam radiotherapy (77–79) and to targeted therapies of the kidney (80). The BED distributions have been corrected to represent a nominal dose per fraction for external-beam radiation or to an extrapolated low dose rate for protracted radiation.
CONCLUDING REMARKS AND FUTURE DIRECTIONS
Quantitative SPECT can be accurate to better than 10% with state-of-the-art SPECT/CT hybrid systems, iterative reconstruction algorithms that incorporate corrections for image-degrading physical factors, and careful system calibration. Such accuracy is possible when quantifying activity in most organs and moderate-sized tumors and even with radionuclides with complex decay schemes (Table 2). However, activity quantification in structures that are small relative to the SPECT resolution (currently 5–25 mm) remains problematic even if CDR compensation is used. When SPECT images are used for 3D dosimetry, the degradation in the images due to the finite spatial resolution is an important consideration. For small volumes, the SPECT-derived activity appears to be distributed over a larger volume, resulting in underestimation of activity on a voxel level. In this case, if explicit 3D methods for dosimetry are applied using SPECT images as the input, the absorbed doses per voxel will be underestimated. Future developments in high-resolution imaging, including voxel-level partial-volume compensation and improved noise suppression methods, should lead to more accurate quantitative imaging of the heterogeneous activity distributions within various tissue structures.
When imaging-based dosimetry results are reported, it is important to include details of the SPECT activity quantification methodology such as acquisition parameters, reconstruction and compensation methods and parameters, and the calibration method, as specified in a recent guidance document issued by the European Association of Nuclear Medicine (81). Before application in patient dosimetry studies, the SPECT quantification methods must be carefully validated by evaluation in phantom studies, as in the examples in Table 2, and, to the extent possible, by in vivo studies. For phantom-based validation, the experimental conditions, such as the imaging geometry, activity distribution, and counting rates, must closely approximate those encountered clinically. Monte Carlo simulation studies can also be used for such evaluation, but the Monte Carlo algorithm must first be validated for the specific SPECT system. Clinically realistic Monte Carlo simulations are possible with voxel-based anthropomorphic phantoms (82) and with dynamic mathematic phantoms that allow accurate modeling of anatomic variation and patient motion (83).
As the use of rapidly clearing therapeutic radiopharmaceuticals increases (e.g., peptides, engineered antibody fragments, and small molecules), it is clear that additional attention will focus on the collection of quantitative 4-dimensional dynamic distributions during critical phases of the biokinetic processes, with rigorous validation of methodologies used to obtain such data. The methods described here will need to be updated and validated using in vivo and ex vivo measurements. New technologies that promise to facilitate fast 3D (i.e., 4-dimensional) acquisition and processing are in varying stages of evaluation. Examples include the use of rapidly rotating γ-cameras with slip-ring gantries to continuously acquire list-mode data that can be reconstructed into 4-dimensional datasets amenable to kinetic analysis. Slowly rotating standard γ-cameras have also been used to reconstruct dynamic time imaging sequences using specialized software (84). In addition, multiple stationary γ-cameras with multipinhole collimation have been used successfully to acquire 4-dimensional images of animals and humans. These are exemplified by several generations of fast nonrotating multidetector SPECT research systems (85), and such systems are now commercially available for use in both animal and patient studies (86–88). Fast electronics for digital signal processing have greatly increased the counting rate capability of γ-cameras to 2–4 million counts per second (89).
Even if the corrections for all physical effects in PET and SPECT were accurate and complete, finite spatial resolution and large voxels will continue to introduce inaccuracies, especially at the microscopic level, given that radiopharmaceuticals are often distributed heterogeneously within a voxel. Therefore, in most cases radionuclide absorbed doses are calculated using an inappropriately large target-region mass. In other words, the energy deposited from emitted particles is not uniformly deposited, and even voxel-level absorbed doses may not reflect the actual biologic effect. This problem potentially can be addressed using subvoxel models of tissue. Such an approach has been applied to bone marrow, for which calculation of the absorbed dose incorporates a macroscopic-to-microscopic conversion (90). Recently, in situations in which SPECT and CT resolutions are inadequate for explicit voxel-based dose calculations, a hybrid voxel and idealized (Monte-Carlo simulation–based) geometry calculation has been implemented to estimate the absorbed dose to normal brain from brain lesions in a thyroid cancer patient (91) and to the vascular wall in lymphoma patients whose tumor encircles aortic vessels (92). This approach assumes a uniform activity distribution in structures below the imaging resolution of SPECT or CT and provides dose–volume histograms in relevant target regions. Further, ancillary procedures such as digital autoradiography and α- and β-camera imaging of excised tissue samples may provide small-scale distribution data (93,94). Although this approach may be limited for tumors (in terms of translation of results among tumors) because of intertumor variability, normal organs such as liver, bone marrow, and kidney may be amenable to this approach for providing more biologically relevant absorbed doses at the subvoxel level than given under the assumption of uniform energy deposition at the whole-voxel level.
Acknowledgments
This work was supported in part by grants R01 EB001994, R01CA109234, and R01CA116477 awarded by the National Institute of Health, U.S. Department of Health and Human Services, and a grant from the Swedish Cancer Foundation. Code using the ESSE scatter compensation and some of the other methods discussed in this paper has been licensed to GE Healthcare for inclusion in a commercial product. Under separate licensing agreements between GE Healthcare and the Johns Hopkins University and the University of North Carolina at Chapel Hill, Dr. Frey is entitled to a share of royalty received by the universities on sales of products described in this article. The terms of this arrangement are being managed by the Johns Hopkins University in accordance with its conflict-of-interest policies. No other potential conflict of interest relevant to this article was reported.
Footnotes
In collaboration with the SNM MIRD Committee: Wesley E. Bolch, Darrell R. Fisher, Roger W. Howell, Ruby F. Meredith, and Barry W. Wessels
Published online Jun. 28, 2012.
- © 2012 by the Society of Nuclear Medicine and Molecular Imaging, Inc.
REFERENCES
- 1.↵
- 2.↵
- 3.↵
- 4.↵
- 5.↵
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.↵
- 14.↵
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.↵
- 41.↵
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.↵
- 61.↵
- 62.↵
- 63.↵
- 64.↵
- 65.↵
- 66.↵
- 67.↵
- 68.↵
- 69.↵
- 70.↵
- 71.↵
- 72.↵
- 73.
- 74.↵
- 75.↵
- 76.↵
- 77.↵
- 78.
- 79.↵
- 80.↵
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.↵
- 93.↵
- 94.↵
- 95.
- 96.↵
- 97.↵
- Received for publication November 2, 2011.
- Accepted for publication March 29, 2012.