Dynamic susceptibility contrast (DSC)-magnetic resonance imaging (MRI) noninvasively measures brain tumor cerebral blood flow (CBF) and cerebral blood volume (CBV), and it has found increasing clinical applications for patient management (1–18). To facilitate multi-institutional comparability and consistency, national initiatives, including National Cancer Institute's Quantitative Imaging Network, Radiological Society of North America's Quantitative Imaging Biomarkers Alliance, and the National Brain Tumor Society's Jumpstarting Brain Tumor Drug Development Coalition, are underway to standardize acquisition and analysis protocols for DSC-MRI (19, 20). A challenge to such efforts is the relative paucity of studies systematically evaluating the influence of DSC-MRI methodology on CBV accuracy. In practice, such validation studies are difficult to perform in patients because of the need for multiple contrast agent (CA) injections and lack of a noninvasive gold standard CBV measure for reference. As an alternative to in vivo validation, in silico digital reference objects (DROs) provide a means for computing synthetic MRI signals and derived kinetic parameters for a range of clinically relevant input conditions. Such a DRO was recently developed for dynamic contrast-enhanced MRI to investigate the biases and variances of algorithms used for image analysis (21).
The goal of this report is to describe the development of a DSC-MRI DRO that recapitulates the heterogeneous signal characteristics measured in glioblastomas. In general, there are two underlying strategies that can be pursued for DROs emulating MRI data. When the primary objective is to establish multisite analysis consistency, synthetic signals can be computed using simple heuristic models approximating the underlying biophysics of signal formation, as the endpoint is to assess the agreement between software estimates of a parameter such as CBV that is explicitly defined by the “ground truth” time course. However, if the intention is to optimize acquisition protocols and CA dosing schemes, such as those used in DSC-MRI, or if the accuracy of the analysis is dependent upon certain physical or physiological assumptions, the synthetic signals should accurately reflect the biophysics of the MRI signal. For the DSC-MRI DRO, we pursued the latter strategy because it enables a more accurate and comprehensive investigation into the DSC-MRI methodology.
In brain tumor DSC-MRI, the acquired signals reflect a complex combination of T1, T2, and T2* changes that depend upon numerous features including CA kinetic parameters (CBF, permeability, intra- and extravascular volume fractions), precontrast T1 and T2*, vascular architecture, cellular microstructure (size, shape, spatial distribution), transvascular and transcellular water exchange, and CA T1 and T2* relaxivity. The sensitivity of the DSC-MRI signal to relaxation time variations is influenced by the acquisition parameters (repetition time [TR], echo time [TE], flip angle [FA], pulse sequence type) and CA dosing scheme (preload and bolus dose and timing). Accordingly, for the DRO to yield realistic signals, its design must reasonably approximate the magnitude and heterogeneity of these physical and physiological parameters in vivo. To that end, we developed a DSC-MRI DRO that is driven by a validated computational strategy to compute MRI signals for realistic 3-dimensional (3D) tissue structures (22); partially constrained by parameter inputs defined from in vivo data; and, for unknown parameters, trained using a public database of DSC-MRI data in glioblastomas.
The computational approach used herein, termed the finite perturber finite difference method (FPFDM) (22, 23), models the effects of water protons diffusing in heterogeneous magnetic field medium based on a 3D tissue structure. The FPFDM computes magnetic field perturbations induced by susceptibility variations between the simulated tissue compartments, and it determines the resulting gradient echo transverse relaxation rates. In addition to a 3D matrix that defines the tissue structure (eg, blood vessels and cells), requisite FPFDM inputs include the static magnetic field strength, the CA concentration in each compartment for determining intercompartment susceptibility differences, the water proton diffusion coefficient, and the DSC-MRI pulse sequence parameters. To ensure clinical relevancy, the DRO derived from these input parameters should replicate the magnitude and heterogeneity of CA-induced T1 and T2* changes during bolus passage through vessels and into the extravascular extracellular space (EES).
The 2-compartment pharmacokinetic model described by Brix et al. (24) was used to simulate concentration–time profiles in plasma (Cp) and the EES (Ce). Inputs to the Brix model include vascular volume fraction, blood flow, CA transfer coefficient (Ktrans), and volume fraction of the EES (ve). Rather than use previously reported mean CBF and CBV values in glioblastoma, our simulated kinetic curves better represented clinical data if the DRO voxels matched the paired, voxel-wise distribution of these parameters across patients (as compared with randomly distributed unpaired parameters). Accordingly, we extracted DSC-MRI data from 23 patients with glioblastoma (>40 000 voxels) in The Cancer Imaging Archive (TCIA) database for characterizing the distribution of paired CBF and CBV values. For this patient cohort, DSC-MRI was acquired at 3T, consisting of General Electric (General Electric Healthcare, Waukesha, WI, USA) (n = 14) and Siemens (Siemens Medical Systems, Erlangen, Germany) (n = 9) scanners using single-echo gradient echo-planar imaging (TR = 1–1.25 seconds, TE = 30 milliseconds, FA = 70–80°, field of view = 240 × 240 mm2, section thickness = 4–5 mm, matrix = 962 or 1282) before, during, and after administration of 0.1 mmol/kg gadopentetate dimeglumine (Gd-DTPA) infusion at 4 ml/s followed by a saline flush. Five minutes before bolus injection, a 0.05 mmol/kg Gd-DTPA preload was administered to minimize T1 leakage effects. Residual leakage effects were corrected using the Boxerman–Schmainda–Weisskoff approach (25). Voxel-wise relative CBV and CBF maps were calculated from the leakage-corrected DSC-MRI data and an automated measure of the arterial input function (AIF), using circular singular value decomposition-based deconvolution (26–29). The voxel-wise distributions of Ktrans and ve were characterized using a retrospective analysis of dynamic contrast-enhanced MRI signals extracted from a dual-echo DSC-MRI data set in 11 glioblastomas (30). Because DSC-MRI data yield relative tumor CBV and CBF measures, their values were scaled using data obtained from dynamic computed tomography perfusion imaging (31). In addition, the AIF used as input for the DRO's kinetic modeling was computed as the average AIF among all patients in the TCIA data. Figure 1, A, B, and C shows the average AIF values, CBV and CBF paired distribution, and Ktrans distribution, respectively.
To define a computationally manageable number of tissue models in the DRO that still accurately reflected the in vivo voxel-wise heterogeneity, the 2-dimensional paired distribution of CBF and CBV was first binned into intervals of 5 ml/100 g/min and 1 ml/100 g, respectively. The resulting distribution was then scaled to yield 100 combinations of CBF and CBV pairs, which were then used to define the number and vascular properties of the tissue structures.
Although the component of DSC-MRI signal associated with CA-induced T1 changes is easily calculated by assuming fast water exchange (32–34), the CA-induced T2* changes depend on vascular and cellular microstructural geometry, precluding use of a simple analytical model. To reflect this complexity, we modeled tissue structures using ellipsoids (cells) (22, 23, 35) packed around randomly oriented cylinders (vessels) (36–45). Previously, we showed that modeling cells as ellipsoids rather than spheres provides a more accurate estimate of the magnitude of T2* changes observed in clinical DSC-MRI studies (22, 23, 35), whereas modeling the vasculature structure as randomly oriented cylinders has been shown to accurately estimate the T2* effects that occur when CA is distributed within blood vessels (36–45). The cylindrical vascular volume fraction was fixed using the in vivo extracted CBV values, and vessel sizes varied from 5 to 30 μm (46). Tumor cell volume fractions were allowed to vary within a physiologically relevant range (45%–65%) (47), and the mean cellular axis radii for a given voxel varied between 4 and 15 μm (46). Figure 2 shows a representative 3D volume rendering of 2 tissue structures, one with homogeneous ellipsoids with a constant aspect ratio (Figure 2A) and one showing ellipsoids with heterogeneous shapes (Figure 2B).
Computation of DSC-MRI Signal
The susceptibility differences between the vascular and extravascular compartments were computed using Δχ = χm·[CA], where [CA] is the compartmental CA concentration (Cp and Ce) and χm is the CA molar susceptibility (0.027 × 106 mM−1) (48). In addition to all the aforementioned input parameters, the FPFDM calculates the DSC-MRI signal as described previously (22) using a water proton diffusion rate (D) of 1.3 × 10−3 mm2/s (49), relevant pulse parameters (TE, B0, FA, TR), and precontrast T10 values ranging from 1 to 2.2 seconds. Figure 3 shows representative simulated Cp and Ce time curves (Figure 3A), and the corresponding gradient echo DSC-MRI signal ratio (S/S0) time curves (Figure 3B) for the 2 tissue voxels are shown in Figure 2.
Given the large number of input parameters and a wide range of potential permutations, it is critical to ensure that the DRO's simulated DSC-MRI signals accurately represent the temporal characteristics, magnitude, and distribution of CA-induced T1 and T2* changes observed across typical glioblastomas. To achieve this, we used the voxel-wise TCIA data described above (>40 000 voxels) for identifying the appropriate combination of input parameters. In particular, all computed signals, for an equivalent preload dosing scheme and pulse sequence parameters to those in the TCIA data set, underwent a selection criteria process based on their percent signal recovery (PSR) and the mean and standard deviation of the signals across the DRO. The PSR is a useful metric for comparison because it reflects the magnitude of the signal drop during bolus passage and the postbolus signal recovery. The DRO's input tissue structure (eg, cell size, shape), kinetic parameters (eg, CBF, Ktrans), and physical parameters (precontrast T1) were systematically permutated until the distribution of PSR values and the mean and standard deviation of signals across the DRO agreed with those found in the voxel-wise TCIA data. The PSR agreement was evaluated using a 2-sample Kolmogorov–Smirnov test. In addition, a 95% agreement between the FWHM and the maximum signal drop was used to determine the agreement between the mean signals. To achieve this level of agreement, the iterative process required a DRO consisting of ∼10 000 unique voxels. The data training based on this selection criterion ensured the removal of computed signals from the DRO, because of an unrealistic combination of tissue parameters. Figure 4A–B shows the agreement between the in vivo and in silico mean and standard deviation of the signal. The distribution of PSR values obtained from the training data set and the DRO is shown in Figure 4C. The 2-sample Kolmogorov–Smirnov test yielded a P-value of .69, indicating agreement between the 2 distributions. Table 1 summarizes the final tissue parameter values that were identified through the DRO training.
To validate the DRO's ability to produce reliable signals for pulse sequences and that the CA dosing schemes are different from those in the training data set, we compared simulated dual-echo signals with those found in an in vivo dual-echo DSC-MRI “validation” data set. The validation data set was acquired in patients with glioblastoma (n = 3) at 3T using a dual gradient echo-planar imaging protocol with the following parameters: TR = 1.5 seconds, TE1/TE2 = 7.0/31.0 milliseconds, field of view = 240 × 240 mm2, section thickness = 5 mm, matrix = 962. Measurements were taken before, during, and after administration of Gd-DTPA (0.1 mmol/kg Gd-DTPA, 4 ml/s infusion rate followed by 20 ml of saline flush). In the simulation, the structural and kinetic inputs derived during the training phase remained the same, but the acquisition parameters and dosing scheme were chosen to match those used in the patient data. The goal of this validation study was to determine whether the DRO fully captures the heterogeneity (eg, magnitude and temporal characteristics such as PSR) of the DSC-MRI signals acquired in this separate (and smaller) cohort of patients. To identify this subset of voxels within the DRO, a correlation analysis was performed between the signals in the in vivo and DRO data. The range of PSR values found in the in vivo and DRO data was compared to ensure that the DRO captured the signal heterogeneity measured in the validation set for both TEs. A parameter termed percent relaxation drop (PRD) was also formulated in a similar fashion as PSR using the derived dual-echo ΔR2* time courses and compared between the in vivo and DRO data.
All simulations were performed using Matlab (MathWorks, Natick, MA) running on a high-performance 32-core system with 2.3 GHz processors and 128 GB of RAM.
Figure 5 compares simulated and in vivo dual-echo DSC-MRI data. The DRO could accurately recapitulate the TE = 7 milliseconds and TE = 31 milliseconds signals and the derived dual-echo ΔR2* time courses, which remove T1 leakage effects but retain T2* leakage effects. The PSR and PRD heterogeneity of the in vivo data was also fully reflected in the DRO. This indicates that the trained DRO can accurately model the underlying CA-induced T1 and T2* effects and the associated DSC-MRI signals for different sets of pulse sequence parameters and CA dosing schemes.
Application 1: Influence of Acquisition and Postprocessing Methods on CBV Accuracy
It is well established that T1 and T2* CA leakage effects confound the reliable measurement of CBV (25, 50). DSC-MRI acquisition strategies have been proposed to reduce T1 leakage effects, including the use of preload CA administration, low FAs, long TEs and TRs, and dual-echo pulse sequences. In addition, postprocessing methods have been developed that eliminate residual T1 and/or T2* leakage effects (25, 51–59). However, validation of these acquisition and postprocessing strategies in vivo has been limited because of the lack of a reliable gold standard reference. A potential application of the population-based DRO is the systematic investigation of the acquisition and postprocessing methods that influence the reliability of CBV measurements.
To this end, we computed the percentage difference between tumor CBV simulated with and without (Ktrans = 0) CA leakage effects for a single-dose bolus injection protocol (no preload), FA = 30° and 90°, TE = 30 milliseconds, and TR = 1 and 2 seconds. We also compared CBV accuracy with and without the application of postprocessing leakage correction using the Boxerman–Schmainda–Weisskoff approach. Results of this analysis are shown in Figure 6. For TR = 1 second, FA = 30° yielded more accurate CBV values than FA = 90°, with and without postprocessing leakage correction (Figure 6A). As expected, the uncorrected 90° FA data yielded substantially underestimated CBV across the DRO voxels, reflecting the strong sensitivity to T1 leakage effects. For TR = 2 seconds, a greater fraction of voxels overestimate CBV, indicating a shift toward T2*-dominated leakage effects due to reduced T1 sensitivity (Figure 6B). Leakage correction improved CBV accuracy across all acquisition parameters. A similar approach could be used to systematically investigate the influence of a range of acquisition and postprocessing methods on CBV accuracy.
Application 2: In Silico Optimization of DSC-MRI for Use in Clinical Trials
The population-based DRO can also be used to optimize DSC-MRI for assessment of treatment response in clinical trials. For example, the influence of acquisition and postprocessing methods on the sensitivity of DSC-MRI to a given CBV change can be used to determine protocols that minimize the sample size needed to power a clinical trial. In this context, the DRO serves as an atlas of possible tumor DSC-MRI signals. By using the correlation analysis discussed in the validation section, a virtual patient DSC-MRI data set can be generated by replacing voxel-wise in vivo tumor signals with an atlas-matched version. This analysis can be propagated across an existing clinical trial database to compute in silico pre- and post-treatment DSC-MRI data. Because the simulated signals for a given voxel originate from a unique set of input conditions, the DSC-MRI signals can be recomputed for any combination of acquisition parameters, such as a new FA or CA dosing scheme. This permits systematic investigation of how acquisition and postprocessing methods influence the inter- and intrasubject CBV heterogeneity, pre- and post-therapy. Alternatively, an assumed effect size distribution (eg, 20% ± 5% decrease in a tumor's CBV) could be applied to the untreated cohort of virtual patients and can be used to identify, within the DRO, the “treated” DSC-MRI signals for each voxel.
Figure 7A–B illustrates a simulated pretreatment CBV map for a virtual patient computed using 2 different CA dosing schemes: a single-bolus dose with no preload (method 1) and a single-dose preload preceding a single-bolus dose (method 2). The corresponding treated CBV maps (modeled as a 20% mean reduction in tumor CBV) for both methods are shown in Figure 7C–D. The pre- and post-treatment CBV distributions across the entire tumor region of interest for both acquisition methods are shown in Figure 7E–F. In this example, CBV estimates derived from method 2 were more sensitive to treatment response compared with those derived from method 1, as indicated by the change in CBV. Similar analyses could be extended to cohorts of virtual patients to identify the most robust and sensitive DSC-MRI acquisition and postprocessing strategies for use in clinical trials.
We have described the development of a DRO that recapitulates the DSC-MRI signal characteristics observed in human glioblastomas. The DRO enables signals to be computed for ranges of physiological, physical, and acquisition parameters. Clinical relevance is ensured through the use of a training data set. Furthermore, we validated the DRO's ability to produce reliable signals for different CA dosing schemes and acquisition parameters. Although in silico models may be limited by the accuracy of the biophysical model used, they provide a feasible and robust alternative to in vivo studies, which, in the case of DSC-MRI, may require multiple contrast injections and MRI scans and often lack a reliable “ground truth” for establishing accuracy.
Two key features of the proposed DRO are instrumental to its ability to provide signals that emulate clinical data. First, the DSC-MRI signals are derived using a validated computational approach that enables the incorporation of realistic tissue structures. Unlike heuristic models of DSC-MRI (34), this approach does not make assumptions regarding the voxel-wise CA T2* relaxivity, a parameter that is highly dependent upon vascular and cellular microstructure. In the proposed DRO, the voxel-wise microstructure determines the compartmental volume fractions and the associated CA relaxivity. Second, the training phase ensures that the range of simulated signals reflects the heterogeneity observed in vivo. Without training, there is the potential to introduce bias into the optimization of acquisition and postprocessing methods, as such methods may have not have uniform accuracy across the range of parameters.
Although we have presented two potential applications for the proposed DRO, there exist numerous opportunities for its use. Studies seeking to characterize and explore the biophysical basis of DSC-MRI data in brain tumors have yielded new biomarkers sensitive to the underlying tumor microstructure (eg, morphological features of vessels and cells) (23, 60–62) and hemodynamics (eg, vascular architectural imaging) (63). For these advanced methods, the DRO provides a tool with which to systematically investigate the sensitivity of DSC-MRI to such features and identify optimal acquisition protocols. Furthermore, the DRO can also be used to assess the accuracy of kinetic parameter estimates derived from newly developed pulse sequences, such as the recently proposed multiecho spin and gradient echo (SAGE) approach (64–69).
Although we trained the DRO with and validated it against in vivo data, any simulation approach that models complex biophysical phenomena has limitations. As described previously (22), the current computational approach does not consider the effects of arbitrary or heterogeneous CA distribution within a given tissue compartment such as the EES. The DRO could also be expanded to include the effects of transvascular water exchange rate, intravascular flow dynamics, atypical cellular geometries, and more heterogeneous vascular tree models. However, increasing the biological complexity of the input tissue structures also increases the number of unknown parameters that would need to be characterized.
The proposed DSC-MRI DRO provides a tool that can be leveraged by groups aiming to optimize and standardize acquisition and analysis methods for prospective clinical studies. It also enables the evaluation of bias and variance introduced by multisite data analysis. Such efforts are critical for establishing comparability of DSC-MRI data and interpreting multisite clinical trial data. To facilitate this effort, a range of DSC-MRI DROs is available for download from The Cancer Imaging Archive (www.cancerimagingarchive.net) under the collection name Barrow-DRO. The provided files contain multiple versions of the DRO, computed across a wide range of pulse sequence parameters and preload dosing schemes, all saved in Digital Imaging and Communications in Medicine (DICOM) and Matlab formats. Table 2 summarizes the range of pulse sequence parameters and CA dosing schemes that, when combined, yield 360 different acquisition methods. Each DRO file is a DSC-MRI time series data set similar to what would be acquired clinically and includes predefined regions of interest for the AIF, normal tissue and tumor voxels. Accordingly, these data may be processed using commercial or customized DSC-MRI analysis packages. The data set summary page details the organization of the files, the regions of interest, and the instructions for use.
|TR (ms)||FA (°)||TE (ms)||B0 (T)||Preload + Bolus|
|(1000, 1500, 2000)||(30, 60, 90)||(20, 30, 40, 50)||(1.5, 3)||(0 + 1, 1/4 + 3/4, 1/2 + 1/2, 1/2 + 1, 1 + 1)|