# Quantification of hepatocellular carcinoma heterogeneity with multiparametric magnetic resonance imaging

Stefanie J. Hectors, Mathilde Wagner, Octavia Bane, Cecilia Besa, Sara

Lewis, Romain Remark, Nelson Chen, M. Isabel Fiel, Hongfa Zhu, Sacha

Gnjatic, et al.

► To cite this version:

Stefanie J. Hectors, Mathilde Wagner, Octavia Bane, Cecilia Besa, Sara Lewis, et al., Quantification of hepatocellular carcinoma heterogeneity with multiparametric magnetic resonance imaging. Scientific Reports, 2017, 7, pp.2452. 10.1038/s41598-017-02706-z. http://dx.doi.org/10.1038/s41598-017-02706-z

HAL Id: hal-01542017

https://hal.sorbonne-universite.fr/hal-01542017

Submitted on 19 Jun 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_2.png" /></figure>

OPEN

# Quantification of hepatocellular carcinoma heterogeneity with multiparametric magnetic resonance imaging

Received: 4 January 2017

Accepted: 18 April 2017

Published online: 26 May 2017

Stefanie J. Hectors<fnref n="1" />,², Mathilde Wagner<fnref n="1" />,²,³, Octavia Bane<fnref n="1" />,², Cecilia Besa<fnref n="1" />,², Sara Lewis<fnref n="1" />,², Romain Remark⁴, Nelson Chen<fnref n="1" />,², M. Isabel Fiel⁵, Hongfa Zhu⁵, Sacha Gnjatic ⁶, Miriam Merad⁴, Yujin Hoshida<fnref n="1" /> & Bachir Taouli<fnref n="1" />,²

Tumour heterogeneity poses a significant challenge for treatment stratification. The goals of this study were to quantify heterogeneity in hepatocellular carcinoma (HCC) using multiparametric magnetic resonance imaging (mpMRI), and to report preliminary data correlating quantitative MRI parameters with advanced histopathology and gene expression in a patient subset. Thirty-two HCC patients with 39 HCC lesions underwent mpMRI including diffusion-weighted imaging (DWI), blood-oxygenation-level-dependent (BOLD), tissue-oxygenation-level-dependent (TOLD) and dynamic contrast-enhanced (DCE)-MRI. Histogram characteristics [central tendency (mean, median) and heterogeneity (standard deviation, kurtosis, skewness) MRI parameters] in HCC and liver parenchyma were compared using Wilcoxon signed-rank tests. Histogram data was correlated between MRI methods in all patients and with histopathology and gene expression in 14 patients. HCCs exhibited significantly higher intra-tissue heterogeneity vs. liver with all MRI methods (P < 0.030). Although central tendency parameters showed significant correlations between MRI methods and with each of histopathology and gene expression, heterogeneity parameters exhibited additional complementary correlations between BOLD and DCE-MRI and with histopathologic hypoxia marker HIF1α and gene expression of Wnt target GLUL, pharmacological target FGFR4, stemness markers EPCAM and KRT19 and immune checkpoint PDCD1. Histogram analysis combining central tendency and heterogeneity mpMRI features is promising for non-invasive HCC characterization on the imaging, histologic and genomics levels.

Quantitative functional multiparametric magnetic resonance imaging (mpMRI) of cancer allows for non-invasive assessment of several tumour characteristics, such as cellularity, perfusion and oxygenation, which can be used for tumour characterization and for assessing treatment response<fnref n="1" />. Apparent diffusion coefficient (ADC) quantification using diffusion-weighted imaging (DWI) has shown to correlate with tumour cellularity². Tissue oxygenation can be indirectly assessed using blood oxygenation level-dependent (BOLD) and tissue oxygenation level-dependent (TOLD) MRI³. Tissue perfusion and flow can be quantitatively measured using dynamic contrast-enhanced MRI (DCE-MRI)⁴.

Many studies that employ mpMRI to assess/predict tumour response use central tendency parameters, such as mean or median, over entire regions of interest (ROIs) to determine longitudinal changes in the tumour tissue after treatment⁵. However, such analysis may not represent the exact tumour status, given the intrinsic heterogeneous tumour composition⁵. Heterogeneity analysis of tumour MRI measurements may provide accurate markers

<footnote n="1">¹Translational and Molecular Imaging Institute, Icahn School of Medicine at Mount Sinai, New York, NY, United States. ²Department of Radiology, Icahn School of Medicine at Mount Sinai, New York, NY, United States. ³Sorbonne Universités, UPMC, Department of Radiology, Hôpital Pitié-Salpêtrière, Assistance Publique-Hôpitaux de Paris, Paris, France. ⁴Immunology Institute, Icahn School of Medicine at Mount Sinai, New York, NY, United States. ⁵Department of Pathology, Icahn School of Medicine at Mount Sinai, New York, NY, United States. ⁶Oncological Science, Icahn School of Medicine at Mount Sinai, New York, NY, United States. ⁷Department of Medicine/Division of Liver Diseases, Icahn School of Medicine at Mount Sinai, New York, NY, United States. Correspondence and requests for materials should be addressed to B.T. (email: bachir.taouli@mountsinai.org)</footnote>
of tumour heterogeneity at the genetic, cellular and molecular levels⁵ and thereby allow for a better understanding of tumour characteristics that may affect treatment decisions.

Hepatocellular carcinoma (HCC) lesions are known to exhibit substantial intra- and inter-tumour heterogeneity, due to a large variety in etiological and genetic backgrounds and the long-time development of the disease⁶. Tumour heterogeneity poses a significant challenge for treatment stratification in HCC. While pathological and genetic heterogeneity in HCC lesions have been described⁷, imaging reports on HCC heterogeneity are extremely limited, with only one study reporting visual assessment of HCC heterogeneity on contrast-enhanced MRI⁸, with no study employing quantitative imaging measurements of HCC heterogeneity.

Tumour imaging phenotypes, including histogram features, may potentially correlate with the underlying genotype. Thus, non-invasive imaging, including computed tomography (CT) and MRI, can potentially be used as a surrogate for genomics and transcriptomics (radiogenomics)⁹⁻¹¹. Radiogenomics could potentially be used for the prediction of certain molecular gene signatures¹². A major advantage of imaging versus genomics or histopathological analysis is that imaging captures the entire tumour non-invasively¹². In addition, combined analysis of imaging, histopathology and genomics data may improve characterization of tumours, because next to mutual information, the separate assessments may also provide additional, independent information on tumour properties, which may help predict outcomes¹².

Immunotherapy has recently revolutionized treatment in cancer¹³, ¹⁴. The clinical experience with immunotherapy in HCC is growing, with recent and ongoing phase I trials¹⁵, ¹⁶. The success of such treatment heavily depends on tumour expression of immunotherapy targets, such as immune checkpoints. Identification of imaging features that correlate with gene expression of immunotherapy targets potentially allows for non-invasive prediction of immunotherapy outcome.

The objectives of our study were to: 1) quantify intra- and inter-tumour heterogeneity of HCC lesions using mpMRI, consisting of DCE-MRI, BOLD-MRI, TOLD-MRI and DWI, and 2) report our preliminary results correlating the quantitative MRI parameters with advanced histopathology and gene signatures and expression levels in a subset of patients.

## Patients and Methods

**Patients.** This single-centre prospective study was compliant with the Health Insurance Portability and Accountability Act and approved by the institutional review board at the Icahn School of Medicine at Mount Sinai. Written informed consent was obtained from all subjects. From June 2013 to June 2016, 41 consecutive patients with HCC were enrolled. HCC was diagnosed based on routine imaging by two radiologists in consensus [observer 1 (MW) and observer 2 (CB), with 5 and 6 years of experience in body MRI, respectively], according to the Organ Procurement and Transplantation Network (OPTN) criteria¹⁷. Of the 41 initial patients, 9 patients were excluded because of severe motion in the DCE-MRI acquisition (n = 5), absence of contrast injection (n = 2) or low quality of BOLD-MRI data (n = 2). Thus, 32 patients were finally included (M/F 26/6, mean age 59 y, range 30–71 y), of which 31 had chronic liver disease [chronic hepatitis C (n = 18), chronic hepatitis B (n = 9), non-alcoholic steatohepatitis (n = 2), cryptogenic cirrhosis (n = 1) and alcoholic cirrhosis (n = 1)]. For one patient without history of liver disease, HCC was confirmed at histopathology. Three patients underwent locoregional therapy with either transarterial chemoembolization (TACE) combined with radiofrequency ablation (n = 1), TACE only (n = 1) or radioembolization (n = 1) before MRI (range 111–178 days). Data of 23 of the included patients has been published in at least one of our previous publications³, ¹⁸, ¹⁹. The aim of the previous studies was mainly to quantify each of the individual MRI methods in HCC lesions and to assess test-retest repeatability. Assessment of tumour heterogeneity was beyond the scope of those studies.

