A non-contrast multi-parametric MRI biomarker for assessment of MR-guided focused ultrasound thermal therapies

A. MRI Acquisition and Clinical Biomarkers

All experiments were carried out in accordance with the approved Institutional Animal Care and Use Committee regulations at the University of Utah (Protocol 17–08012, approved 09/07/2017). Following intramuscular injection of VX2 tumors cells (1×10⁶ cells in 50% media/Matrigel solution) into the quadriceps of four New Zealand white rabbits (2.5–3 kg), tumors were grown to approximately 2 cm in length, then ablated using an MRgFUS system (Image Guided Therapy, Inc.) with a 256-element phased-array transducer (Imasonic, Voraysur-l’Ognon, France; 10-cm focal length, 14.4×9.8 cm aperture, f=940 kHz). Animals were intubated and anesthetized with isoflurane (2–3%) and monitored for vitals throughout the procedure. Ablation procedures were performed inside a 3T MRI scanner (Prisma^FIT Siemens, Erlangen, Germany) for tumor targeting, treatment monitoring, and post-treatment assessment. The MRgFUS system incorporated a single-loop MR receive coil around the targeted quadriceps to improve the MR image signal-to-noise ratio. All MR acquisition parameters are listed in Table I.

The rabbit was positioned on the treatment table as shown in Figure 1. High-resolution (0.5-mm isotropic) non-contrast T1w images of the entire treated leg (marked with ** in Table I) were acquired prior to and immediately after ablation for MR registration purposes. Coronal MPMR images of the lower leg were acquired immediately before (T2w, DWI), during (3D MRTI, PRF method), and 20 minutes after (T2w, DWI) the ablation procedure. The ablation procedure targeted 50% of the VX2 tumor and the surrounding muscle tissue. Sonication details for each subject are outlined in Table II. Approximately 40 minutes following the final sonication, CE-T1w MR images were acquired immediately following intravenous gadolinium injection (0.3 mL/kg ProHance) and animals were recovered. Due to the expected transient effects of ablation such as edema and latent apoptosis, which occurs for up to 72 hrs [4], [14], we defined the final treatment effect as the necrotic volume at 3–5 days after treatment. At this time-point post-ablation, the animal was re-positioned on the MRgFUS table and T1w and CE-T1w images were acquired for longitudinal MR and histological registration purposes. Animals were then immediately euthanized for tissue excision and histology processing. A large region of the quadriceps which fully encompassed the tumor and MRgFUS ablation volume was excised and fixed in formalin for two weeks before sectioning and staining with H&E. Further details of gross tissue processing, sectioning, and staining with H&E are reported in the Appendix.

ADC maps were computed using a mono-exponential fit of the DWI signal at 2 b-values (b=20,500), averaged in three directions. Semi-automatic segmentation of the acute NPV was performed by thresholding the CE-T1w image intensity to encompass the non-enhancing region, followed by manual editing. All acute NPV segmentations were evaluated by an expert radiologist (NW) and scored using a Likert scale of 1–5. Expert-evaluated acute NPV segmentations were considered the clinical NPV prediction of tissue necrosis.

B. MRTI Prior Baseline Reconstruction

MRTI data was reconstructed and zero-fill interpolated to 1-mm isotropic resolution. A prior baseline approach similar to Bitton et al. [29] was implemented to account for localized temperature accumulation resulting from multiple sonications in the targeted tissue. This prior baseline approach is an improvement over the current clinical method for calculating MR thermal dose, which assumes a constant baseline temperature throughout the treatment and underestimates the total dose. In summary, the baseline phase from a prior sonication was used to calculate phase accrual at the start of subsequent sonications. A temporal criterion was applied to ensure valid baseline sonications were chosen to account for potential large ${\text{B}}_0$-field drift and subsequent phase-wrap errors. Therefore, the prior baseline sonication, ($S_p$), whose pre-heating baseline phase image is used as the baseline phase image for calculating the temperature change in the subsequent sonications, was reset to the current sonication when the time between sonications exceeded 10 minutes.Ten minutes is the simulated necessary cooling time for the muscle to return to a 1°C rise after a sonication reaching the maximum temperature elevation observed in this study. The details of the simulation implementation are described in [30], using a 0.5-mm iso-tropic T2w-anatomical segmented model of the rabbit hindlimb (Seg3D2, SCI Institute [31], University of Utah) and tissue properties from [32]). Based on the time between subsequent sonications, $S_p$ was reset 1, 2, 1, and 5 times for animals 1–4, respectively. For each sonication, $S_n$, the local phase accumulation was calculated as:

