Microstructural models for diffusion MRI in breast cancer and surrounding stroma: an ex vivo study

The diffusion signal in breast tissue has primarily been modelled using apparent diffusion coefficient (ADC), intravoxel incoherent motion (IVIM) and diffusion tensor (DT) models, which may be too simplistic to describe the underlying tissue microstructure. Formalin‐fixed breast cancer samples were scanned using a wide range of gradient strengths, durations, separations and orientations. A variety of one‐ and two‐compartment models were tested to determine which best described the data. Models with restricted diffusion components and anisotropy were selected in most cancerous regions and there were no regions in which conventional ADC or DT models were selected. Maps of ADC generally related to cellularity on histology, but maps of parameters from more complex models suggest that both overall cell volume fraction and individual cell size can contribute to the diffusion signal, affecting the specificity of ADC to the tissue microstructure. The areas of coherence in diffusion anisotropy images were small, approximately 1 mm, but the orientation corresponded to stromal orientation patterns on histology.


| INTRODUCTION
Breast cancer screening allows the early detection of cancerous lesions, but improved technology increases the likelihood of detecting small, slow-growing cancers that do not require aggressive treatment.
It is estimated that 10% of women who have mammographically detected cancers would not have required treatment in their lifetime ('overdiagnosis'), 1,2 and that post-surgical radiation treatment may not improve 5-year overall survival in some groups of women (a form of 'overtreatment'). 3 There is therefore a need for further tumour characterisation, i.e. beyond cancer detection, to identify patients in whom overdiagnosis and overtreatment are likely; magnetic resonance imaging (MRI) is particularly appealing because of its non-invasive nature and sensitivity to microstructure. Breast cancers show great variation in microstructure: higher grades tend to have increased cell density and a more disorganised structure; 4 immune cell infiltration and cell differentiation affect the distribution of cell types and sizes; and there are changes in the extracellular matrix related to invasion. 5,6 Diffusion MRI is sensitive to many microstructural features.
Diffusion tensor imaging (DTI) 7 and neurite orientation dispersion and density imaging, 8 for example, have produced maps of brain architecture. In the context of cancer, methods such as VERDICT 9,10 and restriction spectrum imaging 11 estimate tissue features such as those related to cell density. However, the exploration of breast microstructure has been relatively limited in comparison.
Clinical work has focused on data acquisition at a small number of low (≤1000 s/mm 2 ) b values and mono-exponential fitting to give an apparent diffusion coefficient (ADC). The ADC shows a difference between benign and malignant lesions, [12][13][14][15] but there is a large overlap in the ADC values of the two groups in many studies, with an area under the receiver operating curve ranging from 0.72 to 0.97 for ADC alone. 16,17 Attempts to use ADC to distinguish histological grades have produced mixed results, even with large numbers of patients. [18][19][20] Some variation in ADC with molecular subtype has been observed, [20][21][22] although results may be affected by the inclusion of necrotic regions common in triple-negative cancers.
However, ADC assumes that all of the water in a particular voxel can be represented by a single ADC. In reality, intracellular water is at least partly restricted by the cell membrane, extracellular diffusion depends on the extracellular space and organisation of cells, and anisotropic structures in the tissue produce diffusion orientation dependence. The use of an inappropriate model yields ADC values that are dependent not just on physiology, but on the choice of b value itself, 23 making comparisons between different scan protocols and between centres with different scanners difficult. The results from a given patient might also be difficult to interpret: although ADC correlates with cellularity in breast tumours, 24,25 immune responses and changes to the extracellular environment may also affect diffusion. More biologically motivated models have the potential to separate microstructural features related to cell proliferation from those caused by the immune response, invasion or less common cancer subtypes, and provide more specific features for tumour characterisation.
Recent clinical studies have begun to explore models beyond ADC. Intravoxel incoherent motion (IVIM) studies have looked at the vasculature in and around the tumour. 16,26,27 Diffusion kurtosis results suggest that diffusion is non-Gaussian and high-b-value measurements contain additional information. 28,29 Kurtosis can also be combined with IVIM. 30 Diffusion tensor (DT) modelling has produced inconsistent results, generally demonstrating lower fractional anisotropy (FA) in cancers compared with normal tissues, 31,32 but not always showing a distinction from benign lesions; 15,33,34 other anisotropy metrics may be more sensitive. 35 Furthermore, the source of anisotropy is uncertain: a preclinical breast cancer model observed lower FA in hypoxic regions with lower collagen fibre content, 36 but the differences between normoxic and hypoxic FAs were small, approximately 0.03.
Other groups have suggested partial restriction in the breast ducts, 31,35 although the average duct diameter (approximately 90 μm in normal breast and larger in patients with ductal carcinoma in situ 37 ) is much larger than the average distance travelled by water molecules during typical diffusion MRI experiments.
In this article, we examined the microstructure in a small set of cancer-containing, formalin-fixed breast tissue samples ex vivo. This allowed for high spatial resolution and histopathological comparison that might shed light on the source of diffusion signal differences. For example, a previous study found large differences in the ADC of epithelial cell regions compared with surrounding stroma, as well as qualitative differences in anisotropy. 38 The ex vivo approach also permitted longer scan times to obtain data over a broad range of gradient strengths, durations, orientations and diffusion times. This rich dataset was then fitted with a set of candidate models which describe the intracellular and extracellular spaces with different shapes and degrees of restriction. Model parameters were then compared with the histological features. This information can be used to optimise clinical scan protocols 39