**MRI acquisition.** The MRI acquisition was performed on 1.5 T (Siemens Aera, Siemens Healthineers, Erlangen, Germany; n = 19) or 3.0 T [Siemens Skyra (n = 5) or Siemens BioGraph mMR (n = 8)] systems. Each of the systems is equipped with a 32-channel spine and flexible body array coil for signal reception. Subjects were asked to fast for 6 hours to eliminate post-prandial effects on portal venous blood flow²⁰. In addition to axial DWI, BOLD-MRI, TOLD-MRI and DCE-MRI, the MRI acquisition consisted of axial and coronal T₂-weighted imaging, axial dual-echo chemical shift imaging, T₁-weighted imaging before and at a delayed phase after contrast injection. The acquisition parameters for DWI, BOLD-MRI, TOLD-MRI and DCE-MRI are listed in Table 1. For BOLD-MRI and TOLD-MRI, R₂\* and R₁ acquisitions were performed before and at the end of a respiratory oxygen challenge of 10–15 minutes at a flow rate of 15 l/min³. DCE-MRI acquisition consisted of 100 dynamic acquisitions. Eight seconds after the start of the acquisition, a half dose (0.05 mmol/kg) of gadobenate dimeglumine (Multihance, Bracco Diagnostics Inc) followed by a 25 ml saline flush was administered intravenously at a rate of 3 ml/s¹⁹. Half dose of the contrast agent was used to avoid saturation of the signal intensity in the DCE-MRI acquisition, because of the high T₁ relaxivity of the contrast agent.

**Image analysis.** Image analysis was performed on MATLAB platform (version R2015a, MathWorks, Natick, MA, USA), unless specified otherwise. The analysis was performed by two observers in consensus [observer 1 and observer 3 (SH), an MRI physicist with 6 years of experience].

**DWI analysis.** DWI acquisition was performed with 16 b-values for the purpose of the estimation of intravoxel incoherent motion (IVIM) diffusion parameters²¹, ²². However, for this particular study, we only used 3 b-values (60, 400 and 800 s/mm²) to determine the apparent diffusion coefficient (ADC), which was estimated as the slope of a linear fit of the logarithmic signal intensity (SI) curves at these particular b-values.

**BOLD/TOLD analysis.** R₂\* maps were generated by fitting of a mono-exponential decay model to the pixel signal values at all acquired echoes³. T₁ maps were calculated using the Look-Locker recovery equation²³ or the Ernst
<table><caption>**Table 1.** Acquisition parameters of the multiparametric MRI protocol. EPI = echo planar imaging, FA = flip angle, FLASH = fast low angle shot, FOV = field-of-view, IR-LL = inversion recovery Look-Locker, GRE = multi-echo gradient echo, TE = echo time, TI = inversion time, TR = repetition time, VFA = variable flip angle.</caption><thead><tr><th></th><th>DWI</th><th colspan="4">BOLD</th><th colspan="3">TOLD</th><th>DCE-MRI</th></tr></thead><tbody><tr><td>Sequence type</td><td>Single-shot spin-echo EPI</td><td>2D GRE</td><td>2D GRE</td><td>2D GRE</td><td>2D GRE</td><td>3D VFA</td><td>2D IR-LL</td><td>2D-IR-LL</td><td>3D FLASH</td></tr><tr><td>Acquisition plane</td><td>Axial</td><td>Axial</td><td>Axial</td><td>Axial</td><td>Axial</td><td>Axial</td><td>Axial</td><td>Axial</td><td>Axial</td></tr><tr><td>MR system</td><td>All</td><td>Skyra (n = 5)</td><td>Aera (n = 10)</td><td>Aera (n = 9)</td><td>mMR (n = 8)</td><td>Skyra (n = 5)</td><td>Aera (n = 19)</td><td>mMR (n = 8)</td><td>All</td></tr><tr><td>TE (ms)</td><td>74–81</td><td>2.5, 4.9, 7.4, 9.8, 12.3, 14.8, 17.2</td><td>4.8, 9.5, 14.3, 19.1, 23.9</td><td>1.7, 2.9, 4.3, 5.6, 6.9, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0</td><td>1.1, 2.4, 3.8, 5.2, 6.6, 8.0, 10.0, 12.0, 15.0, 20.0, 25.0, 30.0</td><td>1.2</td><td>1.0</td><td>1.2</td><td>1.0</td></tr><tr><td>TR (ms)</td><td>4500 (Skyra); one respiration (Aera/mMR);</td><td>165</td><td>242</td><td>319</td><td>249</td><td>10</td><td>2.3</td><td>35.1</td><td>2.7</td></tr><tr><td>FA (°)</td><td>90</td><td>35</td><td>35</td><td>20</td><td>18</td><td>1, 10, 19</td><td>8</td><td>10</td><td>9.5–11.5</td></tr><tr><td>TI (ms)</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>42–1577 (32 TI's)</td><td>80–1445 (40 TI's)</td><td>—</td></tr><tr><td>b-values (s/mm²)</td><td>0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 175, 200, 400, 600, 800</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>FOV (mm²)</td><td>380 × 255</td><td>340 × 255</td><td>340 × 255</td><td>340 × 255</td><td>340 × 255</td><td>350 × 260</td><td>420 × 290</td><td>380 × 280</td><td>380 × 280</td></tr><tr><td>Matrix</td><td>320 × 240</td><td>512 × 384</td><td>384 × 288</td><td>384 × 288</td><td>384 × 288</td><td>384 × 288</td><td>128 × 88</td><td>384 × 288</td><td>384 × 288</td></tr><tr><td>Slice thickness (mm)</td><td>7–8</td><td>7</td><td>7</td><td>7</td><td>7</td><td>5</td><td>8</td><td>8</td><td>4–5</td></tr><tr><td>Number of slices</td><td>20–30</td><td>4–5</td><td>4–5</td><td>4–5</td><td>4–5</td><td>36</td><td>1–2</td><td>1</td><td>40</td></tr><tr><td>Acceleration factor</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>—</td><td>—</td><td>4</td></tr><tr><td>Acquisition time (min:s)</td><td>8:00</td><td>0:15</td><td>0:15</td><td>0:15</td><td>0:15</td><td>0:14</td><td>0:18</td><td>0:10</td><td>0:02 per dynamic</td></tr></tbody></table>

equation for a spoiled gradient recalled echo (SPGR) sequence for the variable flip angle measurements²⁴. The pixel $T_1$ values were inverted to obtain $R_1$ values. In addition, $\Delta R_1$ and $\Delta R_2^*$ maps were generated by subtraction of the baseline $R_1$ and $R_2^*$ measurements from the measurements after oxygen ($\Delta R_{1/2^*} = R_{1/2^*} \text{ post } O_2 - R_{1/2^*} \text{ pre } O_2$)³.

**DCE-MRI analysis.** Prior to pharmacokinetic analysis, motion correction was applied on the DCE-MR images with 3D rigid registration using open-source image analysis software (Firevoxel, CAI²R, New York University, New York, NY, USA). Free-hand 3D regions-of-interest (ROIs) were drawn in the portal vein on the registered images and in the abdominal aorta at the level of the celiac trunk on the unregistered images. Pixel dynamic SI curves were converted to dynamic contrast agent concentration ([CA]) curves by using the signal equation for a SPGR sequence, pre-contrast $R_1$ values as determined from the baseline TOLD measurements (registered to the DCE-MRI measurements, see below) and the contrast agent's relaxivity²⁵. For the portal vein and aorta, pre-contrast $R_1$ values were based on literature²⁶. A haematocrit value of 0.45 was used for conversion from blood [CA] to plasma [CA]. Contrast agent relaxivity values of 8.1 $mM^{-1} s^{-1}$ and 6.3 $mM^{-1} s^{-1}$ were used at 1.5 T and 3.0 T, respectively²⁷.