where ${\varphi }_{0,n}$ is the baseline phase image for $S_n$, ${\varphi }_{0,p}$ is the baseline phase image for $S_p$, and $\text{Δ}{\varphi }_{BA,n}$ is the bulk phase accrual resulting from ${\text{B}}_0$-field drift and inter-scan resonance frequency adjustments. $\text{Δ}{\varphi }_{BA,n}$ was calculated for each coronal slice of the MRTI volume by taking the mean value of $\left({\varphi }_{0,n}-{\varphi }_{0,p}\right)$ in a 15×15 voxel ROI in non-heated quadriceps muscle. For each sonication $S_n$, the phase change due to FUS heating and $\text{Δ}{\varphi }_{acc,n}$ were converted to temperature using the PRF method ($\alpha =-0.01\text{ppm}{/}^{\circ }\text{C}$) [8]. The absolute temperature profile for $S_n$, was calculated as the sum of these components and the baseline body temperature as measured by rectal fiber-optic probe monitoring animal body temperature at the time of the prior sonication $T_{FO,p}$:

where $T$ is temperature (°C), for $k$ MRTI measurements in sonication $n$. To account for thermal dose accumulation between sonications, all time points between the previous and ${\text{n}}^{\text{th}}$ sonication, $S_n$, were padded with the value given by $\text{Δ}{T}_{{\varphi }_{acc,n}}$. Finally, the CTD in each voxel was calculated using the CEM₄₃ model for each sonication and summing across all sonications [5].

C. Histological Registration

To create a co-registered data set for supervised classification training and validation, there were three registration steps: 1) motion during treatment 2) longitudinal changes in subject pose and position between ablation imaging and follow-up imaging [33], and 3) registration between histology and follow-up imaging. The registration methods used in this study have been previously developed [33] [25] but are briefly summarized here. During the MRgFUS procedure, the rabbit hindlimb may have slowly sunk deeper into the water bath, introducing minor motion errors over the course of the 3-hour treatment. To correct for this bulk motion in MRTI data, the final MRTI sonication magnitude image was registered to the post-treatment non-contrast T1w image using variance equalization and highly constrained (to limit deformation to bulk motion) elastic registration. Each of the prior MRTI sonication images was then registered to the deformed final MRTI sonication image using the same elastic registration method. The motion estimated by this registration process was also applied to the pre-treatment MR images and visual inspection was used to ensure the tissue-water boundary was aligned across the images. The bulk motion estimated during the registration process was less than 5 mm for all subjects.

Between treatment day imaging and 3–5 day follow-up imaging, there were unavoidable changes in the animal positioning in the MR scanner. A volume-conserving longitudinal registration algorithm was used to register follow-up images to treatment day images as tissue volume is preserved under normal physiological loading [33]. High-resolution T1w images (marked with ** in Table I) without contrast were used for registration since the features resulting from treatment are not visible. This registration step provided an invertible 3D diffeomorphism between the treatment day images and the follow-up images, which are registered to histology in the following step.

Finally, direct, voxel-wise comparison between ground truth histology and treatment day MR images required accurate histology registration. Histology necrosis annotations were registered to in vivo MPMR images to generate a volumetric label of histology necrosis. Briefly, a novel workflow was used to correct for deformation from every step of the tissue processing for histology. 3D surface registration was used to generate an invertible 3D diffeomorphism for each histology slice in the in vivo MR images. Stacking the individual histology annotations of treated tissue provided a histology-derived volumetric label of the treated volume, co-registered with the follow-up images. See the Appendix for more details on 3D histology-to-MR registration. The diffeomorphisms from longitudinal registration and histology registration were composed, resulting in a voxel-wise alignment between the treatment day T2w images, ADC maps, MRTI data, and the histology necrosis ground-truth label.

D. MPMR Classification Models