| MRI data analysis
The one-and two-compartment models outlined in Table 2 were fitted to the data voxel-wise, excluding voxels that were predominantly fat or with non-mono-exponential T 2 decay (non-mono-exponentiality was defined as a main peak area comprising less than 90% of the total spectral area in the T 2 spectrum from non-negative least-squares analysis 41 ). Compartment shapes are described in detail in Panagiotaki et al. 42 and are summarised in the Appendix: a Ball describes unrestricted (free or hindered) isotropic diffusion; a Tensor describes anisotropic free diffusion (with diffusion coefficients D 1 , D 2 and D 3 in three orthogonal directions characterised by angles θ and ϕ for the primary diffusion direction and α describing the angle of the secondary diffusion direction in the perpendicular plane); a Zeppelin is a cylindrically symmetrical tensor; and a Sphere describes diffusion restricted isotropically by an impermeable membrane with radius R. In this nomenclature, the conventional ADC model is represented by a Ball and a biexponential fit by Ball-Ball. In addition to the diffusion, shape and orientation parameters for each compartment shown in parentheses in Table 2, two-compartment models have extracellular and intracellular volume fractions f E and f I , respectively, and all models include the equilibrium signal S 0 and the T 2 relaxation time constant as fitting parameters.
Data were fitted using an iterative maximum likelihood procedure that accounts for local minima and Rician noise. 42 The noise was derived from correction of the standard deviation of signal in an empty region 43 in each image. Parameters were constrained as follows: 0.01 < D < 3 mm 2 /s for all diffusion coefficients; f E + f I = 1; 0.1 < R < 20 μm; 0 < S 0 ; 0.001 < T 2 < 3 s. Model selection was performed using the Akaike information criterion (AIC), AIC = − 2 ln L + 2k, 44

| Histology and registration
After imaging, the samples were processed and embedded in paraffin wax, and sections approximately 3 μm thick were cut at every 100 μm through the block. The slides were stained with haematoxylin and eosin (H&E), and digitised using a C9600-01 NanoZoomer Digital Slide Scanner (Hamamatsu, Hamamatsu City, Japan) at 20× magnification (21 708 pixels/cm).
Histological images were stacked into a volume using two-dimensional pairwise registrations between adjacent slices based on a block-matching strategy. 46,47 The transformation model used was rigid body and the similarity measure was the correlation coefficient.
This pairwise registration aligns each slice with the subsequent slice and then concatenates transformations to generate a volume consisting of stacked slices registered with respect to a reference slice in the middle of the image stack. The slice separation was taken as 100 μm.
To register the stacked histology volume to the T 2 -weighted MRI, 9-15 manually selected corresponding landmarks were identified in each MRI/histology volume pair. This provided an approximate initial alignment as a result of the differing volume orientations. For the final registration, the volumes were resampled to isotropic voxels (0.5 mm for MRI and 0.1 mm for histology) to reduce orientational bias and the relevant regions (bright foreground voxels in MRI and non-zero voxels in histology) were selected to restrict the region over which the similarity measure was calculated to the internal tissue contrast. An intensity-based affine registration from ITK 48 was then performed using normalised mutual information as the similarity measure and a regular step, single-scale optimisation with MRI as the target volume. The shear TABLE 1 Diffusion-weighted imaging (DWI) scan parameters. The b values corresponding to the gradient strengths (G) at each gradient duration (δ) and separation (Δ)/TE. The number of averages is given in parentheses for cases in which more than one average was performed m)   40  80  120  160  200  240  280  320  360  400   10/18  3  9  37  83  148  232  334  455  594  752  928   30/45  3  30  120  269  478 748 (4) 1077 (8) Tensor-Sphere component of the affine registration will, to first order, correct for any residual cumulative stacking error of the histology slices. A sample of the three-dimensional histology stack in each of three orthogonal views is shown in Figure S1 with the diffusion MRI slice overlaid. As a result of the orientation of the MRI slice with respect to the histological slicing plane, only a portion of the histology slice corresponds to the diffusion image, which is indicated by outlines of tumour regions in subsequent figures. The transformation was also applied to the primary diffusion vectors to keep alignment consistent with the image orientation.