Pharmacokinetic parameter maps were generated by fitting of the dual-input signal compartment model to the dynamic [CA] curves in each pixel²⁸. The following parameters were estimated in tumours and liver parenchyma: arterial flow ($F_a$), portal flow ($F_p$), total flow ($F_t = F_a + F_p$), arterial fraction (ART = $F_a/F_p$), mean transit time (MTT) and the distribution volume fraction of the contrast agent through the tissue (DV).

**Coregistration of mpMRI data.** Although the orientation of all sequences was in the axial plane, there were acquisition differences related to different matrix size, field of view, offset and subject motion. Therefore, a coregistration algorithm (MATLAB function fitgeotrans) was used to register the mpMRI data. Landmark points were placed on both the raw DCE-MR images and the corresponding BOLD, TOLD or ADC parameter maps. Subsequently, nonrigid coregistration was performed, visually checked and repeated if necessary. The coregistration was only performed in slices that contained HCC lesion(s) and that were common to all acquisitions.

**Region of interest (ROI) analysis.** ROIs in the HCC lesions and liver parenchyma were placed by observer 1 on the registered DCE-MR images. For the HCC ROIs, the entire lesion was included, both viable and, if present, necrotic regions. Subsequently, the ROIs were copied and pasted on the BOLD, TOLD and DWI parameter maps and adjusted if necessary to avoid inclusion of artefacts. Histogram quantification included determination of central tendency (mean and median) and heterogeneity [standard deviation (SD), kurtosis and skewness] parameters. Kurtosis determines the peakedness of the distribution. A higher kurtosis has been associated with increased heterogeneity²⁹, because distributions with high peaks often have fat tails of outlier values. A skewness
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Histopathology processing of HCC sections obtained from a 70-year-old male patient with NASH cirrhosis and HCC. Zoomed microscopy images (×20) of an HCC section stained using Multiplexed Immunohistochemical Consecutive Staining on Single Slide (MICSSS) for CD3, CD31 and CD68 and a separate section stained for HIF1α are shown in the top row. Results of the automatic segmentation algorithm are shown in the bottom row, in which pixels identified as stained are coloured in red. CD3 = cluster of differentiation 3, CD31 = cluster of differentiation 31, CD68 = cluster of differentiation 68, HIF1α = hypoxia-inducible factor 1-alpha.</figcaption></figure>

value deviating from zero indicates that there is tail toward either the left or right side of the distribution, which is also indicative of heterogeneity.

**Histopathology.** Fifteen patients underwent partial hepatectomy after the MRI examination (time range between MRI and resection 0–40 days) and two patients had a biopsy of their index HCC lesion (57 days before and 15 days after the MRI exam). Histological grade³⁰ was available for the HCC lesions in these 17 patients. For one patient that underwent resection, no sections were available since the resection was performed in an outside hospital. For each HCC lesion in the 14 remaining resected liver samples, two five-µm-thick formalin-fixed paraffin-embedded sections were available for advanced histopathological staining. The Multiplexed Immunohistochemical Consecutive Staining on Single Slide (MICSSS)³¹ technique was used to sequentially stain endothelial cells (CD31 monoclonal antibody, clone JC70), macrophages (CD68, clone KP1) and T-cells (CD3, clone 2GV6). For detection of hypoxia, a separate HCC section stained for hypoxia-inducible factor 1-alpha (HIF1α) activity was used for each lesion³². The stained HCC sections were digitally imaged at 20x magnification using a whole-slide scanner (Olympus, Tokyo, Japan). A threshold-based segmentation method was implemented in MATLAB to detect stained pixels (Fig. 1). For the segmentation, the red, green and blue intensity values from the images were retrieved and scaled from 0 to 255. Since the stained pixels generally had a dark red appearance, a pixel was considered stained if (1) the red intensity was higher than 1.1 times the green channel and (2) the average signal intensity was lower than 120. These thresholds were empirically chosen based on inspection of all histological images. From these segmentations, stained tumour fractions were derived.

**Gene expression analysis.** Histologically distinct components of the 14 resected HCC tumours were macro-dissected using H&E staining of serial tissue sections as reference as previously described³³. Isolated total RNA samples were profiled to determine transcriptomic HCC subtypes S1, S2 and S3³⁴ with digital transcript counting technology using the nearest template prediction algorithm as previously described³³. In addition, the following HCC marker genes were simultaneously profiled: liver specific Wnt target overexpressed in S3 subtype (GLUL)³⁵, stemness markers (EPCAM, KRT19)³⁶,³⁷, early HCC markers (BIRC5, HSP70, LYVE1, EZH2)³⁸,³⁹, a pharmacological target currently under clinical testing (FGFR4)⁴⁰, potentially targetable angiogenesis marker (VEGFA)⁴¹ and targetable immune checkpoints (CfFigD274, PDCD1, CTLA4)⁴².

**Statistical analysis.** All values are presented as mean ± SD. For all correlation analyses, the Spearman correlation test was used. For all tests, a P-value lower than 0.05 was considered significant. For the analyses of differences and correlations between techniques, TOLD parameters, except ΔR₁, were not compared with the DCE-MRI parameters, because the baseline R₁ measurements were used as input for the DCE-MRI modelling.