To develop the supervised MPMR biomarkers, the dataset comprised all voxels in 3D anatomical segmentations of the quadriceps muscles (Seg3D2, SCI Institute [31], University of Utah) from all four subjects. The number of voxels per subject is given by ”Quadriceps Vol.” in Table II. The inclusion of the entire quadriceps for supervised learning resulted in 12.3% positive class voxels on average, as assigned by the 3D registered histology label ({C0: viable, C1: necrotic}). The number of positive class voxels (1-mm isotropic) per subject is given by ”Histology Vol.” in Table II. Four non-contrast MR image volumes were used for classifier inputs: 1) CTD, 2) maximum temperature projection in time (MTP), 3) T2w images, and 4) ADC maps (Figure 2). Overlapping in-plane 3×3 neighborhoods around each voxel were extracted from each of the four features, resulting in 36 total input features. For classifier training, input features were normalized to the mean and scaled to unit variance (scikit-learn, StandardScalar), with training and test sets normalized separately. First, optimal classifier hyper-parameters were identified with 5-fold cross-validation in the combined dataset from all subjects. Second, the hyper-parameter tuned classifiers were trained and validated in a leave-one-out (LOO) approach, where classifiers were trained on three subjects and validated on the fourth, for a total of n=4 folds. Details of classifier development are below.

Three binary classifiers were explored for the development of the voxel-wise MPMR biomarker: a logistic regression classifier (LRC; ”liblinear” solver [34]), a support vectors machine classifier (SVMC; SVC implementation, ”rbf” kernel [34]), and a random forest classifier (RFC; [34]). The LRC model was chosen as an option for a model with minimal complexity, while SVMC with radial basis functions and RFC were implemented in order to learn possible nonlinear relationships in the MPMR feature space. Random forests, with the capacity to create models with high complexity, are also prone to over-fitting data. The SVMC algorithm is also capable of higher complexity than logistic regression; however, is less susceptible to over-fitting in small datasets than RFCs.

For hyper-parameter tuning, the full dataset (N=67,235 voxels) was implemented for a 5-fold cross-validation grid search of classifier hyper-parameters. Folds were generated with a 70:30 random split, maintaining the original dataset imbalance (scikit-learn, StratifiedShuffleSplit [34]). To correct for class imbalance in the dataset, a balanced scorer was used for training the LRC and SVMC models, and balanced sampling was used for each decision tree batch for training the RFC model. Optimal hyper-parameters were selected to minimize the classifier complexity to reduce over-fitting while maximizing the 5-fold average Dice score (f1-score) in each classifier.For the LRC, an l1-penalty with C-parameter=0.05 was selected.For the RFC, the selected hyper-parameters were min samples-split = 0.005, n estimators = 50, and max depth = 10. The selected hyper-parameters for the SVMC were C-parameter = 0.01 and gamma = 0.0001.

To test the accuracy and generalizability of the MPMR biomarkers, the classifiers with optimized hyper-parameters were re-trained in a LOO strategy. For the LOO analysis, the training data consisted of voxels from of three subjects and the validation dataset consisted of voxels from the fourth subject. Optimal thresholds for the output probabilities were selected to maximize the Dice score in the training data. Trained models and optimal thresholds were fit to the validation dataset to obtain probabilities for each voxel. LOO training was repeated for each validation subject for a total of n=4 folds. Similarly, the optimized thresholds for the clinical CTD biomarker were chosen to maximize the Dice score in the LOO training data sets. The binary NPV segmentation is dependent on a user-defined threshold between the hypo-intense NPV and the surrounding hyper-intense rim. For comparison to numerically continuous biomarkers, the NPV inter-user variability was simulated by applying a Gaussian blur to the expert-evaluated NPV segmentations (3-voxel kernel, $\sigma =1$ voxel), generating a continuous 3-mm boundary varying from 0 to 1. For the clinical NPV and 240 CEM₄₃ metrics, thresholds of 0.5 and 240 CEM₄₃ were used respectively.

E. Biomarker Evaluation

Precision, recall, and Dice scores were compared for all MPMR and clinical biomarkers for each LOO fold. The Dice metric was chosen to score overall biomarker performance given the severe imbalance of the dataset (11% positive class). This metric provides a harmonic mean of the recall and precision scores. To assess prediction accuracy with more clinically relevant metrics, in each subject the percent difference in volume and the mean-distance-to-agreement (MDA) between the predicted thermal necrosis lesion and the histological necrosis label was calculated for each trained biomarker. See Appendix for MDA definition.

