Evaluation of diffusion models in breast cancer

Purpose: The purpose of this study is to investigate whether the microvascular pseudodiffusion effects resulting with non-monoexponential behavior are present in breast cancer, taking into account tumor spatial heterogeneity. Additionally, methodological factors affecting the signal in low and high diffusion-sensitizing gradient ranges were explored in phantom studies. Methods: The effect of eddy currents and accuracy of b-value determination using a multiple b-value diffusion-weighted MR imaging sequence were investigated in test objects. Diffusion model selection and noise were then investigated in volunteers (n= 5) and breast tumor patients (n= 21) using the Bayesian information criterion. Results: 54.3% of lesion voxels were best fitted by a monoexponential, 26.2% by a stretchedexponential, and 19.5% by a biexponential intravoxel incoherent motion (IVIM) model. High correlation (0.92) was observed between diffusion coefficients calculated using monoand stretched-exponential models and moderate (0.59) between monoexponential and IVIM (medians: 0.96/0.84/0.72×10−3 mm2/s, respectively). Distortion due to eddy currents depended on the direction of the diffusion gradient and displacement varied between 1 and 6 mm for high b-value images. Shift in the apparent diffusion coefficient due to intrinsic field gradients was compensated for by averaging diffusion data obtained from opposite directions. Conclusions: Pseudodiffusion and intravoxel heterogeneity effects were not observed in approximately half of breast cancer and normal tissue voxels. This result indicates that stretched and IVIM models should be utilized in regional analysis rather than global tumor assessment. Cross terms between diffusion-sensitization gradients and other imaging or susceptibility-related gradients are relevant in clinical protocols, supporting the use of geometric averaging of diffusion-weighted images acquired with diffusion-sensitization gradients in opposite directions. C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1118/1.4927255]


INTRODUCTION
2][3][4][5] Highly proliferating malignant tumors result in cell density higher than normal parenchyma and therefore more restricted water diffusion.This leads to decreased values of ADC and can be used for tumor detection, as well as for monitoring and predicting tumor response. 4,5The signal attenuation, however, is not only a result of the random microscopic motion of water molecules influenced by cell density, membrane integrity, and tissue microstructure but it also depends on microperfusion and diffusion heterogeneity within a voxel.][10][11][12][13] Moreover, in the presence of additional "nondiffusion" field gradients, such as imaging or susceptibility-related intrinsic gradients, signal attenuation depends on the square of the sum of gradients involved.This may lead to a miscalculation of the diffusion-sensitization coefficient b.5][16] In addition, eddy currents 17 and lower SNR (Refs.18 and 19)  of high b-value images also compromise the quality of data fitting.
The purpose of this study was to investigate whether the microvascular pseudodiffusion effects resulting with nonmonoexponential behavior significantly affect the diffusionweighted signals in breast cancer taking into account tumor spatial heterogeneity.Additionally, methodological factors affecting the signal in low and high diffusion gradient ranges were explored in phantom studies.In particular, image distortion due to eddy currents, inaccuracy of b-value determination due to cross terms between diffusion-sensitizing and imaging gradients, and noise levels were considered.

METHODS
All imaging was performed on a 3 T MRI scanner (Philips Achieva, Best, Netherlands) using a dedicated seven-channel breast coil.

2.A. Multiple b-value DWI protocol
A diffusion-weighted protocol was implemented with ten b-values and diffusion-sensitization gradients applied in nine different directions.Pairs of diffusion-sensitization gradients parallel and antiparallel to imaging gradients directions and at 45 • to the imaging gradients directions were applied.The range of b-values employed was from 0 to 1150 s/mm 2 (Table I).All diffusion-weighted images were coregistered to the b = 0 image prior to quantitative analysis using a rigid body 3D transformation available from the manufacturer's software.Data were analyzed offline using in-house developed  software (Mathworks, Cambridge MA).

2.
A.1.a.Eddy currents.Image distortion was evaluated using a structured test object (5 mm rods in a 200 mm diameter cylindrical phantom filled with a solution of 10 mM CuSO 4 in distilled water) by comparing the b = 0 reference images with diffusion-weighted images acquired sequentially for three orthogonal diffusion gradient directions (frequency, phase encoding, and slice imaging).Three angulations (0˚, 15˚, and 30˚) with respect to the sagittal imaging plane were tested in order to investigate the influence of gradient system  20 phantoms (200 mm spheres containing distilled water with the addition of 10 mM CuSO 4 or 500 g/l sucrose, Sigma-Aldrich) were performed to assess signal attenuation and measured ADC values.A ROI encompassing 80% of both the water and sucrose spherical phantoms (central slice) was used to evaluate the calculated ADC values avoiding partial volume effects. 21ADC values were first calculated separately for opposite directions of the diffusion-sensitizing gradients and compared with the geometrically averaged ADC value.