**Intra-tumour heterogeneity.** Differences in kurtosis and skewness values between parameters of different MRI techniques were tested for significance using a Friedman test followed by post-hoc Bonferroni-corrected Wilcoxon signed rank tests. Mean, median and SD values were not compared between the techniques, because they directly depend on the actual MRI parameter values, which are inherently different between the MRI techniques. Liver and HCC histogram parameter values were compared using Wilcoxon signed-rank tests. For
<table><caption>**Table 2.** Central tendency (mean and median) and heterogeneity (SD, kurtosis, skewness) MRI parameters in 39 HCC lesions and liver parenchyma in 32 patients. *P*-values from Wilcoxon signed rank tests. Significant *P*-values are bolded. ADC = apparent diffusion coefficient (10⁻³ mm²/s), ART = arterial fraction (%), DV = distribution volume (%), F<sub>a</sub> = arterial flow (ml/100 g/min), F<sub>p</sub> = portal flow (ml/100 g/min), F<sub>t</sub> = total flow (ml/100 g/min), MTT = mean transit time (s), R<sub>1</sub> = longitudinal relaxation rate (s⁻¹), R<sub>2</sub>* = transverse relaxation rate (s⁻¹).</caption><thead><tr><th rowspan="2"></th><th colspan="3">Mean</th><th colspan="3">Median</th><th colspan="3">SD</th><th colspan="3">Kurtosis</th><th colspan="3">Skewness</th></tr><tr><th>Liver</th><th>HCC</th><th>P*</th><th>Liver</th><th>HCC</th><th>P*</th><th>Liver</th><th>HCC</th><th>P*</th><th>Liver</th><th>HCC</th><th>P*</th><th>Liver</th><th>HCC</th><th>P*</th></tr></thead><tbody><tr><td>F<sub>a</sub></td><td>49 ± 45</td><td>247 ± 240</td><td><strong>&lt;0.001</strong></td><td>47 ± 45</td><td>241 ± 250</td><td><strong>&lt;0.001</strong></td><td>16.3 ± 13.0</td><td>83.0 ± 66.7</td><td><strong>&lt;0.001</strong></td><td>4.3 ± 1.8</td><td>5.5 ± 6.0</td><td>0.754</td><td>0.76 ± 0.66</td><td>0.44 ± 1.34</td><td>0.232</td></tr><tr><td>F<sub>p</sub></td><td>145 ± 123</td><td>112 ± 156</td><td><strong>0.048</strong></td><td>146 ± 126</td><td>103 ± 166</td><td><strong>0.010</strong></td><td>39.9 ± 32.5</td><td>78.6 ± 78.6</td><td><strong>0.020</strong></td><td>6.3 ± 14.0</td><td>64.1 ± 183.6</td><td><strong>0.005</strong></td><td>0.54 ± 1.57</td><td>4.00 ± 6.41</td><td><strong>&lt;0.001</strong></td></tr><tr><td>F<sub>t</sub></td><td>193 ± 122</td><td>358 ± 337</td><td><strong>0.023</strong></td><td>193 ± 125</td><td>342 ± 350</td><td>0.092</td><td>43.7 ± 36.2</td><td>144.3 ± 126.3</td><td><strong>&lt;0.001</strong></td><td>4.5 ± 3.1</td><td>17.1 ± 54.9</td><td>0.969</td><td>0.62 ± 0.76</td><td>1.14 ± 3.11</td><td>0.984</td></tr><tr><td>ART</td><td>33.9 ± 27.9</td><td>78.9 ± 21.5</td><td><strong>&lt;0.001</strong></td><td>33.0 ± 29.1</td><td>80.3 ± 24.2</td><td><strong>&lt;0.001</strong></td><td>9.7 ± 5.1</td><td>13.1 ± 9.0</td><td>0.055</td><td>7.6 ± 10.5</td><td>30.1 ± 61.8</td><td>0.256</td><td>0.68 ± 1.71</td><td>-2.49 ± 4.37</td><td><strong>&lt;0.001</strong></td></tr><tr><td>MTT</td><td>20.3 ± 12.3</td><td>21.2 ± 15.8</td><td>0.724</td><td>19.7 ± 12.1</td><td>18.8 ± 16.3</td><td>0.667</td><td>4.2 ± 2.9</td><td>10.9 ± 7.3</td><td><strong>&lt;0.001</strong></td><td>5.2 ± 4.8</td><td>9.1 ± 17.4</td><td>0.938</td><td>0.90 ± 0.87</td><td>1.54 ± 2.06</td><td>0.337</td></tr><tr><td>DV</td><td>37.4 ± 22.7</td><td>34.6 ± 24.7</td><td>0.695</td><td>36.8 ± 22.8</td><td>34.0 ± 26.0</td><td>0.710</td><td>4.9 ± 3.4</td><td>9.9 ± 5.7</td><td><strong>&lt;0.001</strong></td><td>3.7 ± 1.7</td><td>5.5 ± 10.8</td><td>0.240</td><td>0.58 ± 0.61</td><td>0.70 ± 1.59</td><td>0.468</td></tr><tr><td>R<sub>2</sub>* pre O<sub>2</sub></td><td>55.3 ± 34.4</td><td>34.5 ± 15.1</td><td><strong>&lt;0.001</strong></td><td>54.7 ± 34.8</td><td>33.1 ± 15.0</td><td><strong>&lt;0.001</strong></td><td>7.6 ± 7.5</td><td>7.6 ± 5.2</td><td>0.681</td><td>4.0 ± 2.9</td><td>5.2 ± 4.3</td><td>0.433</td><td>0.35 ± 0.89</td><td>0.92 ± 1.09</td><td>0.055</td></tr><tr><td>R<sub>2</sub>* post O<sub>2</sub></td><td>51.8 ± 30.8</td><td>33.7 ± 14.5</td><td><strong>&lt;0.001</strong></td><td>51.5 ± 31.4</td><td>32.7 ± 14.6</td><td><strong>&lt;0.001</strong></td><td>8.3 ± 10.7</td><td>8.0 ± 4.7</td><td>0.347</td><td>3.5 ± 1.8</td><td>4.0 ± 2.7</td><td>0.984</td><td>0.15 ± 0.64</td><td>0.72 ± 0.81</td><td><strong>0.011</strong></td></tr><tr><td>ΔR<sub>2</sub>*</td><td>-3.47 ± 10.0</td><td>-0.78 ± 8.74</td><td>0.281</td><td>-3.12 ± 8.8</td><td>-0.86 ± 7.48</td><td>0.264</td><td>8.8 ± 9.5</td><td>8.6 ± 5.2</td><td>0.638</td><td>3.9 ± 2.0</td><td>4.7 ± 3.6</td><td>0.830</td><td>-0.18 ± 0.72</td><td>-0.16 ± 0.94</td><td>0.799</td></tr><tr><td>R<sub>1</sub> pre O<sub>2</sub></td><td>1.76 ± 0.47</td><td>1.67 ± 0.94</td><td>0.052</td><td>1.76 ± 0.46</td><td>1.58 ± 0.73</td><td><strong>0.033</strong></td><td>0.12 ± 0.14</td><td>0.75 ± 3.31</td><td><strong>0.008</strong></td><td>3.7 ± 2.2</td><td>8.1 ± 17.7</td><td>0.531</td><td>-0.12 ± 0.76</td><td>0.72 ± 2.07</td><td><strong>0.012</strong></td></tr><tr><td>R<sub>1</sub> post O<sub>2</sub></td><td>1.87 ± 0.73</td><td>1.90 ± 1.79</td><td><strong>0.027</strong></td><td>1.77 ± 0.63</td><td>1.54 ± 0.61</td><td><strong>0.026</strong></td><td>2.26 ± 7.73</td><td>1.44 ± 5.48</td><td>0.422</td><td>71.4 ± 257.6</td><td>16.5 ± 34.5</td><td>0.327</td><td>2.13 ± 8.07</td><td>1.41 ± 2.76</td><td>0.153</td></tr><tr><td>ΔR<sub>1</sub></td><td>0.11 ± 0.45</td><td>0.23 ± 1.15</td><td>0.176</td><td>0.01 ± 0.24</td><td>-0.05 ± 0.25</td><td>0.829</td><td>2.26 ± 7.73</td><td>1.06 ± 3.46</td><td>0.196</td><td>72.2 ± 257.4</td><td>17.6 ± 44.6</td><td>0.814</td><td>2.14 ± 8.08</td><td>1.50 ± 3.03</td><td>0.112</td></tr><tr><td>ADC</td><td>1.31 ± 0.58</td><td>1.43 ± 0.68</td><td><strong>&lt;0.001</strong></td><td>1.31 ± 0.58</td><td>1.41 ± 0.70</td><td><strong>&lt;0.001</strong></td><td>0.2 ± 0.1</td><td>0.3 ± 0.2</td><td>0.142</td><td>3.8 ± 2.2</td><td>4.0 ± 6.4</td><td><strong>0.030</strong></td><td>0.10 ± 0.79</td><td>0.33 ± 0.90</td><td>0.337</td></tr></tbody></table>

patients with multiple HCC lesions, the histogram features of the different lesions were averaged for comparison with liver parenchyma.

**Inter-tumour heterogeneity.** The absolute coefficient of variation (CV = 100*SD/absolute mean) between all MRI parameter features across lesions was measured. Absolute CV was calculated because of the presence of negative skewness values.

**Correlation between MRI methods.** The correlation of the mean, median, SD, kurtosis and skewness measures between the parameters of different MRI methods was assessed. The correlation analysis between mean and median parameters was performed both on all lesions as well as separately for lesions measured at 1.5 T and 3.0 T. To assess the effect of field strength on the correlations, differences between the correlations at both field strengths were tested for significance using a Fisher transformation test.

**Correlation with histopathology and gene expression levels.** The MRI histogram features in the tumour ROIs were correlated with histopathology stained tumour fractions and gene expression levels. In addition, differences in MRI parameters between HCC grades and molecular subclasses were assessed for statistical significance using Kruskal-Wallis tests.

## Results

Forty-four HCCs were identified in 32 patients. Five lesions were excluded because they were not covered by TOLD/BOLD measurements (n = 4) or were completely necrotic due to previous treatment (n = 1). Two other treated lesions were included since they exhibited limited necrosis (<30%). The final 39 included lesions had an average size of 4.4 ± 3.3 cm (range 1.5–13 cm). The number of analysed lesions per patient was as follows: 1 (n = 27, 84%), 2 (n = 4, 13%) or 4 (n = 1, 3%). The distribution of histological grade for the lesions with pathologic validation was as follows: well differentiated (n = 2), moderately differentiated (n = 9), poorly differentiated (n = 5) and undifferentiated (n = 1).

