Abstract
Innovative diffusion Magnetic Resonance Imaging (dMRI) models enable in vivo mapping of biologically meaningful properties such as cell size, potential biomarkers in cancer. However, while cancers frequently spread to the liver, models tailored for liver applications and easy to deploy in the clinic are still sought. We tackle this unmet need by delivering a practical and clinically viable liver dMRI modelling framework. Through direct comparison of candidate dMRI approaches in mouse and cancer patients’ data, we select a model of intra-cellular diffusion fitted to highly diffusion-weighted images, as it provides the strongest radiological-histological correlates. We demonstrate the potential application of the proposed model in cancer immunotherapy, stratifying the risk of progression based on baseline cell size and density measurements from dMRI. This result, heretofore unreported and not achievable with standard dMRI indices (e.g., apparent diffusion coefficient), suggests that our approach may become a useful tool for precision imaging in oncology.
Introduction
Routine clinical Magnetic Resonance Imaging (MRI) focusses on visualising macroscopic anatomical features, as presence of tumours. Nonetheless, MRI also offers the possibility of measuring biological properties within each pixel of a three-dimensional (3D) scan – known as voxel. This approach, referred to as quantitative MRI (qMRI)1, involves the acquisition of multiple images, each featuring a different contrast, which are then analysed jointly with a mathematical model. qMRI provides promising metrics, which could become quantitative biomarkers complementing the qualitative assessment by the expert radiologist2. Diffusion MRI (dMRI) is a qMRI approach that sensitises the signal to water diffusion with magnetic field gradients1,3,4. Since diffusion in biological tissues is influenced by the microenvironment where diffusion takes place, dMRI ultimately enables the indirect estimation of properties at the micrometric scale5, such as the size of cells restricting water6,7. dMRI bridges the gap between macroscopic and microscopic imaging, and has found applications in brain5, spinal cord8, prostate6, breast9 imaging and beyond.
Innovative dMRI techniques are also urgently needed in abdominal imaging, as in liver MRI10,11. The liver is a frequent site for cancer metastasisation 12, and liver tumours are common targets for treatment response assessment in oncology. However, current response criteria such as RECIST13 have limitations, in that they rely on MRI or computed tomography (CT) merely to measure tumour size, without accounting for changes under therapy at the cellular level. Novel dMRI metrics could enable the non-invasive characterisation of cancer microenvironments, shedding light on the composition of tumours that cannot be biopsied. The new readouts could also provide information on tumour heterogeneity, relevant in the development of treatment resistance14,15, and could better stratify patients eligible for treatments such as immunotherapy16, given the challenge of predicting which patients can benefit from these innovative drugs17. This would be a major advancement in oncology, as it may allow for personalised treatment planning, reductions in sample sizes in clinical trials, and ultimately improve patient outcomes18.
The most recent biophysical dMRI techniques describe the non-vascular liver tissue signal as the sum of contributions from intra-cellular and extra-cellular water19–21. While these models provide promising readouts22, their practical use in real-world settings is made unfeasible by i) the high number of dMRI images (and hence long scan time) required to support model fitting, and by ii) the requirement for specialised dMRI acquisitions23, beyond default examinations available in the scanner console. In this study we aim to tackle this unmet need by delivering a practical liver dMRI signal model that is truly feasible in hospital settings, i.e., on 1.5T or 3T systems, with scan time that does not exceed 15 minutes, and using vendor-provided dMRI sequences. With this objective in mind, we embraced the latest “histology-informed” dMRI development paradigm, which is based on informing signal model design with co-localised histology. The framework has shown promise in delivering dMRI approaches with unprecedented fidelity to cytoarchitecture24,25, maximising biological specificity26.
In this article, we aimed to identify a practical mathematicl model that maximises the agreement of dMRI estimates of metrics such as cell size, to their underlying histological counterparts. We analysed a rich data set of dMRI scans and hematoxylin and eosin (HE)-stained images from excised mouse livers and patients’ liver biopsies. We used these data to select the model maximising radiological-histological correlations, corroborating results with computer simulations. Afterwards, we demonstrated the clinical utility of the designed approach in one exemplificatory task, assessing response to immunotherapy in patients suffering from advanced solid tumours – an urgent, unmet need, given the lack of robust predictors of treatment response for this class of drugs17. In summary, our study delivers a liver dMRI approach that offers metrics with high fidelity to histopathology, and which is feasible in the clinic. The proposed method, based on a single-compartment model of restricted, intra-cellular diffusion, fitted to highly diffusion-weighted (DW) images, identified which patients progress faster from baseline dMRI scans. These results, while exploratory and requiring further confirmation, suggest that our dMRI framework could provide complementary information to standard-of-care imaging, and thus play a key role in oncology research and practice.
Results
Overview: data set
Fig. 1 illustrates the data used in this study. We will refer to data obtained in fixed mouse livers as preclinical, while to data obtained in cancer patients as clinical.
Our data set consisted of preclinical and clinical data. The preclinical data encompasses dMRI scans of seven fixed livers from mice (six implanted with tissue from biopsies of patients suffering from prostate cancer; one without any implantation). We scanned the livers ex vivo on a 9.4T system, and obtained HE histological sections at known position. The clinical data includes in vivo liver dMRI scans performed on 33 patients suffering from advanced liver tumours. Scans were performed on clinical 1.5T and 3T MRI systems. For 18 patients, HE-stained material from a biopsy taken from one of the imaged liver tumours was also available. 30 out of 33 patients effectively participated in a phase I immunotherapy trial, and clinical outcome was available as PFS. We used a total sample size of N = 25 for dMRI model design and of N = 30 for response assessment. In the figure, PC standard fort prostate cancer and HCC for hepatocellular carcinoma.
Preclinical data consists of pulsed gradient spin echo (PGSE) DW MRI scans of seven fixed mouse livers, performed ex vivo on a 9.4T Bruker system. It also includes whole-organ HE-stained sections, obtained at known radiographic position. We studied the livers of mice sacrificed as part of xenograft model development in prostate cancer. Six had been implanted with biopsies of prostate cancer patients, while one had not had any implantation. While the livers from the implanted mice did not grow any tumours, they feature a variety of pathologies, with three unique histopathological phenotypes (Fig. S1). The liver from the mouse with no implantation features normal liver structures, and we will refer to it as Control. Of the six implanted cases, two also show normal liver tissue, with normal representation of all hepatic structures. We will refer to these cases as PatNA1 and PatNA2 (patient biopsy implantation, but normal appearing). Another case exhibits generalised necrosis and diffuse acute and chronic inflammation surrounding necrotic areas, with presence of occluded thrombotic vessels. This specimen will be identified as Patnec (patient biopsy implantation, with necrosis). Finally, three specimens feature an immature, lymphoproliferative process, with various degrees of infiltration of small, lymphoid, atypical cells with abundant mitosis, which infiltrate portal vessels and sinusoidal capillaries, but without producing tumours. These will be referred to as Patinf1 to Patinf3 (patient biopsy implantation, with lymphoid cell infiltration).
We obtained clinical data on cancer patients suffering from advanced solid tumours, participating in an ongoing imaging study at the Barcelona Vall d’Hebron Institute of Oncology (VHIO, Spain). The study involves the acquisition of MRI data, alongside clinical and biological information, in patients eligible for a phase I immunotherapy trial. We included data from 33 patients with liver malignancies (mean/std of age: 62.91/12.34 year; 16 male, 17 female). dMRI was based on diffusion-weighted (DW) echo planar imaging (EPI) scans performed on a 1.5T Siemens Avanto system (Twice-Refocussed Spin Echo (TRSE) DW-EPI) and on a 3T GE SIGNA Pioneer system (PGSE DW-EPI). We also obtained digitised HE-stained biopsies from one of the liver tumours. Biopsies were available for 18 patients, and were collected after dMRI, but before immunotherapy. Clinical outcome (progression-free survival (PFS)) was available for 30 patients, who effectively entered the immunotherapy trial.
We used dMRI-histology data to design the dMRI signal model (N = 25) and dMRI-PFS data for immunotherapy response assessment (N = 30).
Overview: dMRI signal models
We studied five dMRI biophysical models, accounting for restricted diffusion inside spherical cells and hindered diffusion in the extra-cellular space6,19,22. Models can be grouped into two families (Fig. 2.A; see Material and Methods).
(A), top: cartoon illustrating the two families of dMRI models considered in this study, consisting of 1) models with no assumption of which of intra/extra-cellular ADC is higher, and 2) models where the extra-cellular ADC is hypothesised to be higher than the intra-cellular ADC. (B), bottom: study overview. We analysed dMRI data from fixed mouse livers (preclinical data) and from cancer patients imaged in vivo (clinical data) to derive estimates of intra-cellular fraction and of cell size. In parallel, we processed histological material from the same tissues (whole-liver sections for the preclinical mouse data; biopsies from one of the imaged tumours for the clinical data), and derived the histological counterparts of such dMRI metrics. We compared dMRI and histological cell size and intra-cellular fraction to select the dMRI model featuring the best fidelity to histology. The utility of the model was then demonstrated in immunotherapy response assessment in vivo. In Fig. 2, pictures from Servier Medical Art have been used. Servier Medical Art by Servier is licensed under a Creative Commons Attribution 3.0 Unported License (https://creativecommons.org/licenses/by/3.0).
The first family is more general, in that it does not make any assumption on which is higher between intra-/extra-cellular ADC (ADCI and ADCE from now on). It includes:
Diff-in-exTD: it accounts for restricted IC diffusion within spherical cells6, and hindered diffusion in the extra-cellular space, with diffusion time dependence (TD)27 in both intra-/extra-cellular spaces28. The diffusion time quantifies the time during which water molecules can sense cellular barriers, before the MR image is acquired.
Diff-in-ex: as previous model Diff-in-exTD, but neglecting TD in the extra-cellular space. Popular techniques such as IMPULSED19 or VERDICT20 are essentially implementations of this model.
Conversely, the second family explicitly assumes that ADCE > ADCI, similarly to related dMRI techniques (e.g., Restriction Spectrum Imaging29; power-law axon radius mapping7). It includes:
Diff-in-exTDFast: equivalent to Diff-in-exTD, ensuring that ADCE > ADCI.
Diff-in-exFast: equivalent to Diff-in-ex, ensuring that ADCE > ADCI.
Diff-in: a model where it is hypothesised that due to fast, extra-cellular diffusion, the extra-cellular signal is negligible, and the measured signal is dominated by intra-cellular water.
Overview: dMRI metrics
All models enable the estimation of volume-weighted mean cell size (vCSMRI, expressed in µm) and intra-cellular signal fraction (FMRI, dimensionless), which can be combined into an apparent cell density per unit volume (
, expressed in cell mm-3)20. For reference, we benchmarked these metrics against routine ADC (in µm2 ms-1) and apparent diffusion excess kurtosis K (dimensionless) from diffusion kurtosis imaging (DKI)30, popular dMRI indices sensitive to cancer cellularity, easy to compute from short acquisitions20,31.
We processed HE-stained histological data with automatic cell detection32 to derive histological counterparts of vCSMRI and FMRI at known radiographic location. The histological metrics were: histological volume-weighted mean cell size (vCShisto, in µm), intra-cellular area fraction (Fhisto, dimensionless), and cell density per unit area (CDhisto, in cell mm-2). We compared vCShisto and Fhisto to vCSMRI and FMRI (Fig. 2B).
A one-pool model of intra-cellular diffusion provides the most histologically meaningful metrics
Fig. 3A summarises the different dMRI models, while Fig. 3B reports values of the MRI-Histology Total Correlation Score (TCS) for all models. TCS measures the overall correlation between histological and radiological readouts of cell size and intra-cellular fraction, and is obtained by summing Pearson’s correlation coefficients between vCSMRI and vCShisto, and between FMRI and Fhisto (see Materials and Methods). Higher TCS point towards stronger histological-radiological correlation. Negative correlations reduce TCS, so they are penalised.
(A): panel summarising the salient differences between the biophysical dMRI models compared in this study. Models can be divided in two families, i.e.: i) models where it is hypothesised that the extra-cellular ADC is higher than the intra-cellular ADC, and ii) models with no hypothesis on which, between intra-/extra-cellular ADC is higher. Violet shades are used for the first family (models Diff-in, Diff-in-exFast and Diff-in-exTDFast), while orange shades for the second family (models Diff-in-ex and Diff-in-exTD). (B): values of TCS for all models, as obtained by fitting them on high b-value images (b > 1800 s/mm2 in the fixed mouse livers; b > 900 s/mm2 in vivo). We performed model selection using a sample size of N = 25 (see Fig. 1). In Fig. 3, pictures from Servier Medical Art have been used. Servier Medical Art by Servier is licensed under a Creative Commons Attribution 3.0 Unported License (https://creativecommons.org/licenses/by/3.0).
The bar plot in Fig. 3B highlights that dMRI models where ADCE > ADCI, shown in violet shades, provide consistently higher TCS values than models that do not make such an assumption (orange shades). We observe the highest TCS for model Diff-in. Note that Fig. 3 refers to TCS values obtained by fitting dMRI models only to high b-value images, as this provided the highest TCS figures. Fig. S2 reports TCS for model fitting performed to the whole set of diffusion images. In this case, TCS is lower, but again, models where ADCE > ADCI provide the highest TCS. In Fig. S2, Diff-in-exFast provides the highest TCS, although this is lower than Diff-in TCS in Fig. 3.
Fig. S3 reports rankings according to additional criteria, namely: the Histology Fidelity Criterion (HFC), measuring the sum of absolute errors in F and vCS estimation via dMRI, and the Bayesian Information Criterion (BIC)33. BIC is commonly used in dMRI model development34,35; it quantifies the overall model fitting quality (penalizing model complexity), but without accounting for histological information. Lower HFC and lower BIC imply better model performance. Fig. S3 reports the number of times, in percentage terms, that a model provides the lowest HFC and BIC across our N = 25 MRI-histology cases. Results essentially confirm rankings seen on TCS: models hypothesising ADCE > ADCI are selected more frequently than models that do not do, according to HFC. The model Diff-in is the most selected model according to both BIC and HFC (fig. S3.B; fitting to high b-value images). Fig. S4 splits HFC and BIC rankings depending on the MRI scanner. In all cases, models with ADCE > ADCI (Diff-in, Diff-in-exFast, Diff-in-exFastTD) are selected more frequently than models Diff-in-ex and Diff-in-ex-TD. When fitting is performed only on high b-value images, Diff-in is the most selected model according to both BIC and HFC.
Computer simulations confirm model selection from MRI measurements
We performed Monte Carlo computer simulations to corroborate the model selection performed on ex vivo and in vivo dMRI data. The simulations consisted in generating synthetic dMRI signals according to the three dMRI protocols used in this study. We synthesised signals for a substrate made of packed spherical cells (Fig. S5), a common body dMRI tissue model6,19,22,23, and then performed model selection on the synthetic signals (see Materials and Methods). Supplementary Tables S1, S2, and S3 report TCS, HFC and BIC rankings from simulated signals. Results confirm that model Diff-in enables the best estimation of cell size and intra-cellular fraction in the vast majorities of cases, thus confirming ex vivo and in vivo MRI results.
Our proposed approach: a one-compartment model of intra-cellular diffusion, fitted to high b-values
In view of all rankings reported above, our recommended modelling approach is the fitting a one-compartment model of restricted, intra-cellular diffusion within spherical cells to high b-values images (≳ 1800 s/mm2 ex vivo, ≳ 900 s/mm2 in vivo) – referred to as model Diff-in.
Cell size and density estimates from the proposed dMRI model correlate with histology
We now report on the Pearson’s correlation coefficient r between histology, Diff-in and DKI metrics. We consider correlations to be weak, moderate, and strong when | r | < 0.4, | r | ≥ 0.4 but | r | < 0.6, and | r | ≥ 0.6.
Table 1 shows that the correlation between Diff-in FMRI (intra-cellular fraction), vCSMRI (cell size index) and CDMRI (cell density) with their counterparts Fhisto, vCShistoand CDhisto are respectively weak (r = 0.19 between FMRI and Fhisto), moderate (r = 0.44 between vCSMRI and vCShisto) and strong (r = 0.70 between CDMRI and CDhisto). The weak correlation between FMRI and Fhisto can be explained, at least in part, with the fact that FMRI is a signal fraction, rather than an actual volume/area fraction (unlike Fhisto), i.e., it is influenced by T2/T1 differences between intra-cellular and residual extra-cellular signals36. Moreover, FMRI estimation can be biased by unaccounted exchange between intra-/extra-cellular water21,37, which is not accounted for in our signal models. Conversely, the much higher correlations between vCSMRI and vCShisto and between CDMRI and CDhisto, point towards the biological specificity of vCSMRI and CDMRI.
The table reports Pearson’s correlation coefficients r and corresponding p-values p of dMRI metrics FMRI (intra-cellular fraction), vCSMRI (volume-weighted cell size index) and CDMRI (cell density per unit volume) with their histological pairs (Fhisto, vCShisto and CDhisto respectively) for the selected dMRI model (Diff-in, fitted to high b-value images). The table also reports correlation coefficients between routine ADC and K from DKI and each of Fhisto, vCShisto and CDhisto. The sample size was N = 25, so that p < 0.05 if | r | > 0.3961. When p < 0.05, grey shadowing is used.
Table 1 also reports correlation coefficients for dMRI ADC and kurtosis K. Both ADC and K exhibit significant, moderate correlations with histological properties, i.e., negative/positive correlation of ADC/K with cell density CDhisto (r = –0.47 and 0.43 respectively) – a result entirely consistent with previous studies38,39. Significant correlations are also seen with Fhisto (r = 0.40, p = 0.048 between K and Fhisto). These findings are in line with known literature: ADC and K are sensitive to the underlying tissue microstructure, but they are also unspecific, being surrogate metrics that conflate different histopathological characteristics into a single number.
Fig. S6 and Fig. S7 show Pearson’s correlation coefficients for all possible pairs of metrics, in the form of correlation matrices. Correlations among dMRI metrics are seen, as a strong negative correlation between CDMRI and vCSMRI (r = –0.84 for model Diff-in-exFast fitted at high b-value). This finding, which indicates that tighter cell packings per unit volume are achieved with smaller cells, appears biophysically plausible, being mirrored by histological CDhisto and vCShisto (r = – 0.88 between CDhisto and vCShisto). Other weak-to-moderate correlations are seen, e.g.: between K and vCSMRI (r = –0.47) and K and CDMRI (r = 0.38), which agree with the correlations observed between K and histological vCShisto (r = –0.31) and CDhisto (r = 0.43). In general, metrics from dMRI models where ADCE > ADCI show stronger correlations with their histological counterparts than models Diff-in-exTD and Diff-in-ex. We observe the strongest dMRI-histology correlations for model Diff-in fitted to high b-value images.
Metrics from the proposed dMRI model reveal intra-/inter-tumour characteristics
Fig. 4 shows maps from the proposed dMRI model Diff-in alongside histological metrics in 3 mouse livers, representative of the 3 phenotypes seen in our mouse data (Control, for normal liver structures; Patinf1, for small cell infiltration; Patnec, for necrosis). Visually, we observe excellent co-localisation between MRI slices and histology sections. The histological details reveal higher cellularity in sample Patinf1 compared to Control, due to packing of small cells in between larger hepatocytes, or an alternation of areas with lower/ higher cell density in sample Patnec. These qualitative trends are confirmed in the histological maps Fhisto, vCShisto, CDhisto, with values in physiologically plausible ranges, as for example intra-cellular fractions around 0.75 and cell sizes of the order of 20 µm40,41. Maps FMRI, vCSMRI and CDMRI replicate the contrasts seeing in their histological counterparts Fhisto, vCShisto and CDhisto. Fig. S8 shows standard dMRI metrics, namely ADC and kurtosis excess K, in the same mouse livers. Visual trends highlight that the higher cell density of sample Patinf1 translated to remarkably reduced ADC and increased K compared to the Control. Lastly, Fig. S9 shows FMRI, vCSMRI and CDMRI maps from another model (Diff-in-exFast). Fig. S10 instead shows maps of the other metrics provided by models Diff-in-exFast and Diff-in (intrinsic cytosol diffusivity D0,I and asymptotic ADCE, DE,∞). The figures highlight that overall, spatial trends seen in maps from the selected model Diff-in agree with those seen in Diff-in-exFast, but metrics from the latter appear noisier. Metrics D0,I and DE,∞ show limited between-sample contrast, and are difficult to validate histologically.
The figure reports MRI and histology data for 3 specimens, representative of the 3 microstructural phenotypes observed in our ex vivo data set, namely: normal liver structures (Control case); a proliferative process, characterized by infiltration of small cells (Patinf1 case); necrosis and inflammation (Patnec case). For all specimens, the following is shown. (A), top left: a high-resolution T2-w anatomical scan is shown next to the corresponding HE section, with histological details. (B), bottom left: histological maps warped to the dMRI space (intra-cellular patch area fraction Fhisto; volume-weighted mean cell size index vCShisto; cell density per unit patch area CDhisto). (D), bottom right: dMRI maps FMRI, vCSMRI and CDMRI from the selected dMRI signal model (model Diff-in, fitted to high b-value images, i.e., b > 1800 s/mm2).
Table S4 reports qualitative per-sample mean and standard deviation of all MRI and histology metrics in mice. FMRI slightly underestimates Fhisto, while vCSMRI slightly overestimates vCShisto. We speculate that the discrepancies may be due, at least in part, to unaccounted factors such as variability in intrinsic cell shape/cytosol diffusivity42 or water exchange21, and by the difficulty of relating accurately 2D histology to 3D MRI43.
Fig. 5 shows FMRI, vCSMRI and CDMRI maps in patients, alongside biopsies. Histopathological assessment highlights the variety of characteristics that can coexist within advanced solid tumours, e.g.: areas of fibrosis; localised areas of tightly packed cancer cells, sourrounded by stromal fibres; necrosis. dMRI FMRI, vCSMRI, CDMRI show contrasts that are plausible with such histopathological features. For example, in a breast cancer liver metastasis in Fig. 5, we observe a core of low intra-cellular fraction FMRI and low cell density CDMRI, compatible with necrosis. In a HCC case instead, we see areas of high FMRI and high CDMRI, sourrounded by lower FMRI and lower CDMRI, potentially indicating the alternation of high cell densities with fibrotic tissue. Fig. S11 shows routine dMRI ADC and K in the same tumours. Spatial trends are also compatible with the histology, e.g., high ADC and low K are seen in the necrotic core of the breast cancer tumour. Supplementary Fig. S12 shows FMRI, vCSMRI and CDMRI from model Diff-in-exFast. Image contrasts match visually those seen in the same metrics from model Diff-in (the proposed approach), giving confidence of the overall robustness of the biophysical modelling framework. Nevertheless, maps appear noisier. Fig. S13 shows intra-cellular cytosol diffusivity D0,I asymptotic ADCE (DE,∞) in the same tumours. Their speckled appearance suggests that these metrics are difficult to measure accurately in vivo19,44.
MRI maps are shown in a biopsied liver tumour in two patients for each MRI scanner, arranged along rows. (A): examples of slices from the high-resolution anatomical T2-w image and from a high b-value image, with biopsied tumour outlined. (C): maps from the selected model (Diff-in, fitted to high b-value images b > 900 s/mm2). From left to right: intra-cellular signal fraction FMRI; volume-weighted mean cell size index vCSMRI; cell density per unit volume CDMRI. (C): histological details from the HE-stained biopsy. For the 1.5T Siemens scanner (first and second rows from top) we report: patient 6 (primary hepatocellular carcinoma) and patient 3 (liver metastases from ovarian cancer). For the 3T GE scanner (third and fourth rows from top) we report: patient 24 (primary hepatocellular carcinoma (HCC)) and patient 30 (liver metastases from breast cancer).
Metrics from the proposed dMRI model stratify the risk of cancer progression in immunotherapy
Finally, we demonstrate the potential utility of the proposed liver dMRI model in an exemplificatory response assessment task. Fig. 6 reports on the PFS stratification based on Diff-in metrics. Panels on the left report results from Kaplan-Meier analyses, log-rank tests and Cox regressions performed after binarising dMRI metrics (higher/lower than the median of the cohort). Panels on the right report results from Cox regressions assessing the continuous dependence of PFS on Diff-in metrics. We detect a statistically significant dependence of the risk of progression on baseline vCSMRI (log-rank test: p = 0.047, Fig. 6C) and CDMRI (log-rank test: p = 0.035, Fig. 6E). These differences correspond to statistically significant Hazard Ratios (HRs) from Cox regression (HR = 0.47, p = 0.050 for binarised vCSMRI; HR = 2.36, p = 0.043 for binarised CDMRI). The risk of progression is about twice as high in patients whose baseline cell size vCSMRI is smaller than the median vCSMRI, or whose cell density CDMRI is higher than the median CDMRI. Importantly, we obtained similar results when vCSMRI and CDMRI were not binarised, but rather used as continuous predictors in Cox regressions (right panels in Fig. 6). In this latter case, the HR was statistically significant for vCSMRI (HR = 0.65, p = 0.034, Fig. 6D), and it approached statistical significance for CDMRI (HR = 1.40, p = 0.055, Fig. 6F). These association are not confounded by age, sex or baseline tumour volume (Supplementary Table S5; HR = 0.59, p = 0.02 for vCSMRI; HR = 1.65, p = 0.01 for CDMRI).
This figure reports on the dependence of patients’ progression-free survival (PFS) on the average value of F, vCS and CD within liver tumours at baseline (i.e., before starting immunotherapy), as obtained by fitting model Diff-in at high b-value. Left: Kaplan-Meier (KM) survival curves of two groups obtained by splitting patients based on baseline F (panel A), vCS (panel C) and CD (panel D) (lower/higher than the sample median). The grey panel reports the p-values of a log-rank sum test comparing the KM curves, and of a Cox regression based on the binarised MRI metric (with the corresponding hazard ratio (HR) estimate and 95% confidence interval). The legend reports the Restricted Mean Survival Time (RMST) and Restricted Standard Deviation of Survival Time (RSDST) for each KM curve. Right: results from univariate Cox regression where the baseline F (panel B), vCS (panel D) and CD (panel F) is a continuous predictor of the survival. The panel shows how changes in baseline F, vCS and CD modulate the survival curve, given the HR estimated for each metric. In the grey box, the p-value and HR (with 95% CI) corresponding to the baseline MRI metric are reported. In all panels, the y-axis shows 1 – p, with p being the probability of progression, while the x-axis shows the time to progression (in days). We performed the response assessment using a sample size of N = 30 (Fig. 1).
Fig. 7 and supplementary Fig. S14, Fig. S15 and Fig. S16 report on response assessment based on on all other dMRI metrics (routine ADC and K in Fig. 7; cytosol diffusivity D0,I for model Diff-in in Fig. S14; vascular fraction fv in Fig. S15; all metrics from model Diff-in-exFast in Fig. S16). While the estimated HRs for metrics vCSMRI, FMRI and CDMRI from dMRI model Diff-in-exFast are in the same direction as those from model Diff-in, their association with the probability of progression is weaker (HR = 1.53, p = 0.05 for CDMRI; Table S5). We do not detect any association between baseline ADC and K and the probability of progression (Fig. 7 and Table S5).
This figure reports on the dependence of patients’ progression-free survival (PFS) on the average value of ADC and K within liver tumours at baseline (i.e., before starting immunotherapy). Left: Kaplan-Meier (KM) survival curves of two groups obtained by splitting patients based on baseline ADC (panel A) and K (panel C) (lower/higher than the sample median). The grey panel reports the p-values of a log-rank sum test comparing the KM curves, and of a Cox regression based on the binarised MRI metric (with the corresponding hazard ratio (HR) estimate and 95% confidence interval). The legend reports the Restricted Mean Survival Time (RMST) and Restricted Standard Deviation of Survival Time (RSDST) for each KM curve. Right: results from univariate Cox regression where the baseline ADC (panel B) and K (panel D) is a continuous predictor of the survival. The panel shows how changes in baseline ADC and K modulate the survival curve, given the HR estimated for each metric. In the grey box, the p-value and HR (with 95% CI) corresponding to the baseline MRI metric are reported. In all panels, the y-axis shows 1 – p, with p being the probability of progression, while the x-axis shows the time to progression (in days). We performed the response assessment using a sample size of N = 30 (Fig. 1).
Discussion
The latest liver dMRI models aim to disentangle intra-cellular and extra-cellular water contributions to the total signal6,19,22,23. This powerful approach enables the estimation of innovative tissue property maps, but its clinical deployment is hampered by the high number of unknown tissue parameters to estimate, which requires impractically long dMRI acquisitions19,45,46. With this challenge in mind, this paper delivers a practical implementation of a two-compartment dMRI signal model, tailored for liver imaging, and truly feasible in the clinic. Through histology-informed model selection, we design a compact dMRI framework consisting of fitting a one-pool model of restricted intra-cellular diffusion to highly DW images. The framework provides cell size and density estimates that correlate with histology, and which enable the stratification of the risk of cancer progression under immunotherapy.
To find the optimal dMRI signal implementation, we analysed co-localised dMRI and histology data (N = 25) from fixed mouse livers and from cancer patients. We compared 5 signal models, each fitted according to two distinct strategies, and ranked them for their ability to estimate intra-cellular fraction and cell size, as seen on histology. Rankings unequivocally suggest the highest radiological-histological agreement is obtained by fitting a single-compartment model of restricted diffusion within spherical cells – a model here referred to as Diff-in –, to images acquired with b-values higher than approximately 900 s/mm2 in vivo and 1800 s/mm2 ex vivo. Interestingly, our central result, confirmed by Monte Carlo computer simulations, points towards the fact that simpler models of diffusion can provide the highest fidelity to the underlying histology, if deployed in appropriate measurement regimes. In practice, our approach suggests that focussing on measurement regimes where the signal is dominated by intra-cellular diffusion, may enable the deployment of simpler models, which still suffice to capture salient microstructural features7,47–49. Notably, our model selection results are consistent with the fact recent estimates of the extra-cellular liver ADC, as high as 2.5 µm2/ms19. Such a high ADCE implies that the extra-cellular signal would decay to roughly 5% or less of its non-DW value even for b-values of 1200 s/mm2 (exp(–b ADCE) ≈ 0.05 for b = 1200 s/mm2 = 1.2 ms-1 µm2 and ADCE = 2.5 µm2 ms-1), justifying the use of single-compartment model of intra-cellular diffusion7,44,50.
Importantly, we conducted dMRI model selection using a variety of criteria, namely: TCS (overall correlation between dMRI and histology), HFC (accuracy in histological property estimation via dMRI), and BIC35 (dMRI model quality of fit, penalising model complexity). In general, all criteria point towards the same direction, with model Diff-in fitted to high b-value images being the top-ranking model. We would also like to emphasise that while we used 15 in vivo dMRI scans for both model design and for the clinical demonstration, the same liver dMRI model implementation (i.e., fitting model Diff-in on high b-value images) would have been selected had we only looked at the ex vivo mouse data or at in silico signals (Fig. S4, Tables S1 to S3). This fact gives confidence on the robustness and generalisability of our model design. Another important observation is that both TCS and HFC, which are histology-informed, suggest that models constraining the extra-cellular ADC to be higher than the intra-cellular ADC, outperform models without such a constraint. Such a better performance does not stand out as clearly when looking at BIC rankings instead (e.g., Fig. S3.A, Fig. S4.A, Fig. S4.B, Fig. S4.C). This minor discrepancy can be understood considering that a good fitting quality may not necessary imply accurate parameter estimation in presence of noise26,44. It also stresses the importance of informing dMRI modelling with histology, for improved biological specificity.
After selecting a practical dMRI model implementation, we investigated its utility in an exemplificatory response assessment task. Immunotherapy has shown promise in several cancers16. However, only a small fraction of patients truly benefits from this class of drugs, and their identification prior to treatment is extremely challenging17. With this in mind, we investigated whether it is possible to stratify the probability of progression under immunotherapy given baseline dMRI metrics within liver tumours (N = 30). Several, independent statistical tests consistently point towards the fact that smaller baseline dMRI cell sizes and higher dMRI cell densities are associated to faster cancer progression (shorter PFS). To our knowledge, this is the first time that such advanced dMRI markers have been tested for patient stratification in immunotherapy in vivo. The deleterious impact of higher cell density on PFS is in line with other studies focussing on different treatments, where higher tumour cellularity has been associated with higher cancer aggressiveness or worse prognosis. In CRC liver metastasis resection, for example, high cell density in resected metastases has been associated to shorter disease-free survival51. Conversely, in breast cancer, lower tumour cellularity has been associated to pathologic complete response in chemo-free dual HER2 blockade treatment (for HER2-positive breast cancer)52, as well as longer survivals in neoadjuvant chemotherapy53. The significant association between PFS and dMRI cell size/density reported here is promising, and motivates future studies to confirm it.
We benchmarked the proposed dMRI approach against well-established DKI apparent diffusion and kurtosis coefficients (ADC and K)30. ADC and K are easy to compute with compact dMRI acquisitions, and are sensitive to cancer cellularity46, a fact confirmed by our data. However, in our cohort, neither of ADC and K showed statistically significant associations with clinical PFS. This may due, at least partly, to the fact that ADC and K are semi-quantitative, protocol-dependent metrics, whose value can change as function of factors as the diffusion time46,54. Here, we did not perform inter-scanner harmonisation deliberately, to stress the quantitative nature of our cell size/density mapping approach, which inherently accounts for inter-scanner protocol differences. However, it is possible that better performances for semi-quantitative ADC and K may be obtained by adopting inter-scanner harmonisation55.
We would like to acknowledge the following potential limitations. Firstly, our sample size (N = 25 for model development; N = 30 for the clinical demonstration) is relatively small. This paper provides a first demonstration of the potential utility of the proposed Diff-in approach. The demonstration is unique of its kind, since it reports heretofore undescribed dMRI-based stratification in immunotherapy. Nonetheless, while works proposing related dMRI techniques relied on similar23,36, if not even smaller19,20, sample sizes, we acknowledge that our exploratory findings require further confirmation in larger cohorts.
Secondly, we point out that results from any dMRI-histology comparison should always be taken with care. Here we related dMRI metrics obtained in vivo to histological indices from biopsies. While we were able to identify the tumours from which the biopsies were taken, we could not identify exactly the tumour area that was biopsied. This may imply that the biopsies are not fully representative of the tumour microenvironment in its entirety. Also, and most importantly, histology has its own limitations, since it provides cell property estimates that may not be, per se, fully accurate. For example, routine HE histology is an inherently 2D technique, unlike 3D MRI. Moreover, it is affected by artifacts (e.g., due to dehydration, paraffin embedding, imperfect staining, cutting, etc56), and the automatic processing of large fields-of-view requires trading off between sensitivity and specificity. We took steps to mitigate these issues, e.g., by accounting for biases due to tissue shrinkage. Nonetheless, our histology-derived estimates of cell properties are likely biased versions of the true figures.
We would also like to acknowledge that the proposed dMRI approach neglects other potentially relevant microstructural properties, such as water exchange between intra-/extra-cellular spaces21,37, presence of cell size/cytosolic diffusivity distributions42,57, or intra-compartmental T2 or T136. On the one hand, ignoring these properties may have biased the estimation of FMRI and vCSMRI21,37. On the other hand, properties such as exchange rates, overlooked in our model, may be relevant markers of cellular stress per se. In future, we plan to incorporate these properties in our models, while ensuring the clinical feasibility of the dMRI protocols required to fit them.
To conclude, this study delivers a practical liver dMRI signal model consisting of a single-compartment of restricted diffusion within spherical cells, which should be fitted to b-values higher than, approximately, 900 s/mm2in vivo. This model offers estimates of cell size and cell density that are correlated to the underlying histology, and which may provide complementary information to routine volumetric tumour burden assessment, for example by stratifying the risk of cancer progression in immunotherapy. Striving to bringing precision imaging one step closer to the clinic, we release our approach as an easy-to-use, open-source Python implementation, which will be freely accessible online.
Methods
dMRI models
Common biophysical body dMRI signal models6,19,20,50,58,59 describe the signal as arising from three, non-exchanging proton pools: vascular water; restricted, intra-cellular water; hindered, extra-cellular, extra-vascular water. The dMRI signal for a PGSE measurement at b-value b, gradient duration/separation δ/ Δ, and echo time TE is
Above, s0 is the apparent proton density, fv is the voxel vascular signal fraction, fI is the tissue intra-cellular signal fraction, T2v/T2I/T2E and av/aI/aE are compartment-wise T2 and diffusion-weighting factors.
av captures intra-voxel incoherent motion (IVIM) effects60. In vivo, the IVIM vascular ADC ranges approximately61 [15; 60] µm2 ms-1. For this reason, for b > 100 s/mm2, the vascular signal vanishes (av ≈ 0), and Eq. (1) reduces to19
A common model for aI in Eq. 2 is that of restricted diffusion within spheres of diameter L19,20:
where
is the Gaussian phase distribution approximation of the intra-cellular ADC62. Above, αm is the m-th root of αmR J’3/2(αmR) – 0.5 J3/2(αmR) = 0, J3/2(x) is the Bessel function of the first kind and order 3/2, and J’3/2(x) its first-order derivative. ADCI depends on the intrinsic cytosol diffusivity D0,I and on the cell size L = 2R (R: radius; L: diameter). Noting that dMRI-derived L represents a volume-weighted mean cell size statistics7,42, we will refer to it as volume-weighted cell size (vCS).
Conversely, the extra-cellular, extra-vascular signal may be described in terms of hindered diffusion in a tortuous space19,27,59:
with
In Eq. 6, DE,∞ is the asymptotic59 ADCE for Δ → ∞.
The 5 implementations of the two-compartment model
We investigated 5 implementations of Eq. 2, divided into two families. The first family includes models that do not make assumptions on which of ADCI/ADCE is higher:
Diff-in-exTD: the most general model, relying on the full expression of ADCE in Eq. 6;
Diff-in-ex: a simpler implementation of Diff-in-exTD that neglects extra-cellular TD (β = 0 in Eq. 6).
In the second family of models, we constrain ADCE > ADCI. It includes
Diff-in-exTDFast: equivalent to Diff-in-exTD, with the lower bound for DE,∞ ensuring ADCE > ADCI for any L.
Diff-in-exFast: equivalent to Diff-in-ex, but again ensuring that ADCE > ADCI for any L.
Diff-in: a model where the extra-cellular signal is negligible compared to the intra-cellular one, due to ADCE being much larger than ADCI, so that Eq. 2 simplifies to
In all models we used T2I ≈ T2E ≐ T2T, given the challenge of resolving accurately multiple T2 constants20,36.
Fitting
We fitted the 5 models using custom-written Python routines, based on objective function minimisation initialised by a grid search. The objective function was fobj = −ln(λ), where λ is the offset-Gaussian likelihood34. Fitting provides estimates of vCS and voxel intra-cellular signal fraction
We also combined vCS and F into a cell density per unit volume20
Preclinical data
Animals
We obtained data from 7 fixed livers of NOD.Cg-Prkdcscid IL2rgtm1WjI/SzJ mice. All experimental protocols were approved and monitored by the Vall d’Hebron Institute of Research Animal Experimentation Ethics Committee (CEEA; registration number 68/20) in accordance with relevant local and EU regulations. We studied six livers from mice implanted with cells derived from biopsies of prostate cancer patients, as part of an ongoing study, plus an additional liver from a mouse without any implantation. We implanted one tumor biopsy core with growth factor-enriched Matrigel (Corning) subcutaneously in the flank of each mice. We derived tissue from the following biopsies: iliac bone metastasis biopsy (metastatic castration-resistant prostate cancer, presenting with bone metastasis and Gleason score 3+4 adenocarcinoma); prostate biopsy (patient with metastatic hormone-sensitive prostate cancer, presenting with bone metastasis and Gleason score 5+4 adenocarcinoma); two liver biopsies (patient with metastatic castration-resistant prostate cancer, presenting with bone and visceral metastasis and Gleason score 4+4 acinar adenocarcinoma; patient with metastatic hormone-sensitive prostate cancer, presenting with bone and liver metastasis and Gleason score 4+4 adenocarcinoma). After implantation, we measured tumour size using calipers and monitored mouse weight weekly, sacrificing animals by cervical dislocation under general anesthesia when tumour volume exceeded 2000 mm3. We collected the livers, fixed them overnight in formalin, and transferred them to phosphate-buffered saline (PBS) solution.
MRI
We scanned livers on a 9.4T Bruker Avance system at room temperature. Livers were tightened with sewing thread to a histology cassette and placed into a Falcon® tube, filled with PBS solution. A 1-channel birdcage coil was used (excitation/reception). The protocol included a T2-weighted fast spin echo sequence (resolution: 144 μm × 144 μm × 2.216 mm) and PGSE dMRI (Fig. S17A; TR = 2700 ms; resolution: 386 μm × 386 μm; matrix size: 86 × 86; 4 slices, 2.216 mm-thick, NEX = 1). The protocol featured: δ = 10 ms, Δ = {15, 30} ms, 10 linearly spaced b-values for each Δ (minimum/maximum nominal b: 0/2800 s/mm2). DW images corresponding to Δ = 15 ms were acquired at each of TE = {31, 45, 65} ms, and to Δ = 30 ms at each of TE = {45, 65} ms. We i) denoised dMRI scans with Marchenko-Pastur Principal Component Analysis (MP-PCA) Python denoising63 (kernel: 7×7×3), ii) mitigated Gibbs ringing (MrTrix3 local sub-voxel shift method64), and iii) corrected temporal signal drifts by assessing signal changes in a PBS solution region, accounting for TE (PBS T2: 500 ms).
Finally, we fitted the Diff-in-exTD, Diff-in-exTDFast, Diff-in-ex, Diff-in-exFast and Diff-in models voxel-by voxel (tissue parameter bounds: [0; 1] for fI; [0.8; 2.6] μm2 ms-1 for D0,I; [8; 40] µm for vCS; [0.8; 2.6] μm2 ms-1 for DE,∞ in models Diff-in-ex and Diff-in-exTD and [1.75; 2.6] μm2 ms-1 in models Diff-in-exFast and Diff-in-exTDFast; [0; 10] μm2 for β in models Diff-in-ex-TD and Diff-in-exTDFast). For fitting, we fixed fv and T2T to values obtained through a a two-pool vascular-tissue model65 (fitting bounds: [0; 1] for fv; [5; 80] ms for T2T). Fitting was performed i) on all images with b > 1000 s/mm2 (suppressing vascular signals, referred to as fitting on whole image set); ii) on b > 1000 s/mm2 images (high b-value fitting). In our ex vivo data, the vascular signal captures PBS solution contamination (PBS ADC: roughly 2.4 μm2 ms-1). For this reason, we adopted a b-value threshold of 1000 s/mm2 to achieve acceptable PBS signal suppression. We used instead a minimum b-value of 1800 s/mm2 for high b-value fitting (minimising extra-cellular contributions), given the reduction in intrinsic tissue diffusivity expected ex vivo.
For comparison, we computed ADC and apparent diffusion excess kurtosis K by fitting
to DW images acquired at TE = 45 ms, Δ = 30 ms, with in-house Python code.
Histology
After MRI, samples underwent histology. We cut two 4 μm-thick histological sections at known position, stained them with HE, and digitised them (Hamamatsu C9600-12 slide scanner; 0.227 μm resolution). An experienced pathologist (S.S.) inspected images qualitatively. We then processed them with the automatic cell detection tool of QuPath32, obtaining per-cell area A and diameter 
. Afterwards, we split images into 386 μm × 386 μm patches (matching the MRI resolution), computing patch-wise histological volume-weighted cell size vCShisto, intra-cellular area fraction Fhisto and cell density per unit area CDhisto42. vCShisto, defined as
is a more accurate counterpart of dMRI cell size than the arithmetic mean7,42 aCShisto = < l >. We accounted for biases coming from: i) estimating the size of 3D objects from 2D views (bias 1), ii) tissue shrinkage (bias 2), by rescaling vCShisto and CDhisto. The final vCShisto estimate was 1.4616 times larger than the value obtained from direct image processing (1.4616 = 1.2732×1.148; 1.2732, derived from the theory of spherical caps, accounts for bias 1; 1.148 accounts for bias 2, and corresponds to a plausible shrinkage of 12.9% following dehydration, clearing and paraffin embedding56). The final CDhisto estimate was 1.318 times smaller than the value derived from direct image processing, since 1 mm2 of shrunk tissue corresponds to 1.148×1.148 mm2 = 1.318 mm2 of unprocessed tissue (plausible shrinkage 12.9%56). Lastly, we co-registered histological maps to MRI42 using DiPy66.
Clinical data
Cohort
We obtained data from patients suffering from advanced solid tumours, recruited for an ongoing imaging study approved by the Vall d’Hebron University Hospital Ethics committee (PR(AG)29/2020). Patients, eligible for a phase I immunotherapy trial at VHIO (Barcelona, Spain), provided informed written consent to participate in the imaging study. We included 33 patients with liver malignancies (mean/std of age: 62.91/12.34 year; 16 male, 17 female), of which 3 suffered from primary HCC, while 30 had liver metastases from different primary cancers (10 colon, 8 melanoma, 3 rectal, 2 ovarian, 2 gastric, 2 breast, 1 renal, 1 endometrial, 1 ureteral). We obtained baseline dMRI scans (i.e., acquired immediately before starting immunotherapy), and digitised HE-stained biopsies from one of the imaged liver tumours. We obtained biopsies from 18 patients (6 scanned at 1.5T, 12 at 3T), collected after baseline dMRI. 30 out of 33 patients finally entered the immunotherapy trial after screening; for them, a clinical outcome in the form of PFS was available. PFS represents the lag between therapy starting date and progression or death (whichever occurs first), with progression determined via RECIST13, or in case of established clinical worsening.
MRI
We imaged patients at the level of the abdomen. We scanned 11 patients on a 1.5T Siemens Avanto scanner using the vendor 18-channel body coil for detection. The protocol included a T2-weighted fast spin echo scan (resolution: 1.4 × 1.4 × 5 mm3; 32 slices; TR = 4500 ms; TE = 82 ms; echo train length: 29; NEX = 8; GRAPPA = 2) and fat-suppressed DW TRSE (Fig. S17B) EPI (dMRI scan time: 16 minutes). It featured: resolution: 1.9 × 1.9 × 6 mm3; 32 slices; TR = 7900 ms; bandwidth 1430 Hz/pixel; averaging of 3 orthogonal diffusion directions × 2 signal averages (effective NEX = 6); GRAPPA factor of 2; 6/8 partial Fourier imaging. The dMRI protocol consisted of b = {0, 50, 100, 400, 900, 1200, 1600} s/mm2, each for TE = {93, 105, 120} ms. One additional image (b = 0 s/mm2; TE = 93 ms) was acquired with reversed phase encoding polarity. The gradient timings (Fig. S17B) were: δ1 = 8.9 ms, δ2 = 17.6 ms, δ3 = 20.4 ms, δ4 = 6.0 ms, Δ1,2 = 17.4 ms and Δ1,4 = 63.9 ms when TE = 93 ms; δ1 = 13.2 ms, δ2 = 19.3 ms, δ3 = 24.8 ms, δ4 = 7.7 ms, Δ1,2 = 21.7 ms and Δ1,4 = 74.2 ms when TE = 105 ms; δ1 = 18.9 ms, δ2 = 21.0 ms, δ3 = 30.5 ms, δ4 = 9.5 ms, Δ1,2 = 27.5 ms and Δ1,4 = 87.5 ms when TE = 120 ms. The b-value is
We scanned 22 more patients on a 3T GE SIGNA Pioneer scanner, using the vendor 48-channel torso coil for signal reception, with 32 channels enabled for detection. The protocol included a respiratory-gated T2-weighted fast spin echo scan (resolution: 1.4 × 1.4 × 6 mm3; 32 slices; TR = 4615 ms; TE = 52.86 ms; echo train length: 16) and respiratory-gated, fat-suppressed PGSE (Fig. S10A) EPI (dMRI scan time: 16 minutes). It featured: resolution: 2.4 × 2.4 × 6 mm3; 32 slices; TR = 6000 ms; bandwidth 1953 Hz/pixel; averaging of 3 orthogonal diffusion directions × 2 signal averages (effective NEX = 6); ASSET factor of 2. The dMRI protocol consisted of b = {0, 50, 100, 400, 900, 1200, 1500} s/mm2, each for TE = {75, 90, 105} ms. The gradient timings (Fig. S17A) were: gradient duration δ = {0.0, 3.9, 5.2, 9.2, 15.0, 18.2, 21.0} ms for TE = 75 ms, δ = {0.0, 3.9, 5.2, 9.2, 13.0, 15.8, 18.5} ms for TE = 90 ms and 105 ms; gradient separation Δ = {0.0, 27.8, 29.0, 33.0, 28.7, 31.8, 34.7} ms for TE = 75 ms and Δ = {0.0, 27.8, 29.0, 33.0, 37.0, 39.6, 42.3} ms for TE = 90 ms and TE = 105 ms.
dMRI post-processing consisted of slice-wise Python MP-PCA denoising (kernel: 5 × 5)63; MRTrix3 Gibbs unringing64; motion correction via affine co-registration67; FSL distortion correction68 (1.5T data only). An experienced radiologists (R.P.L.) segmented tumours on the T2-w scan, enabling per-patient tumour volume computation. Afterwards, we warped the tumour mask to dMRI using ANTs69 non-linear co-registration, and fitted the 5 dMRI models, fixing again fv and T2T to previously computed values65 (fitting bounds: [0; 1] for fv; [20; 140] ms for T2T; [0; 1] for fI; [0.8; 3.0] μm2 ms-1 for D0,I; [8; 40] µm for vCS; [0.8; 3.0] μm2 ms-1 for DE,∞ in models Diff-in-ex and Diff-in-exTD, and [1.75; 3.0] μm2 ms-1 in models Diff-in-exFast and Diff-in-exTDFast; [0; 10] μm2 for β in models Diff-in-ex-TD and Diff-in-exTDFast).
We fitted the 5 dMRI models i) on images acquired at a b-value b > 100 s/mm2, to suppress vascular signals (fitting to the whole image set); ii) to b > 900 s/mm2 images, to also minimize extra-cellular contributions (high b-value fitting). For scans performed on the 1.5T Siemens system: i) we used Δ1,2 + δ2 in place of Δ in Eq. 6 (Fig. S10B), ii) we replaced Eq. 4 with a numerical implementation of restricted diffusion within spheres, based on Radial Basis Function interpolation of synthetic signals generated for DW-TRSE with Monte Carlo simulations70.
For both scanners, we also computed ADC and excess kurtosis K by fitting Eq. 10 on b > 100 s/mm2 images (shortest TE), with in-house Python code.
Histology
We performed ultrasound-guided biopsies of one liver tumour at the Barcelona Vall d’Hebron University Hospital (Spain). The biological material underwent standard processing, HE staining and digitalisation (Hamamatsu C9600-12 slide scanner; resolution: 0.454 μm). An experienced pathologist (S.S.) assessed the images and drew a region-of-interest (ROI) outlining the tumours. In parallel, an experienced radiologist (R.P.L.) inspected ultrasound and MR images, outlining the biopsied tumour on the latter. We processed HE data with QuPath and computed per-biopsy vCShisto, Fhisto and CDhisto, as previously described. Finally, we rescaled vCShisto and CDhisto, by multiplying vCShisto from image processing by 1.503 (1.503 = 1.1806×1.2732, where 1.2732 accounts for biases from 2D sectioning, and 1.1806 accounts for a plausible tissue shrinkage of 15.3% following fixation, dehydration, clearing and paraffin embedding56). The final CDhisto estimate was instead 1.3938 times smaller than the value derived from direct image processing, since 1 mm2 of shrunk tissue corresponds to 1.1806×1.1806 mm2 = 1.3938 mm2 of unprocessed tissue for a shrinkage factor of 15.3%56.
Statistical analyses
dMRI model selection
We carried out model selection independently for each of the two fitting strategies. The MRI-histology Total Correlation Score (TCS) selects the model providing the highest Pearson’s correlation between vCSMRI and vCShisto, and between FMRI and Fhisto. It is defined as
where r(vCSMRI, vCShisto) and r(FMRI, Fhisto) are the correlation coefficients of vCSMRI and FMRI with histological vCShisto and Fhisto. The correlation between CDMRI and CDhisto was not included in Eq. 13 since CDMRI is determined analytically from vCSMRI and FMRI. For TCS computation, we pooled together mouse and human data (N = 25).
We also performed model selection using a Histology Fidelity Criterion (HFC), and popular Bayesian Information Criterion (BIC)33,34. HFC rewards the models providing the best accuracy in the numerical estimation of histological cell size and intra-cellular fraction estimation, i.e., minimising
Information on CDMRI and CDhisto was not included in Eq. 14 since CDMRI is not a degree of freedom of the dMRI models (it is determined analytically from vCSMRI and FMRI). BIC selects the model providing the best goodness of fit, penalising complexity, by minimising
λ is the maximised likelihood, and P/N are the number of model parameters/signal measurements. We performed BIC selection voxel-wise, followed by majority voting across voxels.
Simulated dMRI model selection
We synthesised signals via Monte Carlo diffusion random walks for each of the three dMRI protocols of this study, using the MCDC simulator70. We seeded walkers in a substrate made of spherical cells of identical diameter6,19,22,23 (Fig. S5), controlling the intra-sphere fraction F by adding gaps of increasing size in-between abutting spheres, packed in an ideal cubic lattice. We probed four F values (0.197, 0.323, 0.406, 0.523) and four sphere diameters for each F (8, 16, 22 and 30 µm). We varied intra-/extra-sphere diffusivities (10×10 values; [0.8; 2.6] µm2 ms-1 for the ex vivo protocol and [0.8; 3.0] µm2 ms-1 for in vivo protocols), for a total of 1600 synthetic voxels. We corrupted synthetic signals with Rician noise (b = 0 signal-to-noise ratio: 30), and performed model selection according to TCS, HFC and BIC.
dMRI-histology correlation analysis
We computed mean and standard deviation of all metrics i) within the mouse liver samples, ii) within a mask containing all liver tumours in patients, iii) within the biopsied patients’ tumours. We pooled together metrics from mice and patients to calculate Pearson’s correlation coefficients r. In doing so, ADC was normalised to the ADC of the PBS solutions in mice and to the free water diffusivity in patients (3.0 μm2 ms-1 at 37 °C), given the difference in temperature.
Response assessment in immunotherapy
We studied mean values of in vivo baseline dMRI metrics within liver tumours to assess whether these could stratify the probability of progression under immunotherapy, measured by PFS (N = 30). Firstly, we performed the stratification after binarising all MRI metrics as lower/higher than the median of the cohort, creating two groups. We evaluated group-wise survival curves with the Kaplan-Meier estimator, and compared them with a log-rank test. We also fitted a proportional hazard Cox regression, where the binarised MRI metric was the only regressor. Secondly, we stratified PFS without binarising MRI metrics. We fitted a proportional hazard Cox model using each continous dMRI metric as the only regressor (in the form of a z-score), and then fitted the statistical model again, including age, sex and total baseline tumour volume as confounding factors. We performed all analyses in Python, using lifelines.
Data Availability
Raw MRI and histological images from the mouse livers will be released freely online through the Radiomics Group web site following publication in a peer-reviewed journal. Raw MRI and histological images of patients cannot be made freely available due to ethical restrictions at this stage. Python routines enabling the computation of the diffusion MRI metrics presented in this article will be released freely online after publication in a peer-reviewed journal (https://github.com/fragrussu/bodymritools).
Author information
Francesco Grussu and Raquel Perez-Lopez are joint corresponding authors.
Authors and affiliations
Vall d’Hebron Institute of Oncology (VHIO), Vall d’Hebron Barcelona Hospital Campus, Barcelona, Spain. Francesco Grussu, Kinga Bernatowicz, Irene Casanova-Salas, Ignasi Barba, Sara Simonetti, Garazi Serna, Athanasios Grigoriou, Anna Voronova, Rodrigo Toledo, Paolo Nuciforo & Raquel Perez-Lopez.
Cardiff University, Cardiff, United Kingdom. Marco Palombo.
University of Vic - Central University of Catalonia (UVic-UCC), Vic, Spain. Ignasi Barba.
Department of Biomedicine, Faculty of Medicine and Health Sciences, University of Barcelona, Barcelona, Spain. Athanasios Grigoriou & Anna Voronova.
PET/MR Unit, CETIR-ASCIRES, Barcelona, Spain. Valezka Garay.
Department of Radiology, Hospital Universitari Vall d’Hebron, Barcelona, Spain. Juan Francisco Corral, Xavier Merino, Richard Mast, Núria Roson & Manuel Escobar.
Institut de Diagnòstic per la Imatge (IDI), Barcelona, Spain. Juan Francisco Corral, Xavier Merino, Richard Mast, Núria Roson & Manuel Escobar.
Siemens Healthineers, Madrid, Spain. Marta Vidorreta.
GE HealthCare, Madrid, Spain. Pablo García-Polo García.
Medical Oncology Service, Vall d’Hebron Barcelona Hospital Campus, Vall d’Hebron Institute of Oncology (VHIO), Barcelona, Spain. Maria Vieito, Joaquin Mateo & Elena Garralda
Contributions
Conceptualization: F.G., R.P.L., K.B., M.P., E.G., R.T., P.N., J.M. Methodology: F.G., R.P.L., K.B., M.P., A.G. Investigation: F.G., R.P.L., K.B., I.C.S., I.B., S.S., G.S., A.G., V.G., J.F.C., X.M., R.M., N.R., M.E., M.Vie., R.T., P.N., J.M., E.G. Resources: R.P.L., F.G., P.N., J.M., E.G., N.R., M.E., V.G., M.Vid., P.G.P.G, I.B. Formal analysis: F.G. Visualization: F.G. Software: F.G., K.B., A.G. Data curation: F.G., R.P.L., K.B., A.V., G.S., I.C.S., A.G. Project administration: F.G., R.P.L., K.B., E.G., P.N., R.T., J.M., I.C.S. Funding acquisition: R.P.L., E.G., R.T., P.N., J.M., F.G., K.B., I.C.S. Supervision: F.G., R.P.L., E.G., R.T., P.N., J.M. Writing—original draft: F.G., R.P.L., K.B., M.P. Writing—review & editing: all authors
Corresponding authors
Correspondence to Francesco Grussu or Raquel Perez-Lopez.
Ethics declarations
Competing interests
This study received funding from AstraZeneca. M.Vid. works for Siemens Healthineers. P.G.P.G. works for GE HealthCare. K.B. worked as a researcher at the Vall d’Hebron Institute of Oncology (Barcelona), and is now an employee of AstraZeneca. AstraZeneca, Siemens and General Electric did not influence the acquisition and analysis of the data, the interpretation of the results, or the decision to submit the manuscript in its current form for publication.
Ethics
All experimental protocols in animals were approved and monitored by the Vall d’Hebron Institute of Research Animal Experimentation Ethics Committee (CEEA; registration number 68/20) in accordance with relevant local and EU regulations. The imaging study in cancer patients was approved by the Vall d’Hebron University Hospital Ethics committee (PR(AG)29/2020), Barcelona, Spain. Patients provided informed written consent to participate in the study.
Supplementary Information
This article includes Supplementary Figure S1 to S17 and Supplementary Tables S1 to S5.
Supplementary Information
Illustration of the radiological-histological co-localisation on the 7 fixed mouse livers obtained from mice implanted with a biopsy from a prostate cancer patient. (A), left: illustrative slice of the high-resolution anatomical T2-weighted fast spin echo. (B), centre: hematoxylin and eosin (HE)- stained section, taken from the MRI slice shown to the left. (C), right: detail of the microstructure characterising each specimen, as assessed by an experienced pathologist (SS). Different specimens are arranged along different rows. From top to bottom: Control, normal liver structures (no biopsy implantation); PatNA1 and PatNA2, normal appearing normal liver structures after prostate cancer biopsy implantation; Patinf1, Patinf2 and Patinf3: pathology following implantation, consisting of an immature, lympho-proliferative process (infiltration of small cells in sinusoidal spaces); Patnec, pathology following implantation, consisting of necrosis and inflammation.
Values of TCS for all models, as obtained by fitting models on the whole image set (set of images with negligible vascular signal contributions, i.e., b > 1000 s/mm2 in the fixed mouse livers and b > 100 s/mm2 in vivo). We evaluated TCS for histology-informed model selection using a sample size of N = 25.
Frequency of model selection based on the Bayesian Information Criterion (BIC, quantifying how well a model fits the dMRI signals, penalising model complexity) and on the Histology-fidelity Criterion (HFC, quantifying how accurately a dMRI models estimates the intra-cellular fraction and the volume-weighted cell size as seen on histology). (A) reports results when models are fitted to the entire set of measurements with negligible vascular signal constributions (b > 1000 s/mm2 for suppression of PBS fluid within vessels in the fixed mouse livers; b > 100 s/mm2 for IVIM signal suppression in vivo on clinical systems), while (B) reporting results obtained when fitting models only on high b-value images (b > 1800 s/mm2 in the fixed mouse livers; b > 900 s/mm2 in vivo).
Frequency of selection of each of 5 biophysical dMRI models on 3 MRI-histology data subsets. First column: selection on 7 fixed mouse livers scanned ex vivo on a preclinical 9.4T Bruker system (A and D, left); Second column: selection on 6 liver tumours imaged in vivo on a clinical 1.5T Siemens system (B and E, middle); Third column: selection on 12 liver tumours imaged in vivo on a clinical 3T GE system (C and F, right). Plots on top (A to C) refer to dMRI model fitting performed on images where the vascular signal was suppressed (“whole image set fitting”, b > 1000 s/mm2 for suppression of PBS fluid within vessels on the 9.4T; b > 100 s/mm2 for IVIM signal suppression on clinical systems). Plots to the bottom (D to F) refer to dMRI model fitting performed on images where both vascular and extra-cellular, extra-vascular signals were suppressed (“high b-value fitting”, b > 1800 s/mm2 on the 9.4T; b > 900 s/mm2 on clinical systems). Violet: models where extra-cellular ADC is larger than intra-cellular ADC; orange: models with no constraints on which is larger between intra-/extra-cellular ADC. The Bayesian Information Criterion (BIC) selects a model depending on the goodness of MRI signal fitting. The Histology Fidelity Criterion (HFC) selects a model depending on the overall agreement between MRI volume-weighted Cell Size (vCS) and intra-cellular fraction (F) with their histology counterparts.
The synthetic environment consisted of meshed spheres of fixed diameter, representing cells, as this is a common biophysical model used in several dMRI techniques (e.g., VERDICT, IMPULSED). We used the synthetic environment to generate dMRI signals via Monte Carlo simulations for each of the 3 dMRI protocols considered in this study (the PGSE protocol used on the ex vivo mouse livers at 9.4T; the PGSE protocol used in patients in vivo at 3T; the DW TRSE protocol used in patients in vivo at 1.5T). Afterwards, we performed performed dMRI model selection on the synthetic signals, following the same procedures implemented for actual ex vivo and in vivo dMRI data. We controlled the intra-sphere fraction F by adding gaps of increasing size in-between abutting spheres packed in an ideal cubic lattice. We probed 4 different values of F (approximately equal to 0.197, 0.323, 0.406, 0.523; notice that the maximum theoretical value of F for cubic lattice packing is equal to 0.5236). For each value of F, we varied the cell diameter (8, 16, 22 and 30 µm), intra-sphere diffusivity (10 linearly-spaced values in the ranges [0.8; 2.6] µm2 ms-1 and [0.8; 3.0] µm2 ms-1 for the ex vivo and in vivo protocols respectively) and extra-sphere intrinsic diffusivity (again, 10 linearly-spaced values in the ranges [0.8; 2.6] µm2 ms-1 and [0.8; 3.0] µm2 ms-1 for the ex vivo and in vivo protocols respectively), generating a total of 1600 synthetic voxels. Before dMRI signal model fitting, we corrupted synthetic signal with Rician noise at a signal-to-noise ratio (SNR) of 30 on the b = 0 signal s(b=0) (SNR = s(b=0)/σ, where σ is the noise standard deviation).
Matrices illustrating Pearson’s correlation coefficients among all possible pairs of MRI and histology metrics. Histological metrics are: intra-cellular area fraction Fhisto; volume-weighted mean cell size index vCShisto; cell density per unit area CDhisto. MRI metrics are: apparent diffusion coefficient ADC; apparent diffusion excess kurtosis K; intra-cellular area fraction FMRI; volume-weighted mean cell size index vCSMRI; cell density per unit area CDMRI. Metrics FMRI, vCSMRI and CDMRI were obtained by fitting models with no assumptions on which is larger between intra-cellular and extra-cellular ADC (Diff-in-ex and Diff-in-exTD). The 4 panels refer to models Diff-in-ex and Diff-in-exTD fitted according to 2 different strategies. Panel (A): model Diff-in-ex fitted on the whole set of measurements with vascular signal suppression (b > 100 s/mm2 in vivo, b > 1000 s/mm2 ex vivo); panel (B): model Diff-in-ex fitted on high b-value measurements (b > 900 s/mm2 in vivo, b > 1800 s/mm2 ex vivo); panel (C): model Diff-in-exTD fitted on the whole set of measurements with vascular signal suppression; panel (D): model Diff-in-exTD fitted on high b-value measurements. We calculated correlation coefficients using a sample size of N = 25 (Fig. 1).
Matrices illustrating Pearson’s correlation coefficients among all possible pairs of MRI and histology metrics. Histological metrics are: intra-cellular area fraction Fhisto; volume-weighted mean cell size index vCShisto; cell density per unit area CDhisto. MRI metrics are: apparent diffusion coefficient ADC; apparent diffusion excess kurtosis K; intra-cellular area fraction FMRI; volume-weighted mean cell size index vCSMRI; cell density per unit area CDMRI. Metrics FMRI, vCSMRI and CDMRI were obtained by fitting models that assume that the extra-cellular ADC is always larger than the intra-cellular ADC (Diff-in, Diff-in-exFast and Diff-in-exTDFast). The 6 panels refer to models Diff-in, Diff-in-exFast and Diff-in-exTDFast fitted according to 2 different strategies. Panel (A): model Diff-in fitted on the whole set of measurements with vascular signal suppression (b > 100 s/mm2 in vivo, b > 1000 s/mm2 ex vivo); panel (B): model Diff-in fitted on high b-value measurements (b > 900 s/mm2 in vivo, b > 1800 s/mm2 ex vivo); panel (C): model Diff-in-exFast fitted on the whole set of measurements with vascular signal suppression; panel (D): model Diff-in-exFast fitted on high b-value measurements; panel (E): model Diff-in-exTDFast fitted on the whole set of measurements with vascular signal suppression; panel (F): model Diff-in-exTDFast fitted on high b-value measurements. We calculated correlation coefficients using a sample size of N = 25 (Fig. 1).
(A): high-resolution fast spin echo scan acquired in fixed mouse livers scanned ex vivo on the 9.4T Bruker system. (B): standard diffusion metrics, namely ADC (apparent diffusion coefficient) and K (apparent diffusion kurtosis excess). These metrics were obtained by fitting the standard diffusion kurtosis signal representation s = s0 exp(– b ADC + K (b ADC)2/6) to the set of measurements at fixed TE = 45 ms and Δ = 30 ms. From top to bottom, the figure reports maps from 3 specimens, representative of the 3 different microstructural phenotypes seen in our mouse liver data. These are: normal liver structures (illustrated by the Control case, e.g., mouse with no biopsy implantation); pathology following biopsy implantation, consisting of an immature, lympho-proliferative process (infiltration of small cells in sinusoidal spaces, illustrated by case Patinf); pathology following biopsy implantation, consisting of necrosis and inflammation (illustrated by case Patnec).
(A): high-resolution fast spin echo scan acquired in fixed mouse livers scanned ex vivo on the 9.4T Bruker system. (B): metrics from the Diff-in-exFast model fitted to the whole DW image set (b-values with negligible vascular signal contributions, i.e., b > 1000 s/mm2 on fixed ex vivo tissue, to suppress signal from PBS-filled vessels). From left to right: intra-cellular signal fraction FMRI; volume-weighted cell size index vCSMRI index; cell density per unit volume CDMRI. Maps from 3 specimens are reported along different rows. The specimens are representative of the 3 different microstructural phenotypes seen in our mouse liver data. From top to bottom, these are: normal liver structures (illustrated by the Control case, e.g., mouse with no biopsy implantation); pathology following biopsy implantation, consisting of an immature, lympho-proliferative process (infiltration of small cells in sinusoidal spaces, illustrated by case Patinf); pathology following biopsy implantation, consisting of extended necrosis and inflammation (illustrated by case Patnec).
(A): high-resolution fast spin echo scan acquired in fixed mouse livers scanned ex vivo on the 9.4T Bruker system. (B): diffusivity metrics from biophysical model Diff-in-exFast, namely: intrinsic intra-cellular cytosolic diffusivity D0,I; asymptotic extra-cellular diffusion coefficient DE,∞. (C): intrinsic intra-cellular cytosolic diffusivity D0,I from model Diff-in fitted to high b-value images (b > 1800 s/mm2). Maps from 3 specimens are reported along different rows. The specimens are representative of the 3 different microstructural phenotypes seen in our mouse liver data. From top to bottom, these are: normal liver structures (illustrated by the Control case, e.g., mouse with no biopsy implantation); pathology following biopsy implantation, consisting of an immature, lympho-proliferative process (infiltration of small cells in sinusoidal spaces, illustrated by case Patinf); pathology following biopsy implantation, consisting of extended necrosis and inflammation (illustrated by case Patnec).
(A): high-resolution fast spin echo scan as well as a high b-value diffusion image, with biopsied tumour outlined. (B): standard diffusion metrics in the biopsied tumour. Metrics are: apparent diffusion coefficient (ADC) and apparent diffusion kurtosis excess (K). These were obtained by fitting the standard diffusion kurtosis signal representation s = s0 exp(– b ADC + K (b ADC)2/6) to the set of measurements at fixed, minimum TE and b > 100 s/mm2. Maps are shown in four representative patients (two patients for each MRI scanner), along different rows. For the 1.5T Siemens scanner (first and second rows from top): patient 6 (primary hepatocellular carcinoma) and patient 3 (liver metastases from ovarian cancer). For the 3T GE scanner (third and fourth rows from top): patient 24 (primary hepatocellular carcinoma) and patient 30 (liver metastases from breast cancer).
(A): high-resolution fast spin echo scan as well as a high b-value diffusion image, with biopsied tumour outlined. (B): salient metrics of the Diff-in-exFast model fitted to the whole set of images with negligible vascular signal contributions (b > 100 s/mm2). Metrics are shown in the biopsied tumour. From left to right: intra-cellular signal fraction FMRI; volume-weighted cell size index vCSMRI index; cell density per unit volume CDMRI. Metrics are shown in four representative patients (two patients for each MRI scanner), along different rows. For the 1.5T Siemens scanner (first and second rows from top): patient 6 (primary hepatocellular carcinoma) and patient 3 (liver metastases from ovarian cancer). For the 3T GE scanner (first and second rows from bottom): patient 24 (primary hepatocellular carcinoma) and patient 30 (liver metastases from breast cancer).
(A): high-resolution fast spin echo scan as well as a high b-value diffusion image, with biopsied tumour outlined. (B): diffusivity maps from biophysical model Diff-in-exFast in the biopsied tumour. Metrics are: intra-cellular cytosolic diffusivity D0,I and asymptotic extra-cellular diffusion coefficient DE,∞. (C): intra-cellular cytosolic diffusivity D0,I for biophysical model Diff-in fitted being fitted only to high b-value images (b > 900 s/mm2). Metrics are shown in four representative patients (two patients for each MRI scanner), along different rows. For the 1.5T Siemens scanner (first and second rows from top): patient 6 (primary hepatocellular carcinoma) and patient 3 (liver metastases from ovarian cancer). For the 3T GE scanner (first and second rows from bottom): patient 24 (primary hepatocellular carcinoma) and patient 30 (liver metastases from breast cancer).
This figure reports on the dependence of patients’ progression-free survival (PFS) on the average value of the intrinsic intra-cellular cytosol diffusivity D0,I within liver tumours at baseline (i.e., before starting immunotherapy), as obtained by fitting model Diff-in at high b-value. Left (panel A): Kaplan-Meier (KM) survival curves of two groups obtained by splitting patients based on baseline D0,I (lower or higher than the sample median). The grey panel reports the p-values of a log-rank sum test comparing the KM curves, and of a Cox regression based on the binarised MRI metric (with the corresponding hazard ratio (HR) estimate and 95% confidence interval). The legend reports the Restricted Mean Survival Time (RMST) and Restricted Standard Deviation of Survival Time (RSDST) for each KM curve. Right (panel B): results from univariate Cox regression where the baseline D0,I is a continuous predictor of the survival. The panel shows how changes in baseline D0,I modulate the survival curve, given the HR estimated for each metric. In the grey box, the p-value and HR (with 95% CI) corresponding to the baseline MRI metric are reported. In all panels, the y-axis shows 1 – p, with p being the probability of progression, while the x-axis shows the time to progression (in days).
This figure reports on the dependence of patients’ progression-free survival (PFS) on the average value of the vascular signal fraction fv within liver tumours at baseline (i.e., before starting immunotherapy). The vascular signal fraction fv was computed in the initial fitting step, which disentangles the vascular from the non-vascular (tissue) signal, before the latter is split into intra-/extra-cellular contributions in the biophysical model fitting step. The same representation layout as in Fig. S6 was used. Left (A): Kaplan-Meier (KM) analysis, log-rank sum test and Cox regression based on the binarised fv (higher/lower than the sample median). Right (B): Cox regression modelling the probability of survival as a continuous function of baseline fv. In all panels, the y-axis shows 1 – p, with p being the probability of progression, while the x-axis shows the time to progression (in days).
This figure reports on the dependence of patients’ progression-free survival (PFS) on the average value of all Diff-in-exFast metrics within liver tumours at baseline (i.e., before starting immunotherapy). In each row, from top to bottom: PFS based on baseline volume-weighted Cell Size vCS (panels A and B), Cell Density CD (panels C and D), intra-cellular fraction F (panels E and F), intrinsic intra-cellular cytosol diffusivity D0,I (panels G and H), asymptotic extra-cellular diffusion coefficient DE,∞ (panels I and J). The same representation layout as in Fig. S6 was used. Left (A, C, E, G, I): Kaplan-Meier (KM) analysis, log-rank sum test and Cox regression based on the binarised MRI metrics (higher/lower than the sample median). Right (B, D, F, H, J): Cox regression modelling the probability of survival as a continuous function of baseline MRI metrics. In all panels, the y-axis shows 1 – p, with p being the probability of progression, while the x-axis shows the time to progression (in days).
(A): pulsed gradient spin echo (PGSE sequence, also known as Stejskal-Tanner sequence, pulsed-field gradient (PFG), or single linear diffusion encoding) used to acquire data on the 9.4T Bruker system on fixed mouse livers ex vivo and on the 3T GE system on patients in vivo. δ and Δ respectively indicate the diffusion gradient duration and separation, while TE is the echo time. (B): twice-refocussed diffusion-weighted spin echo sequence used to acquire data on the 1.5T Siemens system on patients in vivo. δn and Δn,m respectively indicate the duration of the n-th gradient lobe and the separation time between the n-th and m-th gradient lobes, for n,m = 1, …, 4. TE is again the echo time. In both panels, “Signal readout” corresponds to sampling the center of the k-space (i.e., zero spatial frequency).
We performed model selection on synthetic signals simulated for all the dMRI protocols considered in this study (ex vivo PGSE, used on fixed mouse livers; in vivo PGSE and DW TRSE, used in patients in vivo; see Methods for a full description of the protocols). We fitted the models on protocol subsets obtained with the same b-value thresholds used when analysing actual MRI signals (“Regular fit”: fitting on all b-values with negligible vascular contributions; “High b only fit”: fitting on b-values minimising extra-cellular signal contributions). The Table reports the value of TCS = r(vCSest,vCSgt) + r(Fest,Fgt), where vCS is the cell size, F the intra-cellular fraction, r(x,y) the Pearson’s correlation between variables x and y computed pooling together all synthetic voxels, and where subscripts est and gt respectively indicate estimated and ground truth values. Higher values of TCS point towards better model performance. For each protocol and fitting strategy, the model with the highest TCS is flagged by gray shadowing and bold font.
We performed model selection on synthetic signals simulated for all the dMRI protocols considered in this study (ex vivo PGSE, used on fixed mouse livers; in vivo PGSE and DW TRSE, used in patients in vivo; see Methods for a full description of the protocols). We fitted the models on protocol subsets obtained with the same b-value thresholds used when analysing actual MRI signals (“Regular fit”: fitting on all b-values with negligible vascular contributions; “High b only fit”: fitting on b-values minimising extra-cellular signal contributions). For each model, the table reports the percentage of synthetic voxels where HFC = |vCSest – vCSgt|/vCSgt + |Fest – Fgt|/Fgt was the lowest across all models. Above, vCS is the cell size, F the intra-cellular fraction, and subscripts est and gt respectively indicate estimated and ground truth values. Higher percentages indicate smaller estimation errors, and therefore point towards better model performance. For each protocol and fitting strategy, the model with the highest proportion of synthetic voxels with minimum HFC is flagged by gray shadowing and bold font.
We performed model selection on synthetic signals simulated for all the dMRI protocols considered in this study (ex vivo PGSE, used on fixed mouse livers; in vivo PGSE and DW TRSE, used in patients in vivo; see Methods for a full description of the protocols). We fitted the models on protocol subsets obtained with the same b-value thresholds used when analysing actual MRI signals (“Regular fit”: fitting on all b-values with negligible vascular contributions; “High b only fit”: fitting on b-values minimising extra-cellular signal contributions). For each model, the table reports the percentage of synthetic voxels where the Bayesian Information Criterion (BIC, a standard metric of model fitting quality that penalises model complexity) was the lowest across all models. Higher percentages indicate smaller BIC values across synthetic voxels, and therefore point towards better model fitting quality. For each protocol and fitting strategy, the model with the highest proportion of synthetic voxels with minimum HFC is flagged by gray shadowing and bold font.
The table reports mean and standard deviation (within brackets) of histology and dMRI metrics in the 7 fixed mouse livers that were scanned on a 9.4T Bruker system. Histological maps were computed within patches matching the in-plane MRI resolution and then warped non-linearly to dMRI space. The histological maps are: per-patch intra-cellular area fraction Fhisto, per-patch arithmetic mean cell size aCShisto; per-patch volume-weighted mean cell size vCShisto, cell density per unit patch area CDhisto. dMRI metrics are: apparent diffusion coefficient ADC, apparent diffusion kurtosis excess K, intra-cellular signal fraction FMRI, volume-weighted cell size index vCSMRI, apparent cell density per unit volume CDMRI. Metrics FMRI, vCSMRI and CDMRI are reported for both models Diff-in-exFast and model Diff-in, with Diff-in fitted only to high b-value images (b > 1800 s/mm2). In model Diff-in-exFast, the extra-cellular ADC does not feature diffusion time dependence and is constrained to be larger than the intra-cellular ADC. In model Diff-in, the extra-cellular signal is modelled as negligible compared to the intra-cellular one (i.e., total signal dominated by intra-cellular water). Specimens are: Control (normal liver structures); NA1 and NA2 (normal appearing cases, i.e., normal liver structures despite sub-cutaneous biopsy implantation); Patinf1-3 (cases developing liver pathology following sub-cutaneous biopsy implantation, consisting of small cell infiltration in sinusoidal spaces, in between larger hepatocytes); Patnec (case developing liver pathology following sub-cutaneous biopsy implantation, consisting of necrosis and inflammation). aCShisto, always considerably lower than vCShisto, was included to highlight the impact of the largest cells in the computation of statistics based on weighting by cell volume (vCShisto).
The table reports the Hazard Ratios (HR) for different MRI metrics, with relative 95% confidence interval and p-value, estimated through Cox proportional hazard regressions. The models assessed the dependence of the probability of progression on the baseline mean value of MRI metrics within liver tumours, accounting for sex, age and tumour volume. Results are shown for standard diffusion metrics (apparent diffusion and excess kurtosis coefficients, ADC and K), for the vascular signal fraction fv, for metrics from model Diff-in-exFast (intra-cellular fraction F, volume-weighted cell size vCS, cell density per unit volume CD, intrinsic intra-cellular cytosolic diffusivity D0,I, extra-cellular asymptotic diffusion coefficient DE,∞) and for metrics from model Diff-in (intra-cellular fraction F, volume-weighted cell size vCS, cell density per unit volume CD, intrinsic intra-cellular cytosolic diffusivity D0,I) fitted at high b-value (b > 900 s/mm2). Grey shadowing highlights HRs whose p-value is ≤ 0.05.
Acknowledgments
We thank the whole medical oncology, radiology, pathology, molecular biology, clinical trial, and IT teams at the Vall d’Hebron University Hospital and at the Vall d’Hebron Institute of Oncology in Barcelona (Spain), without whom this study would not have been possible. We are also thankful to the Vall d’Hebron Radiology department and to the ASCIRES CETIR clinical team for their assistance, and to past and present members of the Radiomics group for useful discussion and advice. Finally, we would like to express our sincere gratitude to all patients and their families for dedicating their time to research. VHIO would like to acknowledge: the State Agency for Research (Agencia Estatal de Investigación) for the financial support as a Center of Excellence Severo Ochoa (CEX2020-001024-S/AEI/10.13039/501100011033), the Cellex Foundation for providing research facilities and equipment and the CERCA Programme from the Generalitat de Catalunya for their support on this research. This research has been supported by PREdICT, sponsored by AstraZeneca. This study has been co-funded by the European Regional Development Fund/European Social Fund ‘A way to make Europe’ (to R.P.L.), and by the Comprehensive Program of Cancer Immunotherapy & Immunology (CAIMI), funded by the Banco Bilbao Vizcaya Argentaria Foundation Foundation (FBBVA, grant 89/2017). R.P.L is supported by the “la Caixa” Foundation CaixaResearch Advanced Oncology Research Program, the Prostate Cancer Foundation (18YOUN19), a CRIS Foundation Talent Award (TALENT19-05), the FERO Foundation through the XVIII Fero Fellowship for Oncological Research, the Instituto de Salud Carlos III-Investigación en Salud (PI18/01395 and PI21/01019), the Asociación Española Contra el Cancer (AECC) (PRYCO211023SERR) and the Generitat de Catalunya Agency for Management of University and Research Grants of Catalonia (AGAUR) (2023PROD00178). The project that gave rise to these results received the support of a fellowship from “la Caixa” Foundation (ID 100010434). The fellowship code is “LCF/BQ/PR22/11920010" (funding F.G., A.V., and A.G.) and “LCF/BQ/PI20/11760033” (funding I.C.S). I.C.S. also receives the support of the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 847648. This research has received support from the Beatriu de Pinós Postdoctoral Program from the Secretariat of Universities and Research of the Department of Business and Knowledge of the Government of Catalonia, and the support from the Marie Sklodowska-Curie COFUND program (BP3, contract number 801370; reference 2019 BP 00182) of the H2020 program (to K.B.). M.P. is supported by the UKRI Future Leaders Fellowship MR/T020296/2. A.G. is supported by a Severo Ochoa PhD fellowship (PRE2022-102586).
Footnotes
In this revision, word count has been cut down considerably. We went through all of introduction, methods, results, and discussion carefully, rewriting the text in a more concise way. No changes have been made to the data, analyses, results, and take-home messages - these are exactly the same as the previous version. Other minor changes are related to the formatting of sections such as author contribution statement, etc, as this now better suits a new target journal.
References