A. Histological Registration

The acute NPV segmentations were representative of clinician segmentations, achieving an average score of 4.75 out of 5 in the four subjects by the expert radiologist. The volume-conserving longitudinal MR registration error was 1.26 ± 0.52 mm and the average volume change of tissue during registration was 0.24 ± 0.12 percent relative to the starting volume. The average total registration error between acute MR imaging and histology was 1.00 ± 0.13 mm.

To demonstrate MR-to-histology registration, the “delayed” NPV acquired 10 minutes prior to euthanasia 3 days post-ablation is compared to the registered H&E in Figure 3 in subject 4. The delayed NPV segmentation on CE-MRI (Figure 3b) correlates to the region of coagulative necrosis in the registered H&E section (Figure 3a), which is characterized by pale eosin staining (pallor) of muscle and tumor tissue. A volumetric comparison of the longitudinally-registered 3D histology necrosis label (pink), delayed NPV (green), and acute NPV (blue) is shown in Figure 3c. In 2D cross-section profiles of the volumes in each dimension, the pink necrosis contour aligns well with the in vivo delayed NPV segmentation. As expected, the acute NPV, acquired 40 minutes post-ablation, is larger than the delayed NPV and necrosis label.

B. MPMR Leave-One-Out Classifier Training

MPMR classifiers with optimized hyper-parameters were trained in a LOO approach with four folds. All voxels from one subject served as the test dataset for each fold.The average precision-recall (PR) curves across all LOO folds for the MPMR biomarkers are shown in Figure 4 and compared to the clinical NPV (blue) and CTD metric (red) curves. All MPMR classifiers out-perform the CTD model, with the SVMC and LRC predictors surpassing a Dice score of 0.6.

In Figure 5, the training (dashed) and testing (solid) Dice scores are compared in each subject across MPMR classifiers. All MPMR classifiers achieved higher Dice scores than the optimized CTD in all subjects. Although RFC achieved the highest Dice score in the full dataset, the average RFC Dice score (0.58, Table III) was the lowest of all MPMR classifiers in the LOO training. Additionally, the test Dice score was on average 25.9% lower than the training scores for the RFC indicating potential data over-fitting. In contrast, testing scores were only 8.0% and 4.0% lower than the training scores for the LRC and SVMC models on average. Therefore, trained LRC and SVMC models generalized well to new subjects in each fold.

Precision, recall, and Dice scores for all LOO folds for all clinical and MPMR biomarkers are summarized in Figure 6 and Table III. Acute NPV, the current gold-standard MR biomarker of thermal ablation, achieved the highest average Dice score of 0.68 ± 0.10 (n=4), followed by the SVMC (0.63 ± 0.05) and LRC (0.63 ± 0.04), then RFC (0.58 ± 0.15). MPMR classifiers were slightly more precise than the acute NPV at the cost of reduced recall scores. Interestingly, the RFC achieved the highest average precision of all clinical and MPMR biomarkers of 0.62 ± 0.21.

In subjects 1 and 2, the LRC improved Dice over acute NPV by 2.2 and 8.9%, respectively; and in subject 2, SVMC improved Dice over acute NPV by 1.8%. However, acute NPV had the greatest Dice in subjects 3 and 4. In comparison to the optimized CTD, all three MPMR classifiers improved Dice scores in subjects 1–4 by 941.4%, 7.8%, 46.7%, and 13.7% on average across n=3 classifiers. Furthermore, the intersubject variability in Dice scores, as measured by the standard deviation in Table III was the greatest for the 240 CEM₄₃ and CTD models (0.27 and 0.26, respectively) and lowest for the SVMC and LRC models (0.05 and 0.04, respectively).

C. Optimized CTD Leave-one-Out Training