| RESULTS
Fitting quality, parameter reproducibility and posterior parameter distributions were first examined as part of the model selection process.
In a second section, model parameters were compared with histology, first examining parameters associated with compartment size and restriction, and then those associated with orientational structure.    Zeppelin-Sphere, Tensor-Sphere) for f I , D I , R, θ and ϕ, but differ for related diffusion coefficients (e.g. D 1 from Ball-Sphere is between D 1 and D 2 for Zeppelin-Sphere). Histograms from voxels with lower f I (see Figure S3) showed similar patterns, but with narrower θ and ϕ distributions, probably because of the larger extracellular signal.

| Fitting results and model selection
There was generally good agreement between the posterior distributions and parameter maps (see Figure S4) for the Zeppelin-Sphere and Tensor-Sphere models. Reproducibility ( Figure S5) was also similar, and so subsequent results are presented for the simpler Zeppelin-Sphere model.

|
Model parameters -f I and R Figure 4 shows H&E-stained histology in the top row, parametric maps from ADC (second row) and selected parameters from the Zeppelin- For the mucinous carcinoma (last column), the fitted R parameter hits the 20 μm maximum allowed by the fitting procedure in most voxels. This large R value is equivalent to unrestricted diffusion given FIGURE 4 Parametric maps for the apparent diffusion coefficient (ADC, second row) and selected parameters from the Zeppelin-Sphere model. Regions of low cellularity and correspondingly high ADC are outlined in red and typically correspond to lower f I in the Zeppelin-Sphere model. However, some regions (e.g. cyan outline) have high cellularity relative to their surroundings, but higher ADC, which may be explained by a larger cell radius, R. Additional samples are shown in Figure S6. Regions without colour were excluded either as fat (orange outline in first column) or as non-mono-exponential T 2 (orange outline in second column), which often corresponded to necrotic regions on histology the diffusion lengths probed in this experiment; thus, this finding is consistent with Figure 2 data that the Zeppelin-Ball model is a better choice in this sample.    Although the constraints of clinical scans limit the diffusion data that can be obtained, the results of the rich protocol suggest that there is valuable information that is not captured by most clinical protocols and diffusion models. For example, a clinical protocol with scan parameters producing signal sensitive to the observed R range of 6-9 μm could be designed. This approach has been successfully applied in prostate cancer to distinguish tumour from benign regions, 9 and may improve tumour characterisation in breast.

| Parameter values
Parameter values varied both across and within samples. The MCMC parameter distributions suggested that the fitted parameters were relatively stable, and comparison with histology further supported the hypothesis that parameter variations reflected true microstructural differences. This heterogeneity makes a simple summary of parameter FIGURE 7 Higher magnification histology of two regions (outlined in boxes) with higher fractional anisotropy (FA), demonstrating correspondence between the primary diffusion direction and directional patterns of the fibrous stroma FIGURE 8 Higher magnification of the grade 3 mucinous carcinoma showing correspondence between the primary diffusion direction and directional patterns in the fibrous stroma, even in the areas of low cellularity. Note that the colour scale and arrow length have been adjusted from previous figures to better identify regions of coherence values challenging. The mucinous carcinoma examined had the highest ADC ( Figure 4, third column), as has been reported previously, and is attributed to low cellularity. 13,15 In other samples, the ADC in regions of low cellularity was approx- This is the first study using two-compartment restricted diffusion models to examine ex vivo breast tissue. However, Lasič et al. 50 examined MCF-7 cells that had been grown in vitro and found a median size of 13.2 μm assuming a log-normal distribution and a width of 0.6 μm. This is larger than the radius range observed in most sample regions in the present study (6-9 μm). This may reflect true biological differencesthe ex vivo samples include stromal regions containing smaller cells such as lymphoid cellsor may indicate cell shrinkage as a result of fixation or the use of a single cell radius parameter rather than a cell size distribution. Models incorporating cell size distributions were beyond the scope of this study, but should be examined in the future.
The diffusion coefficients themselves are affected by the room temperature scan, e.g. they are lower than the intracellular (1.5 × 10 −3 mm 2 /s) and extracellular (2.8 × 10 −3 mm 2 /s) values observed for in vitro breast cell samples at 37°C. 50 Fixation may also affect diffusion through cross-linking and decreased water content, but work in prostate suggests that the relative signal fractions in compartments are similar before and after fixation and that changes are unlikely to affect model ranking. 51,52 Studies in brain 53   There is also evidence that collagen reorients in invasive tumours, 6 and is more strongly aligned in malignant samples relative to hyperplastic samples. 56 Of particular note in this study is that, although the mucinous carcinoma had lower anisotropy, small regions of coherence were still observable and corresponded to stromal orientation patterns, whereas mucinous carcinomas have proven difficult to distinguish from normal and benign tissue using conventional ADC methods. 13 Thus, the ability to detect small regions of anisotropy within and around tumours is a potentially valuable biomarker, and may become achievable in the near future with the use of reduced field-of-view sequences, double diffusion encoding sequences 57 or higher order diffusion tensor methods. 17 FIGURE 9 The colour fractional anisotropy (FA) map for the Zeppelin portion of the Zeppelin-Sphere model for the original 0.25 × 0.25 × 0.5 mm 3 data (a) downsampled (by averaging) to 1.0 × 1.0 × 0.5 mm 3 (b) and 2.0 × 2.0 × 0.5 mm 3 (c). At lower resolutions, the anisotropy becomes weaker (colours are less bright) and it is more difficult to discern a coherent direction in the data

