ROCKETSHIP: a flexible and modular software tool for the planning, processing and analysis of dynamic MRI studies

Validation of ROCKETSHIP software with simulation

Precision and accuracy of individual kinetic model fitting

The robustness of ROCKETSHIP was evaluated with simulations. Simulated datasets containing
two (Tofts, Patlak models), three (Extended Tofts) and four parameters (two-compartment
exchange model, 2CXM) were generated using MATLAB. The ability of ROCKETSHIP to recover
the appropriate DCE-MRI parameter values from simulated data at different SNR and
time resolutions was evaluated. To achieve this, curves were generated using specific
parameter values defined in Table 2. For each parameter permutation Rician noise at different SNRs was added to the signal
intensity versus time curves and fitted. This process was repeated with new noise
100 times in a Monte-Carlo manner to estimate the accuracy and precision of the fits.
Resultant curves were fitted with ROCKETSHIP using the same model that was used to
generate the curves. The population AIF published by Parker et al. was used in all
simulations 15]. Default fitting options were used. The accuracy of the fitted parameters were evaluated
with the concordance correlation coefficient (CCC) 45].

Table 2. Fitting parameters for simulation studies

Figure 3 shows plots comparing fitted K_trans values to the values used to generate the simulated curves for different models.
At high SNR and short time resolution, there is good concordance between simulated
and fitted values. This is reflected in the CCC comparisons (Tables 3, 4 and 5). Fitting using ROCKETSHIP was generally able to recover parameter values with high
accuracy and precision, as demonstrated by the small error bars seen in Fig. 3 and the high CCC values for the fitted to actual value comparisons across multiple
models and with different parameters. Results shown in Fig. 3 and Tables 3, 4 and 5 demonstrate that ROCKETSHIP perform on par or better compared to QIBA simulation
fittings using the Tofts model with DCE@UrLAB and DCEMRI.jl both qualitatively (compared
to figures generated by DCE@UrLAB 33]) and quantitatively (CCC???0.9 in general for K_trans and v_e recovery using Tofts model, similar to or better than the reported CCCs generated
by DCEMRI.jl 28]). As with all curve fitting algorithms, the accuracy of the fit will be dependent
on a number of factors 18], including SNR, sampling resolution and the number of parameters being fitted (Fig. 4, Tables 3, 4 and 5).

Fig. 3. K_trans fitting of simulated data. Simulated data with time resolution of 0.5 s and SNR?=?100
were fitted using the same model used to generate the simulation with ROCKETSHIP using
default settings for the Patlak method (a), Tofts (b) and Extended Tofts models (c). K_trans simulated vs. fitted were plotted as a function of v_e and v_p. Dashed line is unity. Error bars denote standard deviation. Given the similar fits,
points for different v_e and v_p may overlap. Concordance correlation coefficients for these (and other model fits)
are shown in Tables 3, 4 and 5

Table 3. Concordance correlation coefficients (CCC) comparing fitted and simulated Ktrans using
different models and dependent parameters

Table 4. Concordance correlation coefficients (CCC) comparing fitted and simulated vp using
different models and dependent parameters

Table 5. Concordance correlation coefficients (CCC) comparing fitted and simulated ve using
different models and dependent parameters

Fig. 4. v_e fitting at different time resolutions. Simulated data using the Tofts model were
generated at SNR?=?5 and at time resolutions of 0.5 s (a) and 6 s (b). Simulated vs. fitted v_e were plotted as a function of K_trans. Dashed line is unity. Error bars represent standard deviation. As expected, lower
time resolution results in a high standard deviation of the curve fits. Given the
similar fits, points for different K_trans may overlap. Concordance correlation coefficients for these (and other model fits)
are shown in Tables 3, 4 and 5

Accuracy of model selection using nested model analysis

The data-driven, nested model analysis function in ROCKETSHIP was also tested. Simulation
curves reflective of each level of nesting (steady-state, Patlak and extended Tofts
models) were generated as above. The resultant curves were fitted using the nested
method.

Simulation results using the nested model fitting are shown in Fig. 5 and Table 6. The nested model fitting was able to recover the appropriate model for the majority
of the fitting curves. Accuracy of the nested model fitting was shown to be dependent
on SNR. Lower SNR led to more model misclassification of voxels.

Fig. 5. Nested model selection from simulated data. a and b show fitting for steady-state model simulated data. c and d show the fitting for Patlak simulated data. All the generated curves at SNR?=?100
converged to the correct model. At lower SNR, some of the curves incorrectly converged
to Model 3 (extended Tofts). e and f show fitting on extended Tofts simulated data. Again, the majority of the curves
converged to the correct model. The percentage of voxels attributed to each model
by the nest model algorithm is shown in Table 6

Table 6. Model selection of simulated data using the nested model method

In vivo examples

Preclinical data set

DCE-MRI was performed on an athymic BLABL/c nude mouse bearing a HER2-expressing BT474
human breast cancer tumor. A 17^?-estradiol pellet (Innovative Research of America) was implanted subcutaneously into
the back of the mouse one day prior to orthotopic inoculation (6^th and 9^th mammary fat pad) of 5?×?10⁶ cells in matrigel (BD Biosciences). The tumor volume was???200 mm³ at the time of imaging. Mouse care and experimental procedures were carried out in
accordance with protocols approved by the Institutional Animal Care and Use Committee
at Caltech.

