Diffusion MRI signal cumulants and hepatocyte microstructure at fixed diffusion time: Insights from simulations, 9.4T imaging, and histology.

PURPOSE: Relationships between diffusion-weighted MRI signals and hepatocyte microstructure were investigated to inform liver diffusion MRI modeling, focusing on the following question: Can cell size and diffusivity be estimated at fixed diffusion time, realistic SNR, and negligible contribution from extracellular/extravascular water and exchange?
METHODS: Monte Carlo simulations were performed within synthetic hepatocytes for varying cell size/diffusivity L / D 0 , and clinical protocols (single diffusion encoding; maximum b-value: {1000, 1500, 2000} s/mm2 ; 5 unique gradient duration/separation pairs; SNR = { ∞ , 100, 80, 40, 20}), accounting for heterogeneity in ( D 0 , L ) and perfusion contamination. Diffusion ( D ) and kurtosis ( K ) coefficients were calculated, and relationships between ( D 0 , L ) and ( D , K ) were visualized. Functions mapping ( D , K ) to ( D 0 , L ) were computed to predict unseen ( D 0 , L ) values, tested for their ability to classify discrete cell-size contrasts, and deployed on 9.4T ex vivo MRI-histology data of fixed mouse livers
RESULTS: Relationships between ( D , K ) and ( D 0 , L ) are complex and depend on the diffusion encoding. Functions mapping D , K to ( D 0 , L ) captures salient characteristics of D 0 ( D , K ) and L ( D , K ) dependencies. Mappings are not always accurate, but they enable just under 70% accuracy in a three-class cell-size classification task (for SNR = 20, b max = 1500 s/mm2 , δ = 20 ms, and Δ = 75 ms). MRI detects cell-size contrasts in the mouse livers that are confirmed by histology, but overestimates the largest cell sizes.
CONCLUSION: Salient information about liver cell size and diffusivity may be retrieved from minimal diffusion encodings at fixed diffusion time, in experimental conditions and pathological scenarios for which extracellular, extravascular water and exchange are negligible.

INTRODUCTION

Diffusion‐weighted (DW) MRI relies on the self‐diffusion of water residing in biological tissues to probe cellular microarchitecture. In classical pulsed gradient spin echo, , two diffusion gradients sensitize the acquisition to diffusion. The first gradient tags spin phases according to spatial position, whereas the second one, played out after a certain interval (known as diffusion time), cancels such tags for stationary spins. In the presence of diffusion, water molecules change their position during the diffusion time, and the tag removal is incomplete. This leads to MRI signal attenuation, which carries a signature of tissue microstructure.

Model‐based methods offer practical solutions to the estimation of microenvironment properties from MRI by adopting geometric models of microstructure. This leads to tractable expressions that parametrize the signal as a function of sequence and microstructural parameters. So far, model‐based methods have found several clinical applications, , , , , , , , in spite of potential biases occurring as a result of modeling oversimplifications. , Modeling has focused on neural , , , , , , and prostate , , tissue characterization, as well on cell‐size measurements, , , , relevant in oncology. However, less attention has been paid to other organs, such as the liver. , , , Biologically specific DW MRI methods are urgently required in liver diseases, such as liver cancer, a leading cause of cancer‐related death. Liver cancer (either primary or metastatic ) shows a variety of microstructural characteristics. Quantitative liver MRI methods offer sensitivity to cancer pathology, but still fail to distinguish key pathological differences (e.g., substitution of either sinusoidal endothelial cells or liver hepatocytes by neoplastic cells , ). There is a pressing need for new clinically viable liver MRI readouts; these could help reduce the use of invasive biopsies, which sparsely sample the tissue, are prone to false negatives, and can result in complications for the patient, and could support diagnosis and treatment selection.

A key step in diffusion MRI development is the identification of microstructural features that can be estimated from clinical‐like (i.e., intermediate b‐values and limited scan time) measurements. To our knowledge, such a characterization for hepatocytes, which account for up to 85% of liver volume, is still lacking. Here we considered realistic hepatocyte sizes and diffusion protocols that could be feasible in the clinic (single diffusion encoding, maximum b‐value up to 2000 s/mm2, fixed diffusion time with gradient separation/duration / in the range of [25; 75] ms and [10; 40] ms, SNR as low as 20 at b = 0). Through Monte Carlo simulations and co‐localized 9.4T ex vivo MRI and histology of fixed mouse livers, we specifically investigated the following question: Can cell size and diffusivity be estimated from signal cumulants at fixed diffusion time and realistic SNR, under the assumption of negligible contributions from extracellular/extravascular water and water exchange? While experiments performed at varying diffusion times are ideal for cell‐size measurement, , , , , techniques providing summary cell‐size indices with minimal acquisitions have the potential of bringing quantitative MRI one step closer to the clinic.

METHODS

We simulated intracellular signals at fixed diffusion time and processed them to estimate cell size and cell diffusivity . The approach was also tested on 9.4T ex vivo MRI scans of fixed mouse livers. All analysis code is made available ( ), and was executed on two Ubuntu 20.04.2 machines (18‐core, 3.00‐GHz Intel® Core i9‐10980XE CPU).

Cell generation