| Limitations
In addition to the use of a single average cell size parameter, and the fixation and temperature issues already discussed, this study has several limitations. The number of samples was small, but variation in the preferred microstructural model and parameters was observed even within this range of grades and histological subtypes of breast cancer. Samples were examined voxel-wise to maximise the information obtained about different microstructural environments, but additional samples are needed to determine whether the findings are generalisable across all breast cancers.
The gradient strengths used were larger than those commonly available clinically, but Figure 1 demonstrates that single-compartment models, such as DTI, diverge from the data even at low b values (e.g. red circles). More limited gradient strengths and diffusion times may result in more uncertainty in model parameters, but a priori information, such as that obtained from ex vivo studies and validated in vivo, may be useful in constraining models applied to more limited clinical data.
All models assumed no exchange of water between compartments during the measurement, although there may be some additional signal decay, particularly at long diffusion times and high b values, arising from exchange effects. T 2 was assumed to be mono-exponential, and a separate sequence ascertained where this assumption failed and excluded these voxels from fitting. The method could potentially be extended to include regions with multi-exponential T 2 , given sufficient data. We assumed spherical cells of uniform size, which is a simplification of the real biological system. In cases in which there is some eccentricity in the cell shape, the radius estimate will represent a volume average of this parameter. Future work could extend the model selection to include compartments with anisotropic restriction; however, fitting both a cell size and compartment eccentricity using a basic pulsed gradient spin echo sequence biases both the radius and eccentricity parameters. 58

| CONCLUSIONS
This is the first study to examine such a broad range of diffusion data in human breast tissue samples and to model the data using both anisotropy and restriction. The data from most cellular cancer regions and the adjacent fibroglandular tissue were best explained using a Tensor-Sphere or Zeppelin-Sphere model, indicating that both restriction and anisotropy are present in breast cancer tissues. There were no voxels in which ADC or DTI were the best models. Although variations in ADC often corresponded with variations in cellularity on histology, there were exceptions in which additional information was provided by the radius parameter R and intracellular volume fraction f I from the Zeppelin-Sphere model. Regions of anisotropy corresponded to extracellular regions with aligned collagen on histology, but directions were only coherent over areas of approximately 1 mm and require high spatial resolution or diffusion techniques sensitive to sub-voxel anisotropy 17,57 for their detection.

SUPPORTING INFORMATION
Additional Supporting Information may be found online in the supporting information tab for this article. The Ball compartment is equivalent to the signal for Gaussian diffusion: S Ball ¼ e −bD1 with b = (γgδ) 2 (Δ − δ/3), where γ is the gyromagnetic ratio, g is the gradient strength, δ is the gradient duration and Δ is the gradient separation.
The Zeppelin is the product of signal along the primary diffusion direction and the direction perpendicular to this: S Zeppelin ¼ e −b cos 2 ψD1 e −b 1−cos 2 ψ ð Þ D2 , where ψ is the angle between the normalised gradient direction ĝand the parallel diffusion direction n̂defined in spherical co-ordinates with θ and ϕ: cosψ ¼ ĝ⋅ n̂, n̂¼ sinθ cosϕ; sinθ sinϕ; cosθ ð Þ .
The Sphere signal is calculated using the Gaussian phase distribution approximation: 59 and J ν is the Bessel function of the first kind, order ν. The summation was carried out over the first 31 roots of the equation.