**Intra-tumour heterogeneity.** Average mean, median, SD, kurtosis and skewness values for all assessed parameters in the liver parenchyma and HCC lesions are listed in Table 2. HCC lesions exhibited more intra-tissue heterogeneity compared to liver parenchyma (Fig. 2), with significantly higher SD (F<sub>a</sub>, F<sub>p</sub>, F<sub>t</sub>, MTT, DV and R<sub>1</sub> pre O<sub>2</sub>; P < 0.020), kurtosis (F<sub>p</sub> and ADC; P < 0.030) and skewness (F<sub>p</sub>, ART, R<sub>2</sub>* post O<sub>2</sub> and R<sub>1</sub> pre O<sub>2</sub>; P < 0.012). Kurtosis and skewness in HCC lesions were significantly different between MRI parameters of different techniques (P < 0.001). Specifically, significantly higher kurtosis was observed for F<sub>p</sub>, R<sub>1</sub> post O<sub>2</sub> and ΔR<sub>1</sub> compared to ADC (P < 0.045). A significantly higher skewness was observed for F<sub>p</sub> vs. ΔR<sub>2</sub>* (P = 0.020), while ART showed a significantly lower skewness compared to R<sub>2</sub>* post O<sub>2</sub>, ΔR<sub>3</sub>*, ΔR<sub>1</sub> and ADC (P < 0.030). Leptokurtosis, i.e. a kurtosis value higher than 3, which is indicative of intra-tumour heterogeneity, was observed in all lesions for at least one MRI parameter. In 9 lesions (23%) leptokurtosis was seen in all MRI methods (DCE-MRI, BOLD, TOLD and DWI), while the remaining lesions showed leptokurtosis in a subset of MRI methods [DCE-MRI, BOLD
a

<figure><img src="image_4.png" /></figure>

b

<figure><img src="image_5.png" /><figcaption>**Figure 2.** 54 year-old male patient with cirrhosis secondary to chronic hepatitis B virus infection and HCC. (a) Representative magnified parametric maps of a large (8.3 cm) HCC. Location of the tumour within the liver is indicated by the white arrow on the T₂-weighted image (bottom row, right). A distinct region in the anterior portion of the tumour of high arterial flow (F<sub>a</sub>) and low R<sub>2</sub>* was observed, reflective of high tumour perfusion and normoxia (grey arrow in F<sub>a</sub> and R<sub>2</sub>* pre O<sub>2</sub> maps). The posterior portion of the tumour displays low F<sub>a</sub> and high R<sub>2</sub>*, suggestive of poor perfusion and hypoxia (white arrow in F<sub>a</sub> and R<sub>2</sub>* pre O<sub>2</sub> maps). (b) Histograms of F<sub>a</sub>, R<sub>2</sub>* pre O<sub>2</sub>, R<sub>1</sub> pre O<sub>2</sub> and ADC in the same lesion. The extensive heterogeneity observed in the parameter maps of F<sub>a</sub> and R<sub>2</sub>* pre O<sub>2</sub> is also reflected in the histograms, as illustrated by the fat tails and pronounced skewness, indicated by arrows. The R<sub>1</sub> pre O<sub>2</sub> histogram also exhibited skewness (black arrow). ADC = apparent diffusion coefficient, ART = arterial fraction, DV = distribution volume, F<sub>a</sub> = arterial flow, F<sub>p</sub> = portal flow, F<sub>t</sub> = total flow, MTT = mean transit time, R<sub>1</sub> = longitudinal relaxation rate, R<sub>2</sub>* = transverse relaxation rate.</figcaption></figure>

and TOLD in 18 lesions (46%); DCE-MRI, TOLD and DWI in 4 lesions (10%); DCE-MRI and TOLD in 4 lesions (10%); BOLD and TOLD in 1 lesion (3%); TOLD and DWI in 1 lesion (3%); DCE-MRI only in 2 lesions (5%).

**Inter-tumour heterogeneity.** CV values of MRI parameters across tumours are shown in Fig. 3. CV values were markedly higher for kurtosis and skewness vs. mean and median values, indicative of higher inter-tumour variation for heterogeneity parameters. SD overall showed similar CV values as the mean and median values, except for a substantial higher CV observed for SD of R₁ pre and post O₂.

**Correlation between MRI methods.** Figure 4 shows heatmaps of significant correlations found between the MRI parameter values in tumours. Numerical values of the correlations are listed in Supplementary Table S1. For the central tendency parameters, significant correlations were found between ADC values and parameters from all other sequences (DCE-MRI, BOLD and TOLD), between DCE-MRI and ΔR₁ and between BOLD and TOLD values. However, no significant correlations were found between BOLD and DCE-MRI central tendency parameter values; while these parameters showed significant correlation with respect to their MRI heterogeneity features, which indicates that these parameters exhibited similar intra-tumour distribution, as illustrated in Fig. 2. Differences in correlations between mean and median MRI parameters at both field strengths were not significant (P-value range mean 0.119–0.949, median 0.098–0.981).

**Correlation with histopathology.** None of the MRI parameters showed a significant difference between histological grades (Kruskal-Wallis P-value range 0.085–0.410). Significant correlations between MRI parameters and histologically stained tumour fractions in 14 lesions (14 patients) are displayed in heatmaps in Fig. 5 (numerical values in Supplementary Table S2). CD3, CD31 and CD68 tumour fractions correlated moderately to strongly
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Inter-tumour variability in MRI parameters. Absolute coefficients of variation (CV) of mean, median, SD, kurtosis and skewness parameter values across 39 HCC tumours (in 32 patients). Kurtosis and skewness showed a substantially higher inter-tumour variability compared to mean and median values. ADC = apparent diffusion coefficient, ART = arterial fraction, DV = distribution volume, F_a = arterial flow, F_p = portal flow, F_t = total flow, MTT = mean transit time, R₁ = longitudinal relaxation rate, R₂* = transverse relaxation rate.</figcaption></figure>

with both central tendency and heterogeneity parameter values. With regards of central tendency parameters, CD3 stained tumour fractions showed moderate significant negative correlations with mean and median R₁ pre O₂ (*P* < 0.029). CD31 staining significantly negatively correlated with mean and median R₂* pre and post O₂ (*P* < 0.009). CD68 stained tumour fractions correlated negatively with mean and median R₁ pre O₂ and median R₁ post O₂ and positively with mean and median ADC (*P* < 0.050). HIF1α tumour fractions did not correlate with any of the central tendency parameters, but did show significant positive correlation with SD of R₂* pre and post O₂ and skewness of R₁ pre O₂ (*P* < 0.048).

**Gene expression analysis.** The distribution of HCC molecular subclasses was as follows: S1 (n = 6, 43%), S2 (n = 3, 21%), S3 (n = 5, 36%). No significant differences between MRI parameters in the different subclasses were observed (P-value range 0.160–0.970). MRI central tendency and heterogeneity features were complementary in terms of correlations with gene expression levels, as illustrated in Fig. 6 (numerical correlation values in Supplementary Table S3). The expression level of CTLA4 only correlated with central tendency parameters, while expression of GLUL, FGFR4, tumour stemness markers KRT19 and EPCAM and immune checkpoint PDCD1 showed significant correlations with MRI heterogeneity parameters only. The expression levels of early HCC
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Heatmaps of correlations between mean, median, SD, kurtosis and skewness of MRI parameters in 39 HCC lesions (in 32 patients). Significant correlations (*P* < 0.05) are coloured according to the scale bar. The correlation between R₁ pre and post O₂ and DCE-MRI parameters was not assessed, because the baseline R₁ measurements were used as input for the DCE-MRI modelling. A combined heatmap of all significant correlations between MRI features is shown on the right, illustrating additional information provided by heterogeneity parameters (SD, kurtosis and skewness). Significant correlations between DCE-MRI and BOLD were for example only seen for heterogeneity parameters and not for central tendency parameters (mean and median). ADC = apparent diffusion coefficient, ART = arterial fraction, DV = distribution volume, F<sub>a</sub> = arterial flow, F<sub>p</sub> = portal flow, F<sub>t</sub> = total flow, MTT = mean transit time, R<sub>1</sub> = longitudinal relaxation rate, R<sub>2</sub>* = transverse relaxation rate.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Heatmaps of correlations between MRI parameters and histopathology in 14 HCC lesions in 14 patients. Significant correlations (*P* < 0.05) are coloured according to the scale bar. A combined heatmap of all significant correlations between MRI features and histopathology is shown on the right, illustrating additional information provided by heterogeneity parameters (SD, kurtosis and skewness). For example, HIF1α tumour fractions correlated with heterogeneity parameters, but not with central tendency parameters (mean and median). ADC = apparent diffusion coefficient, ART = arterial fraction, CD3 = cluster of differentiation 3, CD31 = cluster of differentiation 31, CD68 = cluster of differentiation 68, DV = distribution volume, F<sub>a</sub> = arterial flow, F<sub>p</sub> = portal flow, F<sub>t</sub> = total flow, HIF1α = hypoxia-inducible factor 1-alpha, MTT = mean transit time, R<sub>1</sub> = longitudinal relaxation rate, R<sub>2</sub>* = transverse relaxation rate.</figcaption></figure>

markers *BIRC5*, *HSP70*, *LYVE* and *EZH2*, angiogenesis marker *VEGFA* and immune checkpoint *CD274* significantly correlated with both central tendency and heterogeneity parameters.

## Discussion

In this study we quantified intra- and inter-tumour heterogeneity of HCC lesions using quantitative mpMRI and correlated the mpMRI parameters with histopathological and gene expression analysis in a subset of patients.

Intra-tumour heterogeneity, assessed by the presence of leptokurtosis in the MRI parameters, was observed in all lesions. DCE-MRI and TOLD parameters showed a significantly higher kurtosis compared to ADC, which suggests that intra-tumour heterogeneity in HCC is mainly caused by spatial variation in presence of (functional) vasculature/oxygenation rather than intra-tumoural differences in cellularity.

A higher inter-tumour variation was found for the heterogeneity parameters compared to the mean and median parameter values. This finding suggests that the proposed heterogeneity analysis could potentially allow for a more extensive differentiation between tumours than using the more traditional analysis of central tendency parameter values in ROIs.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Heatmaps of correlations between MRI parameters and gene expression levels in 14 HCC lesions in 14 patients. Significant correlations (*P* < 0.05) are coloured according to the scale bar. A combined heatmap of all MRI features is shown at the figure bottom, illustrating additional information provided by heterogeneity parameters (SD, kurtosis and skewness) compared to central tendency parameters (mean and median) with respect to correlations between MRI and gene expression levels. ADC = apparent diffusion coefficient, ART = arterial fraction, BIRC5 = Baculoviral IAP repeat containing 5, CD274 = cluster of differentiation 274, CTLA4 = cytotoxic T-lymphocyte-associated protein 4, DV = distribution volume, EPCAM = epithelial cell adhesion molecule, EZH2 = enhancer of zeste homolog 2, F<sub>a</sub> = arterial flow, FGFR4 = fibroblast growth factor receptor 4, F<sub>p</sub> = portal flow, F<sub>i</sub> = total flow, GLUL = glutamate-ammonia ligase, HSP70 = 70 kilodalton heat shock protein, KRT19 = keratin 19; LYVE1 = lymphatic vessel endothelial hyaluronan receptor 1, MTT = mean transit time, R<sub>1</sub> = longitudinal relaxation rate, R<sub>2</sub>* = transverse relaxation rate, VEGFA = vascular endothelial growth factor A.</figcaption></figure>

With respect to inter-sequence correlations, we found that DCE-MRI and BOLD MRI central tendency parameters did not correlate, while heterogeneity features of these two sequences showed a significant correlation, indicating that DCE-MRI and BOLD-MRI parameters exhibited a similar intra-tumour distribution. Baudelet *et al*. found that DCE-MRI and BOLD MRI were also largely complementary for the characterization of vasculature in a preclinical fibrosarcoma model<sup>43</sup>. A possible explanation for the lack of a clear correlation between the techniques could be due to the fact that BOLD-MRI has intrinsic low sensitivity in tumour blood vessels with low haematocrit, as opposed to DCE-MRI, in which contrast agent influx is not affected by haematocrit<sup>43</sup>.

Significant correlations were found between BOLD, TOLD and DWI central tendency parameters and histo-pathology. The negative correlations between R<sub>1</sub> pre O<sub>2</sub> and CD3- and CD68-stained tumour fractions and the positive correlation between ADC and CD68 tumour fractions may be explained by a possible increase in interstii-tial fluid associated with inflammation<sup>44</sup>. The negative correlation between R<sub>2</sub>* (pre and post O<sub>2</sub>) and CD31 may be due to the fact that tumours with higher vascularity will most probably be less hypoxic and therefore exhibit a lower R<sub>2</sub>*.

Strong significant correlations were observed between gene expression levels and several MRI parameters. Poor tumour perfusion (low flow/ART and high MTT) was associated with high expression of HCC markers, VEGFA and immune checkpoints CD274 and CTLA4. The relation with HCC markers may be related to tumour progression, which is associated with decreased arterial flow, secondary to increased interstitial pressure and resulting closure of arterial capillaries<sup>45</sup>, and higher expression of HCC markers, including HSP70 and LYVE1<sup>46, 47</sup>. The negative correlation between perfusion and VEGFA expression could potentially be explained by the fact that decreased perfusion may be associated with hypoxia, which triggers expression of VEGFA<sup>48</sup>. However, counterin-tuitively, VEGFA expression was not associated with mean or median BOLD parameters and negatively associated with ΔR<sub>1</sub>, which suggests absence of correlation between MRI-measured hypoxia and expression of VEGFA. The underlying physiology of the negative correlations between tumour perfusion and CTLA4 and CD274 expression is unclear and needs to be verified in a larger study. The expression of immune checkpoints has shown to be positively associated with tumour progression in lung cancer<sup>49</sup>, but no data has been reported on this association in HCC. Interestingly, significant correlations with expression of several genes, including stemness markers and immunotherapy target PDCD1, were only observed with heterogeneity parameters, which further illustrates the complementary properties of central tendency and heterogeneity parameters in terms of tumour characterization.
In general, limited data is available on radiogenomics in HCC. There are a few reports that compare CT imaging traits with HCC genomics analysis¹⁰,¹¹. However, these studies have assessed only qualitative metrics, which are limited by inter- and intra-observer variability. The quantitative assessment of heterogeneity histogram features of functional mpMRI in the current study is probably less prone to inter-observer variability. Nevertheless, the proposed heterogeneity features can potentially be confounded by image noise, which introduces non-biological heterogeneity, and by inaccuracies in parameter estimation of single pixels³⁰. However, we found that the HCC lesions exhibited more intra-tissue heterogeneity than liver parenchyma, which indicates that the observed intra-tumour heterogeneity in the HCC lesions is likely biological in nature.

There are a few published reports on the use of imaging to predict and assess immunotherapy outcome⁵¹. Jajamovich *et al*. found a significant correlation between ADC and a gene signature related to dendritic cell maturation in glioblastoma, which can be used to stratify patients with immunogenic tumours for immunotherapy⁵². Mayerhoefer *et al*. found that DWI is a promising technique for the evaluation of immunotherapy response in lymphoma patients⁵³. The potential value of DWI for immunotherapy was also reported by Qin *et al.*, who found that the decrease or stabilization of the lesion volume on ADC maps, after an initial increase after immunotherapy, was predictive of therapeutic benefit in glioblastoma⁵⁴. Interestingly, we did not find significant correlations between ADC and any of the immunotherapy targets. The role of DWI for the prediction of immunotherapy outcome in HCC needs to be established in future studies. There is no published data assessing HCC response to immunotherapy, given the recent introduction of this therapy in HCC. More studies are needed to determine the exact value of mpMRI for the prediction of immunotherapy efficacy in HCC and other cancer types.

While mpMRI proved promising for the non-invasive assessment of histopathological and genomics properties of HCC lesions in our study, the separate imaging, histopathological and genomics measurements also provided independent information on tumour properties. For example, none of the mpMRI features could distinguish between different histopathological grades. We therefore do not believe that mpMRI could fully replace histopathologic assessment in HCC. The mpMRI analysis could rather be used to improve the characterization of HCC lesions in biopsied or resected samples, by providing mutual and additional information on tumour properties on a whole-tumour level¹². Nevertheless, in the majority of patients HCC is diagnosed based on imaging alone, without histologic confirmation⁵⁵. Treatment stratification in those patients without biopsy could ultimately benefit from knowledge of the correlation of imaging parameters with histopathological and genomics properties of HCC lesions in a training set of resected HCC samples.

The findings described in our study may have substantial clinical implications for characterization and treatment stratification of HCC. Tumours with high intra-tumour heterogeneity generally have a poorer prognosis, which may be secondary to a more aggressive biology or treatment resistance⁵⁶. The observed correlations between MRI features and the expression of potential therapeutic markers for targeted molecular or immunotherapy of HCC⁵⁷ suggest that the proposed mpMRI histogram analysis could be used to noninvasively predict treatment outcome at baseline. Knowledge of the expression of therapeutic targets using non-invasive imaging could aid in the stratification of patient-tailored, personalized treatment. A future study should assess the proposed mpMRI heterogeneity analysis as a predictor of response to immunotherapy. In such study, it could also be assessed which (combination of) MRI features have the highest diagnostic performance for prediction of treatment response.

Our study has several limitations. Assessment of intra-tumour heterogeneity was not feasible for the histopathology and genomics analysis, because only a small part of the lesion was available for these analyses. Another limitation was the relatively small sample size of HCC lesions for histopathology and genomics. Moreover, there was variability between MRI field strengths on which the patients were scanned, which could affect MRI parameter estimation. Nevertheless, the correlations between MRI parameters were not significantly different between both field strengths. Finally, the correlation of the mpMRI features with clinical parameters such as survival and treatment response was not assessed in this initial study.

In conclusion, our results show that central tendency and heterogeneity features of mpMRI data of HCC tumours are complementary in terms of correlations between sequences, with histopathology and with gene expression levels. The proposed histogram analysis is therefore promising for non-invasive HCC characterization on the imaging, immunohistochemical and genomics level.

## References

1. Abramson, R. G. et al. MR Imaging Biomarkers in Oncology Clinical Trials. Magn. Reson. Imaging Clin. N. Am. 24, 11–29, doi:10.1016/j.mric.2015.08.002 (2016).

2. Guo, Y. et al. Differentiation of clinically benign and malignant breast lesions using diffusion-weighted imaging. J. Magn. Reson. Imaging 16, 172–178, doi:10.1002/jmri.10140 (2002).

3. Bane, O. et al. Feasibility and reproducibility of BOLD and TOLD measurements in the liver with oxygen and carbogen gas challenge in healthy volunteers and patients with hepatocellular carcinoma. J. Magn. Reson. Imaging 43, 866–876, doi:10.1002/jmri.25051 (2016).

4. Yankeelev, T. F. & Gore, J. C. Dynamic Contrast Enhanced Magnetic Resonance Imaging in Oncology: Theory, Data Acquisition, Analysis, and Examples. Curr. Med. Imaging Rev. 3, 91–107, doi:10.2174/157340507780619179 (2009).

5. Just, N. Improving tumour heterogeneity MRI assessment with histograms. Br. J. Cancer 111, 2205–2213, doi:10.1038/bjc.2014.512 (2014).

6. Anfuso, B., El-Khobar, K. E., Sukowati, C. H. & Tiribelli, C. The multiple origin of cancer stem cells in hepatocellular carcinoma. Clin. Res. Hepatol. Gastroenterol. 39 (Suppl 1), S92–97, doi:10.1016/j.clinre.2015.05.011 (2015).

7. Friemel, J. et al. Intratumor heterogeneity in hepatocellular carcinoma. Clin. Cancer Res. 21, 1951–1961, doi:10.1158/1078-0432.CCR-14-0122 (2015).

8. Fujita, N. et al. Hepatocellular carcinoma: clinical significance of signal heterogeneity in the hepatobiliary phase of gadoxetic acid-enhanced MR imaging. Eur. Radiol. 25, 211–220, doi:10.1007/s00330-014-3349-9 (2015).

9. Rutman, A. M. & Kuo, M. D. Radiogenomics: creating a link between molecular diagnostics and diagnostic imaging. Eur. J. Radiol. 70, 232–241, doi:10.1016/j.ejrad.2009.01.050 (2009).
10. Kuo, M. D., Gollub, J., Sirlin, C. B., Ooi, C. & Chen, X. Radiogenomic analysis to identify imaging phenotypes associated with drug response gene expression programs in hepatocellular carcinoma. *J. Vasc. Interv. Radiol.* **18**, 821–831, doi:10.1016/j.jvir.2007.04.031 (2007).

11. Segal, E. et al. Decoding global gene expression programs in liver cancer by noninvasive imaging. *Nat. Biotechnol.* **25**, 675–680, doi:10.1038/nbt1306 (2007).

12. Gillies, R. J., Kinahan, P. E. & Hricak, H. Radiomics: Images Are More than Pictures, They Are Data. *Radiology* **278**, 563–577, doi:10.1148/radiol.2015151169 (2016).

13. Pardoll, D. M. The blockade of immune checkpoints in cancer immunotherapy. *Nat. Rev. Cancer* **12**, 252–264, doi:10.1038/nrc3239 (2012).

14. Sharma, P. & Allison, J. P. The future of immune checkpoint therapy. *Science* **348**, 56–61, doi:10.1126/science.aaa8172 (2015).

15. Calabro, L. et al. Tremelimumab for patients with chemotherapy-resistant advanced malignant mesothelioma: an open-label, single-arm, phase 2 trial. *Lancet Oncol* **14**, 1104–1111, doi:10.1016/S1470-2045(13)70381-4 (2013).

16. Sangro, B. et al. A clinical trial of CTLA-4 blockade with tremelimumab in patients with hepatocellular carcinoma and chronic hepatitis C. *J Hepatol* **59**, 81–88, doi:10.1016/j.jhep.2013.02.022 (2013).

17. OPTN policies, policy 9: Allocation of Livers and Liver-Intestines; 9.3.F Candidates with Hepatocellular Carcinoma. (2017).

18. Hectors, S. J. et al. Intravoxel incoherent motion diffusion-weighted imaging of hepatocellular carcinoma: Is there a correlation with flow and perfusion metrics obtained with dynamic contrast-enhanced MRI? *J. Magn. Reson. Imaging* **44**, 856–864, doi:10.1002/jmri.25194 (2016).

19. Jajamovich, G. H. et al. DCE-MRI of hepatocellular carcinoma: perfusion quantification with Tofts model versus shutter-speed model–initial experience. *MAGMA* **29**, 49–58, doi:10.1007/s10334-015-0513-4 (2016).

20. Jajamovich, G. H., Dyvorne, H., Donnerhack, C. & Taouli, B. Quantitative liver MRI combining phase contrast imaging, elastography, and DWI: assessment of reproducibility and postprandial effect at 3.0 T. *PLoS One* **9**, e97355, doi:10.1371/journal.pone.0097355 (2014).

21. Dyvorne, H., Jajamovich, G., Kakite, S., Kuehn, B. & Taouli, B. Intravoxel incoherent motion diffusion imaging of the liver: optimal b-value subsampling and impact on parameter precision and reproducibility. *Eur. J. Radiol.* **83**, 2109–2113, doi:10.1016/j.ejrad.2014.09.003 (2014).

22. Dyvorne, H. A. et al. Diffusion-weighted imaging of the liver with multiple b values: effect of diffusion gradient polarity and breathing acquisition on image quality and intravoxel incoherent motion parameters–a pilot study. *Radiology* **266**, 920–929, doi:10.1148/radiol.12120686 (2013).

23. Deichmann, R. & Haase, A. Quantification of T1 Values by Snapshot-Flash Nmr Imaging. *J. Magn. Reson.* **96**, 608–612, doi:10.1016/0022-2364(92)90347-A (1992).

24. Haase, A. Snapshot FLASH MRI. Applications to T1, T2, and chemical-shift imaging. *Magn. Reson. Med.* **13**, 77–89, doi:10.1002/(ISSN)1522-2594 (1990).

25. Aronhime, S. et al. DCE-MRI of the liver: effect of linear and nonlinear conversions on hepatic perfusion quantification and reproducibility. *J. Magn. Reson. Imaging* **40**, 90–98, doi:10.1002/jmri.24341 (2014).

26. Zhang, X. et al. *In vivo* blood T(1) measurements at 1.5 T, 3 T, and 7 T. *Magn. Reson. Med.* **70**, 1082–1086, doi:10.1002/mrm.24550 (2013).

27. Pintaske, J. et al. Relaxivity of Gadopentetate Dimeglumine (Magnevist), Gadobutrol (Gadovist), and Gadobenate Dimeglumine (MultiHance) in human blood plasma at 0.2, 1.5, and 3 Tesla. *Invest. Radiol.* **41**, 213–221, doi:10.1097/01.rli.0000197668.44926.f7 (2006).

28. Materne, R. et al. Assessment of hepatic perfusion parameters with dynamic MRI. *Magn. Reson. Med.* **47**, 135–142, doi:10.1002/(ISSN)1522-2594 (2002).

29. Yoon, S. H. et al. Tumor Heterogeneity in Lung Cancer: Assessment with Dynamic Contrast-enhanced MR Imaging. *Radiology*, 151367, doi:10.1148/radiol.2016151367 (2016).

30. Edmondson, H. A. & Steiner, P. E. Primary carcinoma of the liver: a study of 100 cases among 48,900 necropsies. *Cancer* **7**, 462–503, doi:10.1002/(ISSN)1097-0142 (1954).

31. Remark, R. et al. In-depth tissue profiling using multiplexed immunohistochemical consecutive staining on single slide. *Science Immunology*, doi:10.1126/sciimmunol.aaf6925 (2016).

32. Ziello, J. E., Jovin, I. S. & Huang, Y. Hypoxia-Inducible Factor (HIF)-1 regulatory pathway and its potential for therapeutic intervention in malignancy and ischemia. *Yale J. Biol. Med.* **80**, 51–60 (2007).

33. Tan, P. S. et al. Clinicopathological indices to predict hepatocellular carcinoma molecular classification. *Liver Int* **36**, 108–118, doi:10.1111/liv.12889 (2016).

34. Hoshida, Y. et al. Integrative transcriptome analysis reveals common molecular subclasses of human hepatocellular carcinoma. *Cancer Res.* **69**, 7385–7392, doi:10.1158/0008-5472.CAN-09-1089 (2009).

35. Goossens, N., Sun, X. & Hoshida, Y. Molecular classification of hepatocellular carcinoma: potential therapeutic implications. *Hepat Oncol* **2**, 371–379, doi:10.2217/hep.15.26 (2015).

36. Villanueva, A. et al. Combining clinical, pathology, and gene expression data to predict recurrence of hepatocellular carcinoma. *Gastroenterology* **140**, 1501–1512 e1502, doi:10.1053/j.gastro.2011.02.006 (2011).

37. Yamashita, T. et al. EpCAM and alpha-fetoprotein expression defines novel prognostic subtypes of hepatocellular carcinoma. *Cancer Res.* **68**, 1451–1461, doi:10.1158/0008-5472.CAN-07-6013 (2008).

38. Llovet, J. M. et al. A molecular signature to discriminate dysplastic nodules from early hepatocellular carcinoma in HCV cirrhosis. *Gastroenterology* **131**, 1758–1767, doi:10.1053/j.gastro.2006.09.014 (2006).

39. Di Tommaso, L. et al. Diagnostic value of HSP70, glypican 3, and glutamine synthetase in hepatocellular nodules in cirrhosis. *Hepatology* **45**, 725–734, doi:10.1002/hep.21531 (2007).

40. Hagel, M. et al. First Selective Small Molecule Inhibitor of FGFR4 for the Treatment of Hepatocellular Carcinomas with an Activated FGFR4 Signaling Pathway. *Cancer Discov.* **5**, 424–437, doi:10.1158/2159-8290.CD-14-1029 (2015).

41. Horwitz, E. et al. Human and mouse VEGFA-amplified hepatocellular carcinomas are highly sensitive to sorafenib treatment. *Cancer Discov* **4**, 730–743, doi:10.1158/2159-8290.CD-13-0782 (2014).

42. Harding, J. J., El Dika, I. & Abou-Alfa, G. K. Immunotherapy in hepatocellular carcinoma: Primed to make a difference? *Cancer* **122**, 367–377, doi:10.1002/cncr.29769 (2016).

43. Baudelet, C., Cron, G. O. & Gallez, B. Determination of the maturity and functionality of tumor vasculature by MRI: correlation between BOLD-MRI and DCE-MRI using P792 in experimental fibrosarcoma tumors. *Magn. Reson. Med.* **56**, 1041–1049, doi:10.1002/mrm.21047 (2006).

44. Baronzio, G., Parmar, G., Baronzio, M. & Kiselevsky, M. Tumor interstitial fluid: proteomic determination as a possible source of biomarkers. *Cancer Genomics Proteomics* **11**, 225–237 (2014).

45. Choi, J. Y., Lee, J. M. & Sirlin, C. B. CT and MR imaging diagnosis and staging of hepatocellular carcinoma: part I. Development, growth, and spread: key pathologic and imaging aspects. *Radiology* **272**, 635–654, doi:10.1148/radiol.14132361 (2014).

46. Joo, M., Chi, J. G. & Lee, H. Expressions of HSP70 and HSP27 in hepatocellular carcinoma. *J. Korean Med. Sci.* **20**, 829–834, doi:10.3346/jkms.2005.20.5.829 (2005).

47. Kitagawa, K., Nakajima, G., Kuramochi, H., Ariizumi, S. I. & Yamamoto, M. Lymphatic vessel endothelial hyaluronan receptor-1 is a novel prognostic indicator for human hepatocellular carcinoma. *Mol Clin Oncol* **1**, 1039–1048, doi:10.3892/mco.2013.167 (2013).
48. Liu, Y., Cox, S. R., Morita, T. & Kourembanas, S. Hypoxia regulates vascular endothelial growth factor gene expression in endothelial cells. Identification of a 5' enhancer. *Circ. Res.* 77, 638–643 (1995).

49. Thommen, D. S. et al. Progression of Lung Cancer Is Associated with Increased Dysfunction of T Cells Defined by Coexpression of Multiple Inhibitory Receptors. *Cancer Immunol Res* 3, 1344–1355, doi:10.1158/2326-6066.CIR-15-0097 (2015).

50. Asselin, M. C., O'Connor, J. P., Boellaard, R., Thacker, N. A. & Jackson, A. Quantifying heterogeneity in human tumours using MRI and PET. *Eur. J. Cancer* 48, 447–455, doi:10.1016/j.ejca.2011.12.025 (2012).

51. Garcia-Figueiras, R., Padhani, A. R. & Baleato-Gonzalez, S. Therapy Monitoring with Functional and Molecular MR Imaging. *Magn. Reson. Imaging Clin. N. Am.* 24, 261–288, doi:10.1016/j.mric.2015.08.003 (2016).

52. Jajamovich, G. H., Valiathan, C. R., Cristescu, R. & Somayajula, S. Integrative analysis of diffusion-weighted MRI and genomic data to inform treatment of glioblastoma. *J. Neurooncol.* 129, 289–300, doi:10.1007/s11060-016-2174-1 (2016).

53. Mayerhoefer, M. E. et al. Evaluation of Diffusion-Weighted Magnetic Resonance Imaging for Follow-up and Treatment Response Assessment of Lymphoma: Results of an 18F-FDG-PET/CT-Controlled Prospective Study in 64 Patients. *Clin. Cancer Res.* 21, 2506–2513, doi:10.1158/1078-0432.CCR-14-2454 (2015).

54. Qin, L. et al. Advanced MRI assessment to predict benefit of anti-programmed cell death 1 protein immunotherapy response in patients with recurrent glioblastoma. *Neuroradiology* 59, 135–145, doi:10.1007/s00234-016-1769-8 (2017).

55. McEvoy, S. H. et al. Hepatocellular carcinoma: illustrated guide to systematic radiologic diagnosis and staging according to guidelines of the American Association for the Study of Liver Diseases. *Radiographics* 33, 1653–1668, doi:10.1148/rg.336125104 (2013).

56. Davnall, F. et al. Assessment of tumor heterogeneity: an emerging imaging tool for clinical practice? *Insights Imaging* 3, 573–589, doi:10.1007/s13244-012-0196-6 (2012).

57. Connell, L. C., Harding, J. J. & Abou-Alfa, G. K. Advanced Hepatocellular Cancer: the Current State of Future Research. *Curr. Treat. Options Oncol.* 17, 43, doi:10.1007/s11864-016-0415-3 (2016).

## Acknowledgements

MW: Fondation ARC Contract grant number: SAE20140601302. SL: RSNA Research Seed Grant #RSD1608. SG: Clinical Strategy Teams Grants from the Cancer Research Institute; NIH Grant P01 CA190174. MM: NIH Grants R01 CA173861, U19 AI 128949, U24 AI 118644, U19 AI 117873. YH: NIH Grant R01 DK099558, EU Grant ERC-2014-AdG-671231 HEPCIR, Irma T. Hirschl Trust, and U.S. Department of Defense W81XWH-16-1-0363. BT: NIH Grants U01 CA172320, 1R01DK087877.

## Author Contributions

Study design: S.H., M.W., B.T. Development of methodology: S.H., M.W., O.B., R.R., S.G., M.M., Y.H., B.T. Acquisition of data: S.H., M.W., O.B., C.B., R.R., N.C., M.I.F., H.Z., Y.H. Analysis of data: S.H., M.W., C.B., Y.H. Writing and review of manuscript: S.H., M.W., O.B., C.B., S.L., R.R., S.G., Y.H., B.T.

## Additional Information

Supplementary information accompanies this paper at doi:10.1038/s41598-017-02706-z

**Competing Interests:** The authors declare that they have no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2017