2.A.2. Clinical studies
Twenty-six women, including 21 with known breast cancer (median age 52, range: 37-81 yr) and 5 healthy volunteers (median age 30, range: 28-33 yr) underwent a breast MRI examination with approval from the Research Ethics Committee and with written informed consent.Patients with biopsy-proven breast cancer had only the tumor-containing breast scanned using the DWI protocol.
In the case of volunteer examinations, both breasts were scanned using T 2 -weighted TSE VISTA (full breast coverage) and multiple b-value DWI sequences (three central slices).

2.C. Image processing and analysis
Regions of interest were drawn manually using postcontrast T 1 -weighted and high b-value images [DWI sequence, b = 1000 s/mm 2 , Fig. 1(B)] for tumor and on fat suppressed T 2 -weighted b = 0 EPI images for normal parenchyma.
Regions with ADC > 2 × 10 −3 mm 2 /s and an absence of contrast uptake on T 1 -weighted images (necrotic or cystic regions) were excluded from the ROIs. 22

2.D. Diffusion signal modeling
In this study, three diffusion attenuation models were considered.The first model assumes a pure diffusion process and is described by a monoexponential decay of the measured signal S 1 as a function of b-value, where S 0 is the initial signal acquired without diffusion weighting (b = 0) and ADC is the apparent diffusion coefficient.The second model is described by a stretched-exponential function. 6,7 in which DDC is the distributed diffusion coefficient and α is the heterogeneity index ranging from 0 to 1.In the case of homogeneous diffusion α = 1, and the function simplifies to the monoexponential decay described by the first model.Lower values of α result either from nonexponential behavior caused by the addition of proton pools with a range of diffusion rates within the imaged voxel or from a process where the motion is intermittent. 7The third model was the biexponential IVIM model, which accounts for a microcapillary perfusion and other flow effects contributing to signal attenuation at low b-values, This model is parameterized by the pseudodiffusion flow fraction f , the pseudodiffusion rate constant D * , and the diffusion rate constant D. [8][9][10][11] Median ROI values were calculated for ADC, DDC, and D.
The penalized-likelihood Bayesian information criterion 24 (BIC) was used to assess how well each of the three models is supported by the data (, Econometrics Toolbox), where L(θ) is the value of the maximized likelihood objective function for a model with k parameters and N data points.Given any two estimated models, the model with the lower value of BIC is the one to be preferred.BIC differences between the models are considered significant, strong, and very strong for ∆BIC = 2-6, 6-10, and >10, respectively. 25or each voxel, the model with the lowest BIC value was chosen and displayed in the form of model selection maps (Fig. 2) and used to calculate overall percentage model preferences.Correlation between measured ADC, DDC, and D was measured for all ROI voxels.Correlations between SNR and BIC values for each model were also used to assess if the signal to noise levels had an influence on model fitting.

2.E. Statistical analysis
The Wilcoxon signed-rank test 26 was used to compare malignant and healthy tissue diffusion parameters.Kendall's tau (τ) 27 was used to test for correlations.The null hypothesis was that there was no correlation between measured parameters.The strength of correlation was tested and the values were considered significant if P < 0.05.Statistical analysis was performed using the  Statistics Toolbox.

RESULTS
Four out of 26 subjects were excluded from the study because of fat suppression failure.

3.A. Image distortion
The eddy-current-related distortion, in the image phase encoding direction, was strongly dependent on the direction of the diffusion gradient (Table II) but not on the imaging plane angulation between 0 • and 30 • .The observed image stretch and displacement varied between 1 and 6 mm.
The sensitization gradient direction resulting in minimum image misregistration was along the phase-encode imaging gradient (sagittal: head-feet), corresponding to the z-axis of the scanner.The largest displacement observed for this direction in the b-value range between 1000 and 1150 s/mm 2 was 2 mm (see Table II).Only the data acquired with the diffusion-sensitization gradients along the phase encoding direction (superior-inferior) in the range of 0-800 s/mm 2 were used to model diffusion behavior in order to minimize the T II.The maximum image misregistration (mm) observed for different directions of the diffusion gradients in a sagittal breast DWI protocol.Frequency/phase/slice corresponds to y/z/x scanner coordinates in the sagittal orientation.

3.B. Accuracy and measured ADC values
The signal attenuation as a function of diffusion gradient b-value in phantoms is monoexponential (Fig. 3) and BIC S 1 << BIC S 2,3 .
ADC values measured with no geometric averaging depend on the direction of the diffusion gradient used [Fig.4(A)].

