Non-invasive intravoxel incoherent motion MRI in prediction of histopathological response to neoadjuvant chemotherapy and survival outcome in osteosarcoma at the time of diagnosis

Background Early prediction of response to neoadjuvant chemotherapy (NACT) is important to aid personalized treatment in osteosarcoma. Diffusion-weighted Intravoxel Incoherent Motion (IVIM) MRI was used to evaluate the predictive value for response to NACT and survival outcome in osteosarcoma. Methods Total fifty-five patients with biopsy-proven osteosarcoma were recruited prospectively, among them 35 patients were further analysed. Patients underwent 3 cycles of NACT (Cisplatin + Doxorubicin) followed by surgery and response adapted adjuvant chemotherapy. Treatment outcomes were histopathological response to NACT (good-response ≥ 50% necrosis and poor-response < 50% necrosis) and survival outcome (event-free survival (EFS) and overall survival (OS)). IVIM MRI was acquired at 1.5T at baseline (t0), after 1-cycle (t1) and after 3-cycles (t2) of NACT. Quantitative IVIM parameters (D, D*, f & D*.f) were estimated using advanced state-of-the-art spatial penalty based IVIM analysis method bi-exponential model with total-variation penalty function (BETV) at 3 time-points and histogram analysis was performed. Results Good-responders: Poor-responders ratio was 13 (37%):22 (63%). EFS and OS were 31% and 69% with 16.27 and 25.9 months of median duration respectively. For predicting poor-response to NACT, IVIM parameters showed AUC = 0.87, Sensitivity = 86%, Specificity = 77% at t0, and AUC = 0.96, Sensitivity = 86%, Specificity = 100% at t1. Multivariate Cox regression analysis showed smaller tumour volume (HR = 1.002, p = 0.001) higher ADC-25th-percentile (HR = 0.047, p = 0.005) & D-Mean (HR = 0.1, p = 0.023) and lower D*-Mean (HR = 1.052, p = 0.039) were independent predictors of longer EFS (log-rank p-values: 0.054, 0.0034, 0.0017, 0.0019 respectively) and non-metastatic disease (HR = 4.33, p < 10–3), smaller tumour-volume (HR = 1.001, p = 0.042), lower D*-Mean (HR = 1.045, p = 0.056) and higher D*.f-skewness (HR = 0.544, p = 0.048) were independent predictors of longer OS (log-rank p-values: < 10–3, 0.07, < 10–3, 0.019 respectively). Conclusion IVIM parameters obtained with a 1.5T scanner along with novel BETV method and their histogram analysis indicating tumour heterogeneity were informative in characterizing NACT response and survival outcome in osteosarcoma. Supplementary Information The online version contains supplementary material available at 10.1186/s12967-022-03838-1.


Background
Neoadjuvant and adjuvant chemotherapy along with surgery have significantly improved the long-term disease-free survival rate of patients with osteosarcoma [1].
Five-year disease-free survival rate for localised osteosarcoma is reported as 60-70%, whereas it is less than 20% in patients presenting with distance metastasis [1], although the survival is relatively lower in resource challenged countries [2,3]. Measurement of necrosis by histopathological evaluation (HPE) of resected lesion is the current gold standard for neoadjuvant chemotherapy (NACT) response assessment and an important prognostic factor for survival in patients with localized osteosarcoma [4]; however, it is only possible after surgery on completion of NACT. There is no robust, non-invasive and early assessment tool for evaluation of prognosis. Early evaluation of treatment response and prediction of survival outcome may help to prevent the patients from undergoing ineffective chemotherapy regimen and triage them to alternate therapeutic options saving time, cost and side-effects [4][5][6][7]. Qualitative evaluation of morphological MR images for intensity variation and tumour size changes, for non-invasive NACT response evaluation in solid tumours, do not always correlate well with the HPE [8,9]. Therefore, the development of non-invasive radiomics based strategies to select patients with poor prognosis who might benefit from alternative therapy is desirable [10].
Recently developed spatial penalty based IVIM analysis method Biexponential model with Total variation penalty function (BE + TV or BETV) [44] has shown qualitatively and quantitatively improved parameter estimation compared to the commonly used IVIM analysis methods like biexponential model [27] and its segmented variants [28][29][30][31][32][33][34][35][36]. Robustness of this novel, BETV method has already been demonstrated in cancer simulations & various clinical applications such as, osteosarcoma [39], Ewing sarcoma [40], lymphoma [41], brain tumour [42] and prostate [43] for characterizing tumour and measuring treatment response. Existing studies on DWI performing prediction of NACT response and treatment outcome in osteosarcoma have shown inconclusive or contradictory findings [11] and there is scope for further research. Therefore, this study evaluates the applicability of this state-of-the-art, BETV, IVIM analysis method for oncological applications using clinical datasets of osteosarcoma.
Further, histogram analysis assessing heterogeneity in tumour microenvironment has shown significant improvement in tumour characterisation and predicting the therapeutic response with a more direct correlation with the underlying structural and pathophysiological changes manifested upon tumour progression [44]. Novelty of this study lies in exploring the role of quantitative IVIM analysis for predicting long-term survival outcome in patients with osteosarcoma that has not yet been performed. In this prospective study, assessment of quantitative IVIM diffusion and perfusion parameters and their histogram analysis has been performed to identify  [46,47]. After completion of treatment, routine follow-up evaluation was performed every 3 months for the first 2 years and every 6 months for the subsequent years. For the purpose of this study follow-up data were collected till 31st December 2020. In addition to clinical evaluation, follow-up imaging examination consisting of chest radiographs alternating with NCCT chest every 3 monthly for a total of 5 years.

