Retrospective analysis of an IRB-approved database of patients treated on our MRg system was done to identify patients with liver lesions treated with MRgSBRT. Patients were included if they had a biopsy-proven primary HCC, known history of a primary malignancy with biopsy-proven metastasis, or consensus to treat after review in a multidisciplinary tumor board. Patients were excluded from analysis if they had synchronous primary malignancies, did not complete MRgSBRT, or did not have post-treatment imaging for restaging. Data abstracted included demographics, performance status, history of cirrhosis, history of hepatitis C infection, prior radiofrequency ablations, prior radiation treatments, prior chemotherapy regimens, prior trans-arterial chemoembolization procedures, pre-treatment complete metabolic panels, pre-treatment complete blood counts, tumor size, tumor location in the liver, histology, radiation dose delivered, and radiation fractions delivered.
Image Acquisition and Target Delineation
Low field strength setup MR images acquired on a 0.35 T split bore unit using a balanced steady-state free procession pulse sequence were retrieved from the MR linac treatment planning system. The image acquisition protocol was similar for all images acquired during treatment; breath-hold images were acquired in 25 seconds with voxel dimensions of 1.5 x 1.5 x 3.0 mm3. Two torso coils with six-phased array elements were placed hemi-circumferentially around the patient. Gross Tumor Volume (GTV) was delineated by a radiation oncologist resident and reviewed by a board-certified radiation oncologist experienced in MRgSBRT and low field MR setup images. MRI safety screening were administered to all eligible patients. Treatment response of restaging scans was assessed using RECIST and mRECIST criteria, depending on initial histology.
Image Analysis, radiomic feature selection, and model construction
Image DICOM files along with the associated GTV structure sets were exported to a local drive and radiomics texture features were extracted using the Texture Analysis Toolbox in MATLAB 2020b (The MathWorks, Natick MA, USA). Binary masks of the GTVs were generated to extract three-dimensional bitmaps from the RTStruct DICOM objects. The dynamic range of intensity values was constrained via “±3σ” Collewet normalization method 23. A total of 39 second-order features were extracted from the GTVs of the daily MR setup images. Feature classes utilized were gray level co-occurrence matrix (GLCM), gray level size zone matrix (GLSZM), gray level run length matrix (GLRLM), and neighborhood gray tone difference matrix (NGTDM) 24–28. The GLCM code/aggregation code for calculation of the GLCM was LFYI/IAZD with the feature name 8ZQL when using the Image Biomarker Standardization Initiative (IBSI) nomenclature 29. The IBSI code/aggregation code for the GLSZM is 9SAK/KOBO and the feature name and formula is 48P8. Delta-radiomic texture features were calculated as the change in radiomic features between the initial setup image (before treatment) and after BEDs of 20 Gy and 40 Gy were delivered, producing BED20 and BED40 feature libraries. These quanta of radiation were selected to account for dissimilar radiation doses and fractionation schedules. For example, a delta-radiomics texture feature value calculated from a hypothetical patient with a BED/Fx = 10 Gy would be calculated for the BED20 with the texture features calculated from the image acquired prior to fraction 1 and the setup image from fraction 3 (delivery of 2 treatments).
Model construction and feature selection
A two-step method was used where the first step was designed to rank features in order of predictive potential, while the second step was used to evaluate the top features. In step one, features were ranked using the Gini Index, calculated during the training phase of a random forest (RF) model. The Gini Index is a measure of how important each predictive variable is for accurate prediction 30. A RF model was built using delta-radiomics as predictors with 500 decision trees considering six features (‘mtry’ parameter) with replacement one at a time. The ‘mtry’ parameter was selected as approximately the square root of the number of delta-radiomics texture features under consideration, thirty-nine total for each timepoint. The RF prediction model selects subsets of the available data for predictors, constructs a simple decision tree, then evaluates the accuracy of predictors included and excluded from each decision tree. This is performed for each tree in the forest. The importance of the variables in the model for prediction is calculated and recorded as the Gini Index. Finally, the two features with the highest Gini Index were selected for further investigation.
Those two features were used to estimate the internal validation area under the receiver operating characteristic curve (AUC) using a bootstrapped logistic regression model. The model was bootstrapped 1,000 times randomly sampling 2/3 of the patient data at each iteration. The internal validation AUC was obtained by feeding in the remaining 1/3 of the patients and predicting response. The AUC was recorded for each iteration. This process was performed for both BED20 and BED40 features that were selected by the Gini Index. The mean AUC, 2.5 percentile, and 97.5 percentile values were calculated using the 1,000 iterations. Two additional random forest models were constructed, a strictly clinical model, and a combined model (delta-radiomic BED20 features with clinical data), to compare the importance of clinical variables and BED20 delta-radiomic texture features using the Gini Index. The clinical model included BED/fraction, maximum tumor dimension, sex, ethnicity, initial Karnofsky performance status score, and existence of prior cirrhosis as variables for predicting treatment response. The AUC estimates produced by the random forest model during training was recorded for evaluation. The process was repeated for the combined model with the top two BED20 delta-radiomics features included to assess the importance of the factors relevant to available clinical information. Both the clinical model and combined models were constructed with 500 trees. The mtry parameter was set to include all data (mtry = 6 for the clinical model and mtry = 8 for the combined model).
All methods were carried out in accordance with relevant guidelines and regulations.
The study was approved by The University of Miami’s institutional IRB committee under IRB#20160817. Informed consent was waived due to the retrospective nature of the study by The University of Miami’s institutional IRB.