3.C. Noise
The SNR distribution calculated for breast ROIs at the bvalue 800 s/mm 2 was positively skewed with median values of 10.09 in tumor and 2.18 in the normal breast tissue.There was no correlation between SNR and the BIC values calculated for the models used (tumor: τ S 1 = 0.05, τ S 2 = −0.03,τ S 3 = 0.07 and parenchyma: τ S 1 = 0.03, τ S 2 = 0.08, τ S 3 = 0.11).

3.D. Model selection
A high correlation was observed between diffusion coefficients calculated using mono-and stretched-exponential models (τ tumor = 0.92, τ parenchyma = 0.9).A moderate correlation was found between the monoexponential model and IVIM (τ tumor = 0.59, τ parenchyma = 0.53).Median values were 0.96, 0.84, and 0.72 × 10 −3 mm 2 /s for ADC (S 1 ), DDC (S 2 ), and D (S 3 ), respectively, in tumor and 1.76, 1.75, and 1.24 ×10 −3 mm 2 /s in normal breast parenchyma (Fig. 5).All diffusion coefficients were significantly higher in the normal tissue than in malignant lesions (P S 1−3 < 0.001).The Bayesian information criteria from the lesion ROI data showed that 54.3% of voxels were best fitted by the monoexponential, 26.2% by the stretched-exponential, and 19.5% by the biexponential IVIM model (N voxels = 4601).In 83% of voxels, the BIC difference between the models was considered significant 25 with 66% positive, 9% strong, and 8% very strong evidence.Figure 2 shows an example of a model selection map and the corresponding alpha parameter map.The distribution of the stretched-exponential alpha parameter in the tumor is negatively skewed (γ 1 = −0.9)with a median value of α tumor = 0.75, the pseudodiffusion rate constant D * tumor = 21.5×10−3 mm 2 /s, and the pseudodiffusion flow fraction f tumor = 0.15.
In the normal parenchyma, 58.3% of voxels were best fitted by the monoexponential, 24.6% by the stretched-exponential, and 17.1% by the biexponential IVIM model (N voxels = 5022).In 83% of voxels, the BIC difference between the models was considered significant with 68% positive, 8% strong, and 7% very strong evidence.The distribution of the stretchedexponential alpha parameter in the parenchyma is negatively skewed (γ 1 = −0.92)with a median value of α parenchyma = 0.8, the pseudodiffusion rate constant D * parenchyma = 16.8 × 10 −3 mm 2 /s, and the pseudodiffusion flow fraction f parenchyma = 0.08.