The optimal CTD thresholds calculated in each fold were 119.1, 120.0, 26.3, and 218.5 CEM₄₃, which are all lower than the accepted 240 CEM₄₃ threshold. The 240 CEM₄₃ and optimized CTD thresholds achieved similar average Dice scores of 0.43 ± 0.27 and 0.42 ± 0.26, respectively. LOO training for optimized CTD led to variable results, where Dice increased by 0.02 on average in 3/4 subjects and decreased by 0.09 in subject 3 (Table III). CTD optimization resulted in higher recall and lower precision as the CTD thresholds in the LOO training datasets were reduced. Inconsistent train-test score differences across folds (Figure 5) indicate that optimized thermal dose thresholds in training subjects do not generalize consistently to new ”unseen” subjects.

D. Spatial Accuracy

Spatial performance was assessed in the test data of each subject in the corresponding LOO fold. Representative examples of clinical and MPMR model predictions in subject 4 are shown in Figure 7. The volume of intersection between the label and prediction is depicted by the opaque green volume, which is positively correlated to the Dice score. The average absolute value of the percent differences between predicted and labeled necrosis volumes were 59.5, 54.7, and 50.0% for the acute NPV, 240 CEM₄₃ and CTD biomarkers, respectively. The absolute differences were 34.3, 30.4, and 28.8% for the LRC, SVMC, and RFC biomarkers, respectively (Table III).

High recall of the acute NPV is due to the consistent overestimation of the histological necrosis boundary (Figure 7a). Conversely, the 240 CEM₄₃ and optimized CTD predictions did not cover the full extent of thermal damage (Figure 7b), or the optimized CTD threshold severely overestimated the necrotic volume in the LOO fold for subject 3. The LRC and SVMC prediction contours (Figure 7c–d) both align more closely with the optimized CTD contour than the RFC (Figure 7e), albeit with higher sensitivity than CTD. The RFC provides the most conservative prediction of the necrotic volume and the least variance in volume difference errors across folds (± 4.5%).

The 5-mm grid in Figure 7 demonstrates the magnitude of contour errors for all predictions. These errors were quantified with the mean-distance-to-agreement (MDA) metric in the whole volume in each subject and are summarized in Table III. The 240 CEM₄₃ and CTD metrics had the greatest MDA error across all subjects, ranging from 1.4–6.5 mm. Acute NPV MDA error was lowest overall with a range of 1.2–1.3 mm for all subjects. Of the MPMR classifiers, MDA was lowest for RFC (1.4–1.7 mm). The SVMC and LRC classifiers were similar with MDA errors ranging from 1.2–2.0 mm.

A. Limitations

The volumetric histological processing utilized in this study was costly and intensive, limiting the number of subjects and amount of data for this analysis. However, a quantified comparison of delayed NPV to registered histological outcomes may allow follow-up CE-MRI to serve as a surrogate label for supervised learning in the future. Due to the limited data acquired for this study, a voxel-wise approach to machine learning was implemented in this study. A 3×3 voxel neighborhood patch was used for each input feature to provide the models with contextual information and promote robust learning despite small registration errors (<3 mm). Although a 2D patch was used, the prediction from the classifiers was only for a single voxel. Voxel-wise predictions tend to be noisy and cannot fully leverage spatial connectivity and context to improve the prediction. For small sample sizes, predictions could be improved using connected components; however, this method requires knowing how many different locations have been treated, which may vary depending on near- and far-field heating. However, larger datasets from more subjects would allow for CNN or deep learning applications which could produce more connected 2D or 3D maps that are less susceptible to image noise. These complex algorithms may also learn post-treatment differences in multiple tissue types, increasing the likelihood of applying transfer learning to several applications.

Sources of error related to MR registration to histological outcomes are multi-fold. Registration errors of MR and histological registration were on the order of an MRTI voxel (1-mm isotropic); however, errors can vary across and within subjects and directly impact classifier scores. The longitudinal MR registration has been rigorously validated in Zimmerman et al. [33]. Although localized tissue swelling after ablation was not directly accounted for, anatomical landmarks near the region of the ablation, such as blood vessels, were utilized to determine target registration error. Single-voxel shifts in pre- and post-ablation registered ADC and T2w maps can introduce errors in difference map calculations. However, treatment effect classifiers often rely on pre- and post-treatment differences as inputs. Pre-processing with context-based radiomic filters or implementation of a CNN may reduce the impact of noisy or imperfectly registered MPMR data.