Imaging was performed using a Biospec (Bruker-Biospin Inc. Billerica, MA) 7 T MRI
scanner and a custom-built birdcage coil. The mouse was anesthetized during the session
with 1.3-1.5 % isoflurane/air mixture. Body temperature was maintained at 36-37 °C
with warmed air through the bore. A variable flip angle method was used to generate
T1 maps. Gradient echo images (FLASH, FA?=?12°, 24°, 36°, 48°, 60°, matrix size?=?140?×?80,
voxel size?=?0.25?×?0.25 mm², slice thickness?=?1 mm, TR/TE?=?200/2 ms) centered on the tumor (3 slices) and the
left ventricle (1 slice) were acquired. Next, DCE-MRI was acquired (FLASH, FA?=?35°,
TR/TE?=?25/2 ms, geometry the same as the T1 maps, time resolution?=?2 s, duration?=?22 min).
After a baseline of 2.5 min, Gd-DTPA (0.1 mmol/kg, Magnevist, Bayer) was injected
intravenously via a tail vein catheter using a powered-injector (New Era Inc.) at
0.5 mL/min. Imaging data were processed using ROCKETSHIP and analyzed with the nested
model method.

Parametric maps and associated concentration vs. time curves are shown in Fig. 6. For this tumor, the extended Tofts model was deemed the most appropriate for the
majority of the tumor. Generated K_trans values were similar to those observed in prior studies using the same tumor cell
line 46]. Consistent with other studies, there is a heterogeneous distribution of K_trans, v_e and v_p highlighting a leaky and vascularized rim with a necrotic core.

Fig. 6. Nested model fitting of DCE-MRI data on a murine breast cancer tumor model. Parameters
for K_trans (a), v_e (b), and v_p (c) are shown. As shown in d, the majority of the voxels fitted best to the extended Tofts model, with some edge
voxels fitting to the Patlak method. e shows the AIF used for the fit (taken from the left ventricle). f shows a sample time curve from the edge of the tumor (denoted by arrow) with corresponding
fit (blue denotes the fit, red lines denote the 95 % prediction bounds for the fitted
curve). Rod phantoms on either side of the mouse were present to allow for signal
drift correction (not used in this case)

Clinical data set

The clinical data set was obtained from a wider study set being accrued by the University
of Southern California (USC) Alzheimer’s disease Research Center. The study was approved
by the USC Institutional Review Board. The imaging dataset used here was obtained
from one of the participants from the no cognitive impairment cohort, as determined
from medical examination and neuropsychological evaluation. All imaging was performed
at the Keck Medical Center of USC. Participants underwent a medical examination, neuropsychological
evaluations and blood draw to ensure appropriate kidney function for CA administration
prior to imaging.

The imaging protocol performed was developed to study BBB changes in patients with
cognitive impairment and is detailed in 6]. Briefly, all images were obtained on a GE 3 T HDXT MR scanner with a standard eight-channel
array head coil. Anatomical coronal spin echo T2-weighted scans were first obtained
through the hippocampi (TR/TE 1550/97.15 ms, NEX?=?1, slice thickness 5 mm with no
gap, FOV?=?188 x 180 mm, matrix size?=?384 x 384). Baseline coronal T1-weighted maps
were then acquired using a T1-weighted 3D spoiled gradient echo (SPGR) pulse sequence
and variable flip angle method using flip angles of 2°, 5° and 10°. (TR/TE?=?8.29/3.08 ms,
NEX?=?1, slice thickness 5 mm with no gap, FOV 188 x 180 mm, matrix size 160 x 160).
Coronal DCE MRI covering the hippocampi and temporal lobes were acquired using a T1-weighted
3D SPGR pulse sequence (FA?=?15°, TR/TE?=?8.29/3.09 ms, NEX?=?1, slice thickness 5 mm
with no gap, FOV 188 x 180 mm, matrix size 160 x 160, voxel size was 0.625 x 0.625
x 5 mm³). This sequence was repeated for a total of 16 min with an approximate time resolution
of 15.4 s. Gadobenate dimegulumine (MultiHance®, Braaco, Milan, Italy) (0.05 mmol/kg),
a gadolinium-based CA, was administered intravenously into the antecubital vein using
a power injector, at a rate of 3 mL/s followed by a 25 mL saline flush, 30 s into
the DCE scan. Imaging data were processed by ROCKETSHIP. The AIF, which was extracted
from an ROI positioned at the internal carotid artery, was fitted with a bi-exponential
function prior to fitting analysis with 2CXM.

Parametric maps and associated concentration vs. time curves are shown in Fig. 7. Compared to the murine tumor that contains leaky vasculature 41], K_trans values for the human brain are lower, as expected given the intact blood brain barrier
in normal human subjects 6]. As shown in Fig. 7e, the bi-exponential function estimated the AIF well and allowed for a good fit of
concentration time curves within the brain parenchyma (Fig. 7f).

Fig. 7. 2CXM fitting of a normal human brain. Parameters for K_trans (a), v_e (b), v_p (c) and F_p (d) are shown. e shows the AIF used (taken from internal carotid artery). The AIF was fitted with
a bi-exponential curve (blue) prior to tissue fitting. f shows a sample time curve from the brain parenchyma (denoted by arrow) with corresponding
fit (blue denotes the fit, red lines denote the 95 % prediction bounds for the fitted
curve). Fold over artifact is seen on the lateral brain edges due to truncation by
the field of view bounding box