DISCUSSION
Diffusion-weighted imaging is increasingly being used for cancer detection and diagnosis.However, it is still an open question as to which models are appropriate for DWI in various body regions and cancer types characterized by different microscopic environments.
Initial implementations of IVIM in the breast have been reported indicating the feasibility of perfusion fraction and pseudodiffusion coefficient measurements 13 in addition to the diffusion coefficient.A biexponential characterization of the diffusion attenuation signal has been reported for lesions, in contrast to a monoexponential dependence in normal tissue.These measurements were performed using the mean signal over ROIs and therefore disregard lesion heterogeneity, which is clearly observed in practice (Fig. 1).It has also been shown that the diffusion attenuation is no longer monoexponential at large b-values (b > 2000 s/mm 2 ) [28][29][30][31][32] including breast cancer preclinical studies, [33][34][35] and this effect is associated with restricted intracellular water diffusion rather than the properties of restricting boundaries in the extracellular space.Another confounding factor is fat signal, which depends on the choice of fat suppression technique, especially relevant in the case of higher b-value (>600 s/mm 2 ) breast DWI images. 36Although clinical systems are increasingly capable of delivering higher b-values, it is likely that noise will dominate breast images acquired with clinical parameters and b > 2000 s/mm 2 .Another concern is an assumption of the IVIM model, which states that time after diffusing particles change their direction is several times shorter than the total duration of diffusion gradients (i.e., 7 times-Le Bihan 9 ).A recent study shows that this assumption is not met in liver and pancreas, 16 suggesting that a further development is needed to appropriately calculate pseudodiffusion coefficient D * .There is, however, lack of evidence for similar shortcomings in breast parenchyma and lesions.
One of the methodological DWI caveats which is also a potential limitation of this work is a choice of b-values, which affects calculated diffusion coefficient values.The range of bvalues used should be adjusted for expected diffusion values in the examined tissues, observed SNR levels, and eddy currents.][13] In our work, the diffusion signal was analyzed for each pixel of ROIs defined in test objects, lesions, and normal tissue.Test object images acquired with diffusion-sensitization gradients of opposite polarity produced significantly different ADC values, suggesting that other gradients contribute toward diffusion weighting.This highlights an essential issue in modeling diffusion behavior: cross terms between diffusionsensitization gradients and other gradients were shown to be relevant (Fig. 4) and may induce error.In this work, geometrical averaging of images produced with diffusionsensitization gradients of opposite polarity was used to minimize the effect of other gradients in order not to compromise the accuracy of the modeling.In addition, phantom work showed that geometrical distortion of high b-value images is significant (Table II) and can introduce errors in the ADC evaluation.In this work, only the data acquired with the diffusion-sensitization gradients along the phase encoding direction (superior-inferior) in the 0-800 s/mm 2 b-value range were used to model diffusion behavior, minimizing the error.Other practical solutions to reduce eddy current distortion are to use the bipolar diffusion gradient scheme; 15 however, the use of temporally asymmetric diffusion gradients leads to the presence of significant concomitant fields and undesired signal voids. 15,16tatistical model selection analysis suggests that the monoexponential model describes majority of the data most accurately.There was, however, evidence of regions with signal best described by stretched-exponential and IVIM models in both tumor and normal breast tissue.This might suggest a presence of increased structural heterogeneity or microperfusion and might be of clinical use.The parameter alpha may be used to identify regions showing significantly non-monoexponential behavior [lower alpha, Figs.2(B) and 2(C)].Median diffusion coefficients calculated using all three models show significant differences between tumor and normal breast parenchyma.Data distribution representing intratumoral heterogeneity (Fig. 5) reveals greater overlap between malignant and normal tissues in the case of IVIM model.The median IVIM diffusion coefficients are markedly lower than the analogous parameters from other models.This could be a consequence of a signal fit bias for smaller b-values in the presence of pseudodiffusion and as a result higher DDC and ADC values. 37n important feature of the Bayesian information criterion is that more complex models are appropriately penalized relative to simpler models that fit the data equally well, given the data SNR. 24Correlations are not observed between the SNR and the penalized-likelihood BIC for the models considered.This indicates the preference for stretched-exponential model is unrelated to noise levels.In the case of insufficient SNR, a pure diffusion process characterized correctly by a monoexponential decay would appear artificially stretched in the higher b-value region due to a noise-floor effect, which would lead to a correlation between SNR and BIC values.

CONCLUSIONS
Cross terms between diffusion-sensitization gradients and other gradients (either imaging gradients or susceptibilityrelated gradients) are relevant in clinical protocols, supporting the use of geometric averaging of diffusion-weighted images acquired with diffusion-sensitization gradients in opposite directions.Microcirculatory pseudodiffusion and intravoxel heterogeneity effects were not observed in approximately half of normal and malignant breast tissue.This result indicates that stretched and IVIM models should be utilized in local voxel analysis rather than global tumor assessment.
Parametric maps F. 1. Example of MRI diffusion images [(A) b = 0 s/mm 2 , (B) b = 1000 s/mm 2 , and (C) ADC map] of 51 yr-old patient with invasive ductal carcinoma in the right breast. of ADC (Refs.1-5) and SNR (Ref.23) were calculated pixel by pixel for ROIs in the tumor and in healthy breast tissue (patients and volunteers) using .

F. 2 .
(A) An example of lesion ROI overlaid on an anatomical T 2 -weighted image.(B) Model selection map showing regions where monoexponential (S 1 ), stretched-exponential (S 2 ), and biexponential (S 3 ) models were best supported by the data (smallest Bayesian information criterion).(C) An alpha parameter map for the stretched-exponential model.A high alpha value close to unity shows a preference for the monoexponential model.(D)-(F) Examples of measured signal (log S) as a function of diffusion weighting (b-values) for voxels with different models selected using BIC [(D) monoexponential, (E) stretched-exponential, and (F) biexponential IVIM].Medical Physics, Vol.42, No. 8, August 2015

F. 3 .
The b-value dependence on normalized DWI signal for water and sucrose (500 g/l) at 21 • C. displacement of individual voxels with increasing b-values, as observed in the phantom work.
shows a clinical example of low [Fig.1(A)] and high [Fig.1(B)] b-value DWI-EPI images (multiple F. 4. Histograms of ADC values for diffusion-sensitization gradient in two opposite directions (plus and minus) measured for water (21 • C) (A) and a histogram for geometrically averaged data (B).Medical Physics, Vol.42, No. 8, August 2015 F. 5. A comparison of patient median ADC, DDC, and D in malignant and normal breast tissue.All diffusion coefficients were significantly higher in the normal tissue than in malignant lesions (P S 1−3 < 0.001).