Histopathological response evaluation
Pathologist, blinded to the clinical status and MRI results, analysed resected lesions, described tumour-size, extent and amount (in percentage) of necrosis relative to the whole tumour volume. Response to NACT was assessed histologically according to the six-grade scale of Salzer-Kuntschik et al. [43]. Patients were categorised into two groups-good-response (patients with ≤ 50% viable-tumor) combining Grade I-IV and poor-response (patients with > 50% viable-tumour) combining Grade V-VI patients as Salzer-Kuntschik grading [43].

Survival outcome evaluation
Endpoints studied were event free survival (EFS) and overall survival (OS). An event was defined as the elapse of secondary tumours or distant metastasis or local recurrence, or death from any cause. EFS was defined as the time interval from the first day of chemotherapy to any of the events or to the last date of follow-up, whichever is first. OS was defined as the time interval from the first day of chemotherapy until death. Patients who were alive without any event at the time of the last follow-up were censored and they were included both in EFS and OS.

MRI acquisition protocol
MRI was acquired at three time-points-baseline or pre-NACT (t0), after the 1st NACT cycle (t1, 2-3 weeks) and after the 3rd NACT cycle (t2, 8-9 weeks). MRI acquisition was performed using a 1.5T Philips Achieva ® MR scanner with phased-array surface coil or an extremity coil. Conventional T1-weighted, T2-weighted and IVIM-DWI sequences were acquired according to the standard MRI acquisition protocol. T1-weighted and T2-weighted images were acquired using the Turbo-Spin-Echo sequence with TR/TE = 528/10 ms and 3797/60 ms respectively, matrix-size = 512 × 512 and 384 × 384 respectively. Acquisition of IVIM-DWI was performed using free-breathing spin-echo echo-planar imaging (SE-EPI) with a variation of gradient strengths of

Quantitative imaging parameters
Tumour-volume (in cc) at three time-points was determined using region of interest (ROI) drawn manually by an expert radiologist (D.K., > 12 years of experience in cancer imaging) covering whole tumour on b = 800 s/ mm 2 DWI images with reference to the morphological T1W and T2W images. Quantitative IVIM parameters Diffusion coefficient (D), Perfusion coefficient (D*) and Perfusion fraction (f) and ADC were evaluated in whole tumour-volume at three time-points t0, t1 and t2. IVIM parameters were evaluated using the state-of-the-art IVIM analysis method BE model with adaptive Total Variation (TV) Penalty function (BETV method) [39]. BETV method applies non-linear least-square (NNLS) optimisation for data fitting with adaptive penalty function Total Variation for reconstruction with good SNR and reduces the non-physiological spatial inhomogeneity in estimated parametric images. The product of relative micro-vascular flow and volume (D*.f) was also calculated voxelwise in tumour-volume and analysed as it provides the information about vascular changes in terms of relative microvascular perfusion or blood flow, analogous to the blood-flow as measured in perfusion imaging [48]. ADC was calculated by a mono-exponential model using b-values ≥ 200 s/mm 2 assuming perfusion effect is negligible at higher b-values [27]. Goodness-of-fit (R 2 ) and Coefficient-of-variation (CV) were calculated as the measure of precision and reproducibility in IVIM parameter estimation respectively. Histogram analysis of imaging parameters was performed in tumour-volume at three time-points. Eleven (n = 11) histogram parameters: mean, standarddeviation(SD), skewness, kurtosis, energy, entropy, 90th, 75th, 60th, 50th and 25th percentiles and their relative percentage changes between time-points t0-t1(ΔI) and time-points t0-t2(ΔII) were calculated for each patient. Quantitative parameters evaluation and histogram analysis was performed using an in-house built toolbox in MATLAB ® (MathWorks Inc., v2017, Philadelphia, USA).