We simulated hepatocytes (polygonal cells ) by perturbing regular prisms with square/pentagonal/hexagonal bases. Prisms were described by triangular meshes and featured a characteristic length L (base‐to‐base height and diameter of the circumcircle relative to each base). We considered 33 values of L in [11; 60] μm (increment: 1.5 μm), obtaining 15 unique cell shapes for each value of L. The cells at fixed L were obtained by perturbing each prism base shape 5 times, displacing vertices at random (displacements drawn from a normal distribution, σ = 0.1L). The range for L covers sizes seen in healthy mammal livers (e.g., 20–30 μm in humans, 30–40 μm in mice ) and in pathology (e.g., swollen hepatocytes in steatosis ; hepatocyte substitution by smaller cancer cells ). Supporting Information Figure S1 shows synthetic cells.

Intracellular spin dynamics

We generated random walks with the MCDC simulator, distributing = 1000 spins uniformly inside each cell (elastic reflection at walls; impermeable walls). We simulated = 140 ms (3000 steps) and varied the intrinsic cell diffusivity in [0.20; 2.40] (45 values; increment: 0.05 ).

Magnetic resonance imaging signal synthesis

For each fixed (D 0, L) value, we pooled together spin trajectories simulated within a neighborhood of (i.e., . This introduces heterogeneity expected in realistic voxels, , leading to 1189 pairs. Above, = 1, …, is the index of a spin within a cell; = 1, …, is the cell‐shape index for fixed cell size; and enumerates the elements of , with . For MRI signal synthesis, we considered single diffusion encoding gradient waveforms , with five unique clinically realistic gradient duration/separation δ/Δ ([10 ms, 50 ms], [20 ms, 25 ms], [20 ms, 50 ms], [20 ms, 75 ms], and [40 ms, 50 ms]). For any fixed (δ, Δ), we synthesized measurements corresponding to seven nonzero b‐values, uniformly spaced in (;), where = 100 s mm–2, a value used to suppress intravoxel incoherent motion (IVIM)–like components , ; and = (1000, 1500, 2000) s mm–2, as the volume‐weighted sum:

For each b‐value, we generated signals for three mutually orthogonal gradients (as common in liver MRI , ), averaged them, and introduced random slow‐flow (intravoxel incoherent motion, or IVIM) , contamination as follows:

In this equation, controls the IVIM contamination, and . , Rician noise was injected at an SNR of (, 100, 80, 40, 20), where denotes no noise added; , denotes the noise variance.

Analysis

We estimated apparent diffusion/kurtosis coefficients / for any , diffusion protocol, and SNR by fitting , through constrained nonlinear least‐squares fitting initialized by linear fitting (; ; ). The value of is the non‐DW signal.

We tested whether can be estimated from and when the contribution of extracellular, extravascular water and transcytolemmal water exchange are negligible. To this end, we related to using color‐coded scatter plots, and studied paired and observations estimating smooth function (, ) mapping to at a fixed protocol and SNR. The estimation was based on the following polynomial functions: where m indicates and in turn. We refer to the estimation of and via Equation 4 as PolyMap. Coefficients a were estimated on 700 randomly selected training pairs out of 1189, and then deployed to predict the 489 unseen values. We assessed the quality of the prediction by visualizing errors against ground‐truth values. For reference, PolyMap was compared with fitting of a biophysical model of the intracellular DW signal (SigFit estimation) as follows:

Equation 5 relies on an approximate expression of the apparent diffusion coefficient for spins diffusing in a bounded medium (wide‐pulse limit). , , The values of and are constants that depend on the geometry: Analysis of intracellular diffusion coefficients from 400 unique coefficients provides , for our synthetic cells. Note that and in Equation 5 are fitted jointly to sets of signal measurements performed at varying b‐value (but fixed and ). This implies that we do not get a single number for the apparent diffusion coefficient value first, and derive and from it afterwards.

Finally, we tested whether it is possible to resolve cell‐size contrasts with the minimal protocols considered here, being that and are difficult to disentangle. We discretized as small (), medium (), large (), and fitted a multinomial logistic regression model (Python statsmodels) in the same form of Equation 4 (: discretised ). The model was fitted to the training set for all diffusion protocols and SNR = 20, and deployed on the validation set. We calculated classification accuracy and estimated 95% accuracy ranges compatible with chance by training on 1000 random permutations of the labels.

Magnetic resonance imaging ‐ histology comparison

Two formalin‐fixed NOD.Cg‐Prkdcscid IL2rgtm1WjI/SzJ mouse livers from an approved, ongoing study (wild‐type [WT] and patient‐derived xenograft (PDX), subcutaneous implantation of prostate cancer bone biopsy) were scanned in phosphate‐buffered saline on a 9.4T Bruker Avance system (room temperature) to test whether our approach can detect histologically meaningful cell‐size differences due to pathology. The DW spin‐echo scans ( = 30 ms; = 10 ms; TE = 45 ms; TR = 2700 ms; 10 b‐values in [0; 4500] s/mm2; two slices, 1‐mm thick; 349 × 273 μm2 resolution) were acquired and preprocessed. , , Images acquired at b > 1700 s/mm2 (i.e., with negligible phosphate‐buffered saline contamination) were analyzed with PolyMap and SigFit. For PolyMap computation, the mapping was learned on signals synthesized for the specific protocol used ex vivo, and corrupted at an SNR equal to the sample median SNR at b = 0, estimated through Marchenko and Pastur principal component analysis. , SigFit fitting was instead performed by either (1) estimating jointly and , or (2) fixing to {0.5, 0.75, 1.0, 1.25, 1.50} in turn to all voxels and then estimating , as in some model‐based approaches.

One 4‐μm‐thick histological section was obtained for each MRI slice, stained with hematoxylin and eosin, and digitized (Hamamatsu C9600‐12 scanner; resolution: 0.227 μm). Cells were segmented with QuPath, obtaining cell‐wise diameters . These were analyzed within patches matching the in‐plane MRI resolution, deriving per‐patch histological cell size

Equation 6 is justified by noting that the total intracellular MRI signal is approximately proportional to , being that is the volume‐weighted sum of individual cell signals and that , implying that was warped to MRI (symmetric diffeomorphic registration of specimens’ manual outlines), and metric distributions were evaluated.

RESULTS

The computation time required to process one MRI protocol was approximately 700 seconds for each pair on one CPU.

Plots in Figure 1 scatter against . Points in the plane correspond to a unique combination, and are colored according to (top) and (bottom). The figure refers to maximum b‐value of 2000 s/mm2, SNR = 20, and multiple combinations of (δ,Δ). The values of and exhibit a wide range of variation (e.g., negative values are seen), depending on protocol δ and Δ. Nonetheless, a trend in the / coloring can be seen (more apparent as SNR increases). In absence of noise, a non‐monotonic relationship between and both and is seen, with points distributed according to complex patterns in the domain (Supporting Information Figure S2). For some specific combinations of , no experimental points are observed. The position of the points in the plane changes depending on (δ,Δ). For example, when = 2000 s/mm2, the median/95% range of are 0.76/[0.27; 1.57] for δ/Δ = 20/25 ms and 0.52 [0.16; 1.20] for δ/Δ = 20/75 ms. For the same gradient timings, median/95% ranges of are 0.44 [0.29; 1.53] and 0.40 [0.07; 2.99]. Results for smaller maximum b‐values follow similar trends, although numerical values of depend on , both in absence or presence of noise (e.g., SNR = 20) (Supporting Information Figure S3 for = 1000 s/mm2). In absence of noise and when = 1000 s/mm2, the median/95% range of values are 0.79 [0.28; 1.61] for δ/Δ = 20/25 ms and 0.55 [0.17; 1.25] for δ/Δ = 20/75 ms, whereas it is 0.58 [0.29; 3.12] and 0.78 [0.14; 6.33] for , larger than what is seen when = 2000 s/mm2.

FIGURE 1

Scatter plots of color‐coded by the underlying average intrinsic cell diffusivity (top, A‐E) and cell size (bottom, F‐J), as obtained when noise is added to the synthetic MRI signals at an SNR at b = 0 of 20. From left to right: Different diffusion times (δ/Δ = 20/25 ms in [A] and [F]; δ/Δ = 40/50 ms in [B] and [G]; δ/Δ = 20/50 ms in [C] and [H]; δ/Δ = 10/50 ms in [D] and [I]; and δ/Δ = 20/75 ms in [E] and [J]). The figure refers to a minimum/maximum protocol b‐value of = 100/2000 s/mm2. Noise‐free intracellular diffusion‐weighted (DW) signals are contaminated by intravoxel incoherent motion (IVIM)–like partial volume

Figure 2 shows scatter plots color‐coded by and for observations belonging to the validation set ( = 2000 s/mm2, δ = 20 ms, Δ = 75 ms, SNR = 20). It also shows and predicted in correspondence of the same values with both PolyMap and SigFit. Similar plots for the noise‐free case and different protocols (e.g., = 1000 s/mm2, δ = 20 ms, Δ = 25 ms) are reported in Supporting Information Figures S4 and S5. Overall, PolyMap does not necessarily predict accurately the values of and on unseen data, especially at high noise level. Nonetheless, it captures the salient characteristics of the and relationships, which appear unique to each diffusion‐encoding protocols and SNR. The SigFit estimation also captures and contrasts, although predictions are less smooth than those from PolyMap.

FIGURE 2

Examples of predictions of intrinsic cell diffusivity and cell size on the validation set. (A–C) Scatter plots colored by cell diffusivity . (D–F) Scatter plots colored by cell size . Left: Signal cumulants at fixed diffusion time colored by underlying ground truth and . Middle: Signal cumulants at fixed diffusion time colored by predictions of and as obtained with PolyMap. Right: Signal cumulants at fixed diffusion time colored by predictions of and as obtained with SigFit. The figure refers to the case when the minimum/maximum protocol b‐values is equal to = 100/2000 s/mm2 and the diffusion gradient duration/separation is δ = 20 ms/Δ = 75 ms, for SNR of 20 and in presence of IVIM contamination

Figure 3 plots PolyMap and prediction errors (prediction – ground truth) against ground‐truth and for different protocols ( = 2000 s/mm2, SNR = 20). The same plots corresponding to SigFit are reported in Supporting Information Figure S6. Further PolyMap and SigFit prediction errors for the noise‐free case and for = 1000 s/mm2, SNR = 20, are included in Supporting Information Figures S7 and S8. The charts reveal that and are overestimated/underestimated at the lower/upper end of their ranges. This trend is observed for different gradient timings and in absence of noise, although to a lesser extent. Higher SNR and longer diffusion times lead to smaller errors. The PolyMap errors are slightly smaller/larger than those from SigFit.

FIGURE 3

and prediction errors for the PolyMap estimation method, scattered against ground‐truth values of and , for different diffusion gradient timings at a fixed maximum b‐value of 2000 s/mm2 and SNR of 20. From left to right: Different gradient timings (δ/Δ = 20/25 ms in [A] and [F]; δ/Δ = 40/50 ms in [B] and [G]; δ/Δ = 20/50 ms in [C] and [H]; δ/Δ = 10/50 ms in [D] and [I]; and δ/Δ = 20/75 ms in [E] and [J]). Plots on top (A‐E) refer to ; plots on the bottom (F‐J) refer to . For each fixed value of (on top, or to the bottom), median errors with interquartile ranges for varying (on top, or varying to the bottom) are also reported

Table 1 reports validation‐set accuracies for the cell‐size classification task. Accuracy values can be as high as almost 70%, such as when = 1500 s/mm2, δ = 20 ms, and Δ = 75 ms, corresponding to 86%, 46%, and 61% correctly classified small, medium, and large cells. Accuracies are above accuracy ranges compatible with chance.

TABLE 1

Accuracies obtained on the validation set for the three‐class cell‐size classification task performed using multinomial logistic regression at SNR = 20 and for all diffusion‐encoding protocols

		b max = 1000 s/mm2	b max = 1500 s/mm2	b max = 2000 s/mm2
δ = 20 ms, Δ = 25 ms, Δ – δ/3 = 18.3 ms	Accuracy	0.54	0.54	0.55
	95% random interval	[0.23; 0.42]	[0.22; 0.43]	[0.22; 0.43]
δ = 40 ms, Δ = 50 ms, Δ – δ/3 = 36.7 ms	Accuracy	0.60	0.61	0.61
	95% random interval	[0.21; 0.46]	[0.21; 0.46]	[0.19; 0.48]
δ = 20 ms, Δ = 50 ms, Δ – δ/3 = 43.3 ms	Accuracy	0.58	0.56	0.60
	95% random interval	[0.21; 0.45]	[0.21; 0.45]	[0.20; 0.46]
δ = 10 ms, Δ = 50 ms, Δ – δ/3 = 46.7 ms	Accuracy	0.57	0.56	0.54
	95% random interval	[0.22; 0.43]	[0.21; 0.45]	[0.23; 0.42]
δ = 20 ms, Δ = 75 ms, Δ – δ/3 = 68.3 ms	Accuracy	0.61	0.67	0.63
	95% random interval	[0.21; 0.45]	[0.19; 0.47]	[0.20; 0.48]

The table also includes the estimated 95% interval of accuracies that can be expected due to chance. An accuracy of 1.00 implies that all validation observations have been correctly classified; an accuracy of 0.00 implies instead that none have been correctly classified.

Figure 4 reports MRI histology results. Unlike the WT, the PDX features widespread infiltration of smaller cells (likely leukocytes) in between larger hepatocytes. This leads to between‐sample contrast ( higher in WT than PDX), replicated in PolyMap and SigFit . The value of is lower in PDX than WT in PolyMap, whereas no differences are seen for SigFit. The SigFit metrics feature salt‐and‐pepper variations and are less smooth than PolyMap. Distributions (Supporting Information Table S1) confirm that agrees well with for both PolyMap and SigFit in PDX. In WT, is larger than , especially for PolyMap. The value of SigFit is more variable than PolyMap in both specimens. Supporting Information Figure S9 reports signal predictions from fitted parameters for both WT and PDX livers, highlighting that both SigFit and PolyMap provide good quality of fit. Supporting Information Figure S10 reports alternative SigFit cell‐size estimates obtained when is fixed to a specific value for all voxels, and not fitted. The value of obtained at fixed is highly dependent on the value used for : For some specific values, the between‐sample cell‐size contrast is even reversed, with larger in the PDX than in the WT, a finding that disagrees with histology.

FIGURE 4

Estimation of intrinsic cell diffusivity and cell size from the 9.4T ex vivo MRI scans of fixed mouse livers, with co‐localized hematoxylin and eosin (HE) histology. Top: Wild‐type (WT) case. Bottom: Patient‐derived xenograft (PDX) case (subcutaneous implantation of bone biopsy from metastatic prostate cancer). From left to right: Non‐DW image; co‐localized HE; example of cell segmentation on HE; histology‐derived cell size index ; MRI cell size estimates through PolyMap and SigFit estimation; MRI cell diffusivity estimates through PolyMap and SigFit estimation

DISCUSSION

Summary and key findings

We performed simulations to relate DW signal features (i.e., apparent diffusion/kurtosis coefficients, and ) to cell microstructure (cell diffusivity/size, and ) at fixed diffusion time, under the hypothesis of negligible sensitivity to extracellular/extravascular water and exchange. We also used cell‐size mappings learned from simulations on 9.4T ex vivo MRI of fixed mouse livers, comparing results to histology. Our work is motived by the fact that estimating summary cell‐size contrasts with minimal protocols may be useful in hospital settings, where scan time is limited and the latest technologies are not available.

Our main finding is that and offer sensitivity to and even when computed at realistic SNR levels, so it appears feasible to establish a mapping . Although the mapping does not estimate accurately and for the studied range, it captures salient cell‐size contrasts at fixed diffusion time. On the 9.4T MRI data, mappings provide cell‐size contrasts that are confirmed by histology, but overestimate , especially for larger cells. The overestimation of is less strong when this is estimated through biophysical models of restricted diffusion, which were considered as a potential alternative to mappings, at the price of more variable parametric maps (especially ).

Simulations

We used state‐of‐the‐art Monte Carlo simulations to study DW MRI protocols that could be implemented in the clinic (i.e., intermediate b‐values, fixed diffusion time, short scan time). Our results demonstrate that associations between and from such protocols and cell diffusivity and size , exist. The relationship is complex and non‐monotonic, with relatively small changes in and causing large variations of and . This may imply that biophysical liver models may benefit from intra‐compartmental kurtosis in the hepatocyte compartment, to better capture departures from Gaussian diffusion.

We also used paired and to compute polynomial functions that estimate and from and (PolyMap). Such functions offer sensitivity to the underlying and , even when computed on noisy data (SNR = 20). Although the estimates are not accurate for the smallest and largest values of and , they may suffice to characterise large cell‐size variations, such as distinguishing discrete cell‐size contrasts, as demonstrated through multinomial logistic regression. For reference, PolyMap was compared with fitting and based on a biophysical model of the DW signal (SigFit approach). Results from SigFit are in line with those from PolyMap. Although SigFit enables slightly more accurate estimation than PolyMap for the low‐intermediate values , PolyMap estimates are closer to ground‐truth values for of the order of 40 µm to 50 µm, plausible in pathological processes such as hepatocyte ballooning. Moreover, PolyMap exhibits higher precision (smoother and ) and better resolves . These results suggest that the relative performances of PolyMap and SigFit depend on the diffusion regime, and, more importantly, that overall it may be feasible to obtain summary descriptors of cell size from clinical acquisitions at fixed diffusion time, if analyzed with appropriate techniques. Such an approach could have application in high‐risk populations, such as patients with a history of hepatitis (at risk of hepatocellular carcinoma ) or primary colorectal cancer (at risk of liver metastases ), and in contexts in which implementing rich acquisitions is not possible. Moreover, mappings tuned for specific diffusion encodings may help mitigate interscanner variability. In that respect, they may prove useful in the retrospective analysis of multicenter clinical data featuring a variety of acquisition protocols. Nonetheless, we remark that acquiring prospective data at varying diffusion weightings and times should be the preferred way to perform cell‐size estimation, when possible.

Magnetic resonance imaging and histology

We estimated and on 9.4T ex vivo DW images of two formalin‐fixed mouse livers (one WT, one PDX), acquired at fixed diffusion time ( = 30 ms; = 10 ms). The MRI indices were related to co‐localized histological cell size , confirming results from simulations. Both PolyMap and SigFit provide good signal quality of fits, suggesting that they both are good representations of the diffusion MRI signal. Moreover, they both detect diffuse cell size alterations in the PDX liver that are confirmed by , despite overestimating actual cell‐size values ( is consistently lower than from MRI). This finding agrees with the overestimation seen in simulations for ground‐truth sizes of up to 35 μm to 40 μm, and may also result, at least in part, from histological tissue shrinkage and biases from neglected extracellular/extravascular water. Nonetheless, we acknowledge that the overestimation of is higher for PolyMap than for SigFit, especially for the WT liver. For PolyMap, we used a mapping evaluated at a single, fixed SNR. It is possible that more accurate results could be obtained learning a mapping for each voxel, tailored to spatially variant noise, or using more sophisticated mapping strategies beyond polynomial fitting (e.g., random forests). These are likely to outperform PolyMap, while also providing clearer biological interpretations than polynomial expansions, whose optimal degree is challenging to determine.

Notably, PolyMap detects PDX‐WT differences in , unlike SigFit. While it is challenging to verify this on the type of histological data at hand (routine hematoxylin and eosin staining), we speculate that it is possible that some of these between‐sample differences may exist. Supporting Information Figure S11 provides examples of the strikingly different cellular composition characterizing the two livers. On visual inspection, hepatocytes in the WT liver contain more fat than those in the PDX. Moreover, the PDX liver is characterized by a nonspecific, lymphoma‐like process, in which cells that are much smaller than normal hepatocytes invade vascular and extravascular spaces. Such cells may feature a distinct intracellular microenvironment as compared with normal hepatocytes, resulting in per‐cell heterogeneity. Taken as a whole, these findings suggest that differences in terms of intrinsic intracellular diffusivity between the two specimens cannot be ruled out a priori. In future work, we aim to perform richer immunohistochemical analyses to gain insight into the tortuosity of the intracellular space, and thus derive histological counterparts of to confirm our MRI findings.

Regarding our 9.4T diffusion MRI acquisition, we used a maximum b‐value of 4500 s/mm2. This is considerably higher than in simulations, where it never exceeds 2000 s/mm2, as in some clinical studies. This can be justified by considering that reductions of up to three times of the average apparent diffusion coefficient can be expected when scanning fixed liver tissue, as compared with the in vivo case. Therefore, = 4500 s/mm2 is expected to cause a signal attenuation somewhat comparable to approximately 1500 s/mm2 in vivo. Also, on the ex vivo data we perform PolyMap and SigFit analyses using a minimum b‐value of 1700 s/mm2. This is done to suppress partial‐volume effects with vessels and capillaries, which are filled at least in part with phosphate‐buffered saline. The diffusivity of phosphate‐buffered saline (1.8–2.0 ) is at least 8–10 times lower than the pseudo‐diffusion coefficient of the IVIM water pool in vivo (15–60 ), justifying the use of a minimum b‐value of 1700 s/mm2 against 100 s/mm2 as done in simulations.

We acknowledge that in this study we tested whether mappings learned on simulated MRI signals could be deployed on actual MRI measurements, performed on fixed ex vivo tissue at 9.4 T. In future work we aim to test such mappings on actual clinical MRI scans of the human liver, and investigate the performance of the approach in the presence of lower SNR, motion, and perfusion.

Methodological considerations

We used a simple geometric model based on perturbations of regular prisms to capture restricted diffusion. Although it sufficed to introduce variability in cell shape and to avoid overly simplistic representations (e.g., cubes), different models (e.g., meshes from histological images) could have been used. We plan to explore them in future work.

Another aspect is that our simulations focused on hepatocytes. We included heterogeneity in cell size/diffusivity, and accounted for partial volume with incoherent perfusion, , , , effectively relying on a two‐compartment model, under the hypothesis that the sensitivity to extracellular, extravascular water and transcytolemmal exchange are negligible. This assumption may be reasonable in the healthy liver, as hepatocytes are tightly packed within hepatic lobules, and account for 70%–85% of the liver volume. They are surrounded by networks of fluid‐filled conduits (sinusoidal capillaries, whose walls embed endothelial, stellate, dendritic, and Kupffer cells; and bile ducts ), whose signal fraction is expected to be on the order of 10%–20%. Interestingly, this two‐compartment model may capture the essence of the DW signal even in some pathological tissues, such as metastases. Nonetheless, extracellular, extravascular water may be relevant in the presence of other pathological processes, such as in fibrosis. In those cases, an additional compartment may be needed : While may still retain sensitivity to , they would not be specific. Finally, we neglected transcytolemmal water exchange. Known intracellular water residence times for hepatocytes and cancer cells of [40 ms; 150 ms] , imply that neglecting exchange may be reasonable in the short/intermediate diffusion times considered here. Nonetheless, further biases may be expected for longer diffusion times. Our work represents a first exploratory characterization of the main components of the liver parenchyma and in specific measurement conditions. In future work, we will generalize our analysis to more complex tissue models.

We explored relationships between and cell microstructure , and tested whether information derived from Monte Carlo simulations enables a mapping . Linking cumulants to microstructure is a powerful approach that has shown promise in the brain. , , Nonetheless, depend strongly on the diffusion‐encoding protocol used for acquisition. Therefore, one would need to learn a mapping for the specific diffusion protocol at hand (i.e., δ, Δ and b‐values). Moreover, may be difficult to measure accurately on noisy data (e.g., can be unstable when is low, being computed by dividing the second cumulant by ). In the future, more advanced signal‐to‐microstructure mappings will be explored (e.g., machine learning , , ).

Moreover, we limited our analysis to clinical single diffusion encoding with moderate b‐values. We acknowledge that more advanced encodings may provide more accurate cell‐size figures, such as combining pulsed/oscillating gradients, b‐tensor encoding, , and power law modeling. In particular, fitting biophysical models of restricted diffusion on measurements performed at varying diffusion time is likely to outperform cell‐size estimation at fixed diffusion time. However, we note that considering such protocols goes beyond the scope of this paper: Our main focus is on simple diffusion encodings at a fixed diffusion time. Our results quantify how much information on cell size can be retrieved with such minimal schemes, being these routine in hospital settings. However, when cell‐size estimation is sought in prospective studies, we recommend that diffusion protocols probe multiple diffusion times—scan time and hardware allowing.

We compared mappings (PolyMap) against fitting a biophysical model of intracellular restricted diffusion (SigFit) on protocols including a single diffusion time. We acknowledge that analyses such as SigFit would normally be performed on measurements performed at variable diffusion times, , , given the challenge of resolving and . A common way to reduce the number of tissue parameter unknowns in such model‐based approaches is to fix to a specific value across all voxels, and estimate only . While this would likely stabilize the fitting, it may lead to unphysical solutions if inappropriate values are used for . This is demonstrated in Supporting Information Figure S10 (the PDX‐WT cell‐size contrast can even be reversed depending on ), warning against the risks of using overly simplified analytical models in conjunction with minimal diffusion encodings.

Finally, our simulated MRI protocols were based on averaging over three gradient directions, common in liver MRI, , and included seven nonzero b‐values, corresponding to a tolerable 5/10‐minute scan. Additional analyses (Supporting Information Figure S12) show that three‐direction averaging suffices to account for anisotropy, and provides that are consistent with mean diffusivity/kurtosis from tensor fits on richer directional schemes. Supporting Information Figure S13 suggests that using seven nonzero b‐values may be a reasonable compromise between accuracy/precision and scan time.

CONCLUSIONS

In experimental conditions for which extracellular, extravascular signal sources and transcytolemmal exchange can be neglected, salient but potentially relevant information on liver cell size and diffusivity may be retrieved from simple diffusion encodings at a fixed diffusion time, provided that these are analyzed with appropriate computational techniques.

DISCLOSURES

FG was supported by PREdICT, a study funded by AstraZeneca at the Vall d’Hebron Institute of Oncology (Barcelona, Spain). The authors report no conflicts. AstraZeneca was not involved in any aspect concerning this study; it has not influenced the analysis of the data, the interpretation of the results, and the decision to submit the manuscript for publication.

FIGURE S1 Synthetic hepatocytes used in this study for Monte Carlo simulations (water diffusion was simulated within such synthetic cells), obtained by perturbing the position of the vertices of triangularly meshed regular prisms. Top to bottom: Different shapes of the prism bases (square, pentagonal, hexagonal). Left to right: Different unique perturbations. The figure also illustrates , which is equal to the base‐to‐base height as well as the diameter of the circumcircle relative to each base

FIGURE S2 Scatter plots of color‐coded by the underlying average intrinsic cell diffusivity (A–E) and cell size (F–J), as obtained when no noise is added to the synthetic MRI signals. From left to right: Different diffusion times (δ/Δ = 20/25 ms in [A] and [F]; δ/Δ = 40/50 ms in [B] and [G]; δ/Δ = 20/50 ms in [C] and [H]; δ/Δ = 10/50 ms in [D] and [I]; and δ/Δ = 20/75 ms in [E] and [J]). The figure refers to a minimum/maximum protocol b‐value of = 100/2000 s/mm2. Noise‐free intracellular diffusion‐weighted (DW) signals are contaminated by intravoxel incoherent motion (IVIM)–like partial volume

FIGURE S3 Scatter plots of color‐coded by the underlying average intrinsic cell diffusivity (A–E) and cell size (F–J), as obtained when noise is added to the synthetic MRI signals at an SNR at b = 0 of 20. From left to right: Different diffusion times (δ/Δ = 20/25 ms in [A] and [F]; δ/Δ = 40/50 ms in [B] and [G]; δ/Δ = 20/50 ms in [C] and [H]; δ/Δ = 10/50 ms in [D] and [I]; and δ/Δ = 20/75 ms in [E] and [J]). The figure refers to a minimum/maximum protocol b‐value of = 100/1000 s/mm2. Noise‐free intracellular DW signals are contaminated by IVIM‐like partial volume

FIGURE S4 Examples of predictions of intrinsic cell diffusivity and cell size on the validation set. (A–C) Scatter plots showing prediction of average intrinsic cell diffusivity . (D–F) Scatter plots showing prediction of average cell size . Left: Signal cumulants at fixed diffusion time colored by underlying ground‐truth and . Middle: Signal cumulants at fixed diffusion time colored by predictions of and as obtained with the PolyMap approach, which relies on using smooth functions and from polynomial interpolation. Right: Signal cumulants at fixed diffusion time colored by predictions of and as obtained with the SigFit approach, which relies on the estimation of and via routine nonlinear least‐squares fitting on the MRI signal. The figure refers to the case when the minimum/maximum protocol b‐values are equal to = 100/2000 s/mm2 and the diffusion gradient duration/separation is δ = 20 ms/Δ = 75 ms for SNR (no noise injected to the DW measurements) and in presence of IVIM contamination

FIGURE S5 Examples of predictions of intrinsic cell diffusivity and cell size on the validation set. (A–C) Scatter plots showing prediction of average intrinsic cell diffusivity . (D–F) Scatter plots showing prediction of average cell size . Left: Signal cumulants at fixed diffusion time colored by underlying ground truth and . Middle: Signal cumulants at fixed diffusion time colored by predictions of and as obtained with the PolyMap approach, which relies on using smooth functions and from polynomial interpolation. Right: Signal cumulants at fixed diffusion time colored by predictions of and as obtained with the SigFit approach, which relies on the estimation of and via routine nonlinear least‐squares fitting on the MRI signal. The figure refers to the case when the minimum/maximum protocol b‐values are equal to = 100/1000 s/mm2 and the diffusion gradient duration/separation is δ = 20 ms/Δ = 25 ms for SNR = 20 and in the presence of IVIM contamination

FIGURE S6 SigFit prediction errors for and scattered against ground‐truth values of and for different diffusion gradient timings at a fixed maximum b‐value of 2000 s/mm2 and SNR = 20. From left to right: Different gradient timings (δ/Δ = 20/25 ms in [A] and [F]; δ/Δ = 40/50 ms in [B] and [G]; δ/Δ = 20/50 ms in [C] and [H]; δ/Δ = 10/50 ms in [D] and [I]; and δ/Δ = 20/75 ms in [E] and [J]). Plots on top (A–E) refer to ; plots on the bottom (F–J) refer to . For each fixed value of (on top, or on the bottom), median errors with interquartile ranges for varying (on top, or varying on the bottom) are also reported

FIGURE S7 PolyMap and SigFit prediction errors for and scattered against ground‐truth values of and , for different diffusion gradient timings at a fixed maximum b‐value of 2000 s/mm2 and SNR (no noise injected to the data). Top (A‐J, rows one and two): PolyMap results (estimation from cumulants via smooth polynomial functions), with errors on row one and errors on row two. Bottom (K‐T, rows three and four): SigFit results (direct fitting on the MRI signal), with errors on row three and errors on row four. From left to right: Different gradient timings (δ/Δ = 20/25 ms in [A] and [F]; δ/Δ = 40/50 ms in [B] and [G]; δ/Δ = 20/50 ms in [C] and [H]; δ/Δ = 10/50 ms in [D] and [I]; and δ/Δ = 20/75 ms in [E] and [J]). In rows one and three, median errors with interquartile ranges for varying and fixed are reported. In rows two and four, median errors with interquartile ranges for varying and fixed are reported

FIGURE S8 PolyMap and SigFit prediction errors for and scattered against ground‐truth values of and for different diffusion gradient timings at a fixed maximum b‐value of 2000 s/mm2 and SNR = 20. Top (A‐J, rows one and two): PolyMap results (estimation from cumulants via smooth polynomial functions), with errors on row one and errors on row two. Bottom (K‐T, rows three and four): SigFit results (direct fitting on the MRI signal), with errors on row three and errors on row four. From left to right: Different gradient timings (δ/Δ = 20/25 ms in [A] and [F]; δ/Δ = 40/50 ms in [B] and [G]; δ/Δ = 20/50 ms in [C] and [H]; δ/Δ = 10/50 ms in [D] and [I]; and δ/Δ = 20/75 ms in [E] and [J]). In rows one and three, median errors with interquartile ranges for varying and fixed are reported. In rows two and four, median errors with interquartile ranges for varying and fixed are reported

FIGURE S9 Examples of DW images obtained ex vivo on the two fixed mouse livers. (A) Images from the wild‐type (WT) liver, alongside image predictions based on fitted model parameters for PolyMap and SigFit. (B) Similar information as in (A) but for the patient‐derived xenograft (PDX) liver. (C,D) Examples of MRI measurements (i.e., logarithm of measured signals) from one representative voxel alongside PolyMap and SigFit fittings for the WT (C, left) and PDX (D, right) livers

FIGURE S10 SigFit cell size map in the two fixed liver samples scanned at 9.4 T: WT (top) and PDX (bottom). From left to right: Full SigFit estimation (cell size and cell diffusivity are estimated jointly at fixed diffusion time); estimation when is fixed and not estimated (values used for : 0.5, 0.75, 1.0, 1.25, 1.25, and 1.5 μm2/ms, as shown from left to right). Median values of across the entire samples are reported for each specimen and SigFit configuration

FIGURE S11 Image patches illustrating the different microstructural environments observed in the two fixed mouse livers studied in this paper. Top: Patches from the WT liver, showing healthy hepatocytes surrounded by stellate cells and sinusoidal capillaries. Bottom: Patches from the PDX liver. On visual inspection, hepatocytes in the PDX appear to contain less fat than in the WT. Moreover, the PDX liver is characterized by a nonspecific, lymphoma‐like process, in which cells that are much smaller than hepatocytes invade vascular and extravascular spaces

FIGURE S12 Investigation on the impact of the number of gradient directions used to compute directionally averaged signals. The figure shows results obtained for fitting performed on seven nonzero b‐values in the range ; δ = 20 ms, Δ = 75 ms; intrinsic cell diffusivity and cell size ; {3, 9, 21, 30, 61} isotropically distributed gradient directions for each nonzero b‐value; no IVIM contamination. (A) Directionally averaged DW signals obtained at SNR for {3, 9, 21, 30, 61} gradient directions. (B,D) Distribution of apparent diffusion coefficient and apparent kurtosis coefficient calculated by fitting Equation 3 to directionally averaged signals over 1000 random noise instantiations with 20 SNR 100, for {3, 9, 21, 30, 61} gradient directions per b‐value. C,E, Distribution of diffusion tensor mean diffusivity and kurtosis tensor mean kurtosis calculated by fitting a full diffusion kurtosis tensor representation to all measurements over 1000 random noise instantiations with 20 SNR 100, for {3, 9, 21, 30, 61} gradient directions per b‐value. For full kurtosis tensor fitting we used the freely available DiPy package (https://dipy.org/documentation/1.4.1./examples_built/reconst_dki/#example‐reconst‐dki). Gradient directions were generated according to Cauryer et al (Magn Res Med. 2013; doi: 10.1002/mrm.24736; free download from http://www.emmanuelcaruyer.com/q‐space‐sampling.php)

FIGURE S13 Investigation on the impact of the number of b‐values used for apparent diffusion coefficient and apparent kurtosis coefficient computation from Equation 3. The figure shows results obtained for {19, 10, 7, 4, 3} nonzero b‐values and δ = 20 ms, Δ = 75 ms; maximum b‐values of 1000 s/mm2 and 2000 s/mm2; ; no IVIM contamination. (A,C) Distribution of ADC over 1000 random noise instantiations with 20 SNR 100, {19, 10, 7, 4, 3} nonzero b‐values, maximum b‐value of 1000 s/mm2 (A), and 2000 s/mm2 (C). Each plot also indicates the value of obtained with 19 nonzero b‐values and SNR for reference. (B,D) Distribution of apparent kurtosis coefficient over 1000 random noise instantiations with 20 SNR 100, {19, 10, 7, 4, 3} nonzero b‐values, maximum b‐value of 1000 s/mm2 (B), and 2000 s/mm2 (D). Each plot also indicates the value of obtained with 19 nonzero b‐values and SNR for reference

TABLE S1 Median and interquartile range (within brackets) of histology‐derived cell‐size index and MRI metrics from the PolyMap and SigFit estimation procedures investigated in this work (average cell size and average intracellular diffusivity )

Click here for additional data file.