Clinical parameters
Tumour-volume, alkaline phosphatase (ALP), lactate dehydrogenase (LDH) at baseline, and classical prognostic factors like primary tumour site (axial/pelvic vs peripheral), presence of metastatic disease at diagnosis were evaluated for their potential impact on chemotherapy response and survival outcome as EFS and OS.

Statistical analysis
Inter-group (between good-response and poor-response) statistical significance (p < 0.05) of clinical parameters and absolute histogram parameters of ADC, D, D*, f & D*.f and their relative percentage changes (ΔI & ΔII) were evaluated using independent sample t test. Intra-group significant (p < 0.05) changes in parameters across timepoints were tested using paired t-test. Predictive performance of statistically significant (p < 0.05) parameters for NACT responsiveness was assessed using receiver-operating-characteristic-curve (ROC) analysis at time-points t0, and t1.
Univariate Cox regression analysis was used to assess the effects of the statistically significant (p < 0.05) clinical and imaging parameters on EFS and OS using hazard ratio (HR). Significant histogram parameters derived from ADC, D, D* & f were tested separately for multicollinearity using variance inflation factor (VIF), while VIF ≥ 8 indicated high collinearity. Using Harrells's c-index [49] and the corresponding generalization of Somers' D rank correlation [50] (SDRC), the parameters with the most discriminative ability for EFS and OS was selected to develop the multivariate cox proportional hazard model in combination with significant clinical parameters. Higher values of C-index and SDRC indicated better discriminative ability. Final multivariate survival model(s) was tested for assumption of proportionality using Schoenfeld test. Kaplan-Meier curves were evaluated for the parameters showing statistical significance (p < 0.05) after multivariate analysis, and differences were assessed by using a log-rank test. Before inclusion in the Cox analysis, the distributions of the continuous parameters were examined for normality and a multiplicative transformation (× 10 3) was applied on extreme observations (mean, SD, energy, 90th-25th percentiles of imaging parameters having values in the order of 10 -3 ) as reported earlier [51]. Analysis for EFS and OS was performed in all patients at time-point t0.
At t1, 90th percentile of ADC and 90th-25th percentile values of D were significantly (p < 0.02) lower among good-responders than poor-responders. Histograms of ADC and D became more negatively skewed among good-responders than poor-responders (ADC-skewness: Illustrative example of IVIM parametric maps and corresponding histograms in tumor volume of representative patients with tumor involving different anatomical regions like femur, tibia and humerus are presented in Figs. 2, 3 and 4 respectively. No significant qualitative or quantitative reginal differences was observed in the IVIM parametric maps evaluated in different anatomical regi ons. "

Chemotherapy response prediction
Clinical parameters like tumour-volume, ALP, LDH, primary tumour site and metastasis were not found to be statistically significant (p > 0.13) between good-response and poor-response groups. Statistically significant imaging parameters and their ROC curve analysis for chemotherapy response prediction is presented in Table 3.  with Sensitivity = 86% and Specificity = 77% in predicting poor-response at t0 (Fig. 5a). At t1, mean of imaging parameters ADC,D,D*,f &D*.f jointly showed AUC = 0.92 in predicting poor-response; while in combination with statistically significant (p < 0.05) histogram parameters ADC-90th-percentile and D-90th-25th-percentiles produced AUC = 0.96 with Sensitivity = 86% and Specificity = 100%; in predicting poor-response to NACT (Fig. 5b).

Survival outcome prediction
Baseline parameters from univariate and multivariate Cox regression analyses which had significant effects on EFS and OS are presented in Table 4. Univariate analyses showed, clinical parameters such as metastasis, tumourvolume and ALP were found to be significantly associated with EFS and OS.  Table 2 Average Apparent Diffusion Coefficient (ADC), Diffusion coefficient (D), Perfusion Coefficient (D*), Perfusion fraction (f) and D*.f in good-response (GR) and poor-response (PR) groups at baseline (t0), after 1st cycle of chemotherapy (t1), and after 3rd cycle of chemotherapy (t2) and the relative percentage changes of parameters across time-points t0 & t1 (ΔI) and t0 & t2 (ΔII)    Fig. S6a, b). Both the models OS-Model-1 and OS-Model-2 produced comparable

Discussion
This prospective study evaluates the role of non-invasive quantitative IVIM imaging in characterizing tumour microenvironment, predicting chemotherapy response and long-term survival outcome of osteosarcoma at baseline and early in the course of treatment. Histogram analysis of IVIM-DWI parameters showed, microvascular perfusion and its heterogeneity in tumour were significantly higher among poorresponders than good-responders characterizing higher angiogenic changes among poor-responders, thus effectively predicted poor-response to NACT at baseline (t0, AUC = 0.87). While after the first cycle of NACT, diffusion parameters along with D*.f predicted poor-response to NACT with high AUC (t1, AUC = 0.96).
For survival outcome, multivariate Cox regression analysis showed smaller tumour volume, higher levels of ADC-25th percentile, D-Mean and lower D*-Mean as the independent predictor of EFS; while nonmetastatic disease, smaller tumour volume, lower levels of D*-Mean and higher D*.f-skewness were independent predictors of OS. These findings suggest quantitative IVIM parameters and their histogram analyses may be useful for characterizing and quantifying heterogeneity in tumour micro-environment and thereby predicting chemotherapeutic response and outcome in osteosarcoma.
IVIM imaging is influenced by the diffusion of free water molecules in the intra & extra-cellular compartments and micro-circulation of water molecules in micro-capillaries. The observations in this study support that cell density and abnormal/immature micro-vessels may decrease in osteosarcoma after chemotherapy resulting an increase in ADC and D and reduction in D* and f similar to the previous studies in osteosarcoma [11-15, 37, 38] and other tumours like colorectal [30], head and neck [31], cervical [32], nasopharyngeal [33,34], breast [35], hepatic [36] and other tumours [28,29]. Table 3 Statistically significant (independent sample t test, p < 0.05) histogram parameters of ADC, D, D*, f and D*.f among goodresponse and poor-response groups at baseline (time-points t0) and after 1st cycle of chemotherapy (time-point t1) and their ROC curve analysis for predicting poor-response to chemotherapy  At baseline, similar mean ADC and D values were observed among both the response groups possibly due to high heterogeneity of cellularity in osteosarcoma, similar to the previous studies [11,14,15,37]. After 1st NACT cycle, mean ADC and D values in tumour were significantly higher among poor-responders than the good-responders; however, after an initial increase in ADC and D at t1 (ΔI≈29%↑), diffusion did not increase any further (ΔII≈ΔI≈29%↑) among poor-responders. Necrotic tumours with large extracellular space (with higher mean and percentile values of ADC&D) are often associated with poor response to therapy [12,31,34] due to hypoxia and tissue acidosis that leads to resistance to chemotherapy [4,52]. Whereas, for good-responders, ADC and D both were observed to be increased throughout all NACT cycles (ΔI = 20%↑, 22%↑ respectively and ΔII = 31%↑, 33%↑ respectively) indicating possible increase in cell death resulting in an increase in diffusion of water molecules in the tumour. At baseline, higher levels of diffusion parameters were observed to be associated with improved EFS in osteosarcoma patients similar to earlier studies in other tumours [28,29,36].
IVIM perfusion related parameters (D*, f, D*.f) correlate with the process of angiogenesis and reflect changes in relative microvascular perfusion, perfusion volumefraction and flow in tumour respectively [48]. Lee et al. showed that, D* and f were significantly correlated with the micro-vessel density score in murine model colorectal cancer; providing information about tumour perfusion and angiogenesis [53]. In this study, a higher heterogeneity in micro-perfusion pattern in osteosarcoma was observed as the markers of therapeutic poor response similar to previous studies [18,19,25,26,44]. On the other hand, comparatively a higher reduction in D* and D*.f during chemotherapy and lower heterogeneity in tumour microvasculature among good-responders indicated relatively lower angiogenic progression. Lower levels of perfusion parameters were associated with improved EFS and OS in osteosarcoma patients similar to earlier studies in other tumours [28,29,36]. Analysis showed, perfusion parameters (D* & f) had higher predictive values (AUC = 0.6-0.8) than diffusion parameters (ADC & D) (AUC = 0.6-0.7) at baseline and after the 1st cycle of chemotherapy in predicting NACT response.
ADC value has been widely used for assessing chemotherapy response and survival outcome in osteosarcoma [11][12][13][14][15][16][17]. However, in this study, IVIM parameters in combination with ADC showed improved prediction performance for chemotherapy response than ADC alone (t0, AUC = 0.77 vs. 0.61; t1, AUC = 0.92 vs. 0.71). Measurement of true tissue diffusion was observed to be necessary and useful for characterizing chemotherapeutic changes in osteosarcoma, similar to other IVIM studies in literature [27,28,30,31]. The ability of IVIM to characterise early changes in microvascular perfusion along with true diffusion in tumour is highly relevant in this era of anti-angiogenic chemotherapy drugs which can be performed without the use of exogenous contrast agent; however, reliability and reproducibility should be ensured. Widely used BE model [27,32] and segmented-BE techniques [30,31,33,35] evaluates IVIM parameters at each voxel independently, overlooking the spatial context that may lead to unreliable solutions resulting in noisy parametric image reconstruction, especially for perfusion related parameters [29,54]. Thus, adding a physiologically plausible spatial constraint to the existing BE model is expected to provide a reliable parametric estimation as shown by Baidya Kayal et al. [39]. This recently developed BETV method [39] incorporates gradient based penalty TV [55] with NNLS optimisation of the BE model to preserve the desired spatial homogeneity in the parametric images. Robustness of this method has been shown earlier in both cancer simulations & clinical cohorts of osteosarcoma [39], Ewing sarcoma [40],  6 Kaplan-Meier survival curves demonstrate differences in patient outcome groups at a cut-off value by log-rank test for event free survival (EFS) and overall survival (OS). a Tumor volume for EFS; difference was p = 0.064 at a cut-off value of 240.9 cc tumor volume. b ADC-25th percentile # for EFS; difference was p = 0.0034 at a cut-off value of 1.01. c D-Mean # for EFS; difference was p = 0.0017 at a cut-off value of 1.11. d D*-Mean # for EFS; difference was p = 0.0019 at a cut-off value of 27.24. e Metastasis for OS; significant difference was p < 10 -3 by log-rank test. f Tumor volume for OS; the significant difference was p = 0.07 at a cut-off value of 240.9 cc. g D*-Mean # for OS; difference was p < 10 -3 at a cut-off value of 33.95. h D*.f-skewness # for OS; difference was p = 0.019 at a cut-off of 3.32. # Parameter values were transformed by multiplying with 10 3 before analysis lymphoma [41], brain tumour [42] and prostate [43]. In this study, the-state-of-the-art BETV method was effectively used for analysing IVIM-DWI acquired with a 1.5T scanner and able to provide potential imaging biomarkers for NACT response and survival outcome in osteosarcoma with satisfactory results.
There are a few limitations of our study. First, among the available pathological scales to categorise chemotherapeutic response groups, in this study, patients were categorised into good-response (≥ 50% HPE-necrosis) and poor-response(< 50% HPE-necrosis) groups with reference to the earlier study by Picci et al. [56]. Secondly, perfusion or functional imaging, as DCE MRI or FDG-PET/CT, would have been beneficial to characterise angiogenesis changes after treatment and validate the findings from IVIM perfusion parameters; however, we could not do contrast MRI due to financial and time constraints considering the fact that our study protocol involved three time-point imaging for each patient. Thirdly, advanced texture analysis might be helpful, but that would expand the scope of current study considerably. Texture analysis and its usefulness can be dealt separately in future studies. Fourthly, different IVIM analysis methods like Bayesian based or stretched exponential methods that may hold relevance for clinical assessment of osteosarcoma; however, these methods were not evaluated in the current study. Finally, as osteosarcoma is a rare tumour, only a limited number of patients (35 patients) could be analysed in this prospective study. Thus, future studies analysing much larger cohort and multi-centric data with standardized MRI protocol could be useful.

Conclusions
In conclusion, clinical parameters such as tumour volume, nonmetastatic disease and ALP can be independent predictors of survival outcome. IVIM diffusion (D) and perfusion-related parameters (D*, f) and their histogram analysis (skewness, entropy, percentile) indicating heterogeneity in micro-vasculature in tumour are useful imaging markers to predict survival outcome at presentation and non-response to NACT before and early in the course of treatment. Therefore, quantitative IVIM analysis with advanced analysis methods can serve as a surrogate marker for characterizing chemotherapeutic response, which can be used for non-invasive monitoring and evaluation of chemotherapy response and treatment outcome in patients with osteosarcoma.