1 INTRODUCTION

Diffusion-weighted MRI (DWI) shows promise for breast lesion classification without exogenous contrast,1 with the apparent diffusion coefficient (ADC) being clinically useful in distinguishing between malignant and benign breast lesions.2-4 However, the underlying assumption of free, or Gaussian, diffusion is a known simplification of the underlying diffusion process in tissues.5 Improved lesion classification through the use of more complex DWI signal representations has been reported. 6, 7

One such representation, the intravoxel incoherent motion (IVIM) model proposed by Le Bihan et al,8 explicitly attempts to account for perfusion influence; increased lesion perfusion is a known marker of more aggressive disease, and the associated IVIM parameter ƒ is a potential biomarker for breast lesion characterization and treatment response.9 The IVIM model describes the diffusion-weighted signal as the weighted sum of signal components originating from two non-exchanging compartments; one pseudo-diffusion and one “true” diffusion component, modelling the motion of water molecules in the blood vessels and the surrounding tissue, respectively. The IVIM model signal dependence on b value is given as (1)where S(0) is the signal intensity at b = 0 s/mm2 (b0), D describes the diffusion coefficient in tissue, D∗ describes the pseudo-diffusion coefficient within the blood vessels, and ƒ represents the volume fraction of the pseudo-diffusion component.

Despite much work with IVIM in different organs and pathologies, standardization of IVIM-centered acquisition protocols remains an unmet need.10-13 In particular, a dependence on imaging echo time (TE) for derived IVIM parameters arising from relaxation effects, which must be accounted for to ensure accurate estimation, has been demonstrated.14, 15 With TE often minimized (not standardized) for DWI acquisition, this potentially confounds studies using data acquired across different systems.

The standard IVIM model described above neglects any weighting from differential relaxation rates between compartments. Explicit inclusion of two distinct transverse relaxation variables gives an extended IVIM model: (2)where S(0, 0) is the signal intensity at b = 0 s/mm2 and  TE = 0 ms, ƒc is the corrected pseudo-diffusion signal fraction, and T2p and T2t are transverse relaxation times in the pseudo- and true diffusion compartments.

By explicitly modelling transverse relaxation times, TE-independent values for the pseudo-diffusion coefficient D∗ and signal fraction ƒ can be estimated. In the healthy liver at 1.5 T, Jerome et al found that the estimations of standard-model IVIM ƒ increased with longer TE, resulting in overestimations of more than ≈ 50% arising from longer transverse relaxation times in blood than in tissue.16

Transverse relaxation times in both blood and tissue depend on the magnetic field strength B0,17 where T2 generally decreases for increasing B0 and the dependence is stronger in blood than in tissue. At the current standard clinical field strength of 1.5 T, T2 in blood is considerably longer than in tissue, resulting in a significant effect on the estimated pseudo-diffusion fraction, with ƒ > ƒc. Conversely, at the ultra-high field of 7 T, T2 is significantly shorter for venous blood than for tissue, giving ƒ < ƒc and consequently a smaller IVIM effect. Values of ƒ and ƒc are therefore expected to be more similar at 3 T than at 1.5 T, but significant differences may still exist.

The effect of echo time dependence on IVIM measurements in the breast remains an open question. Inclusion and modelling the effect of distinct transverse relaxation times with IVIM analysis would allow for more general comparison of IVIM data across varying protocols and systems, albeit at the expense of additional scan time (and patient discomfort), and so should be critically assessed.

The aim of this study was to investigate the echo time influence on IVIM-derived parameters in breast lesions at 3 T. Simulated data and experimental DWI data acquired from breast cancer patients were examined to assess whether the effect results in significantly different ƒ and ƒc. Increased accuracy of the IVIM measurements could provide improved biomarkers for breast cancer characterization and treatment planning, and would enable comparisons of IVIM parameters measured at different echo times in multi-center studies.

2 METHODS 2.1 Simulations

Synthetic MRI data were generated using Equation (2) with scan parameters TE, TE, and b values matched to the values of the clinical DWI protocol described below, and with breast IVIM and T2 properties at 3 T taken from literature:  T2t = 71 ms, 18 ƒ = 9.8 %, D = 1.15 μm2/ms and D* = 15.1 μm2/ms.19 The value of T2p depends heavily on oxygenation level; we therefore selected two values: T2p= 100 ms and T2p = 150 ms, which represent oxygenation levels of 80 % and 95 % respectively.20 To simulate noise, random numbers were added to the underlying signal value, using a Gaussian distribution with mean value μ = 0 and standard deviation σ = k•S(TE = 77 ms, b = 0 s/mm2), with the scale factor k determining the level of noise and thus the signal-to-noise ratio (SNR). The following values of k were used in the calculation: 0.01, 0.02 and 0.04, corresponding to SNR values of 100, 50 and 25 for the b0 image at TE = 77 ms, respectively. The procedure was repeated 500 times for each value of k, thereby simulating the signal from 500 repeated measurements, with the parameter estimation performed on each synthetic measurement in order to evaluate the effect of measurement noise on the uncertainty in the resulting model parameter values.

2.2 In vivo experiments 2.2.1 Patient cohort

Patients included in this study were recruited as part of a multi-center prospective phase II trial (PETREMAC, Clinicaltrials.gov #NCT02624973) approved by the regional committee for medical and health research ethics (REC) in western Norway (2015/1493); all patients gave informed consent prior to inclusion. A total of 24 patients (median age 53 years, range 37 to 74), with biopsy-proven invasive breast cancers (4 cm inclusion criterion) or locally advanced breast cancer, were scanned with an extended MRI protocol at one of the sites. Patients subsequently underwent individualized neoadjuvant therapy, determined by tumor estrogen receptor (ER), progesterone receptor (PgR), human epidermal growth factor-2 (HER2), and TP53 mutation status. Data from these patients were previously reported in a study that focused on the segmentation of breast lesions using simultaneously-acquired 18F-FDG-PET (fluorodeoxyglucose positron emission tomography) data, where the monoexponential diffusion model was applied to a subset of the DWI acquired.21

In the current study we report results from baseline scans (prior to treatment) in a subset consisting of 15 patients for whom successful DWI data including acquisition with the short echo times, described below, were collected. Nine of the 24 original patients were excluded due to DWI image artefacts (n = 7) and incomplete data acquisition (n = 2). A summary of demographic characteristics of the patients included in this study is presented in Table 1.

TABLE 1. Clinical characteristics of the patient cohort Characteristics All patients Histological type NST 13 ILC 1 MD 1 Histological grade 2 7 3 8 Estrogen receptor status (%) Negative 5 ≥1-10% 1 > 10%-50 7 > 50% 2 Progesterone receptor status (%) Negative 10 ≥1-10% 3 >10%-50 2 HER2 status Positive 7 Negative 8 Ki67 (%) <30 5 ≥30 10 Pathological characteristics are determined based on histopathologic analysis of pre-treatment core needle biopsies. NST—invasive carcinoma of no special type; ILC—invasive lobular carcinoma; MD—carcinoma with medullary features; HER2—human epidermal growth factor receptor 2. 2.2.2 MRI acquisition and image preprocessing

The MRI protocol included a multiple-b-value DWI sequence for conventional analysis, using six b-values: 0, 50, 120, 200, 400, and 700 s/mm2. Axial bilateral images were acquired using a single-shot echo-planar imaging (EPI) sequence, with TR = 9000 ms, TE = 77 ms, GRAPPA 2, resolution 2x2x2.5 mm3, 60 slices. Additionally, a geometry-matched sequence with b-values b0 = 0 s/mm2 and b50 = 50 s/mm2 was acquired at the minimum echo time allowed with this smaller maximum b value, TE = 53 ms. Correction for geometric distortion for each diffusion-weighted image was performed using a phase-reversed b0 image (acquired at TE = 77 ms), as demonstrated by Teruel et al for use in the breast.22 Other MRI sequences acquired included clinical standard Dixon, T2-weighted, and dynamic contrast-enhanced (DCE) imaging (axial 3D FLASH, TR/TE 5.88/2.21 ms, 0.7x0.7x2.5 mm3, 72 slices, flip angle 15°, eight dynamic images (one baseline) with 1 min time resolution.

Volumes of interest (VOIs) for analysis were delineated on the enhancing solid tumor on an early post-contrast DCE image, ignoring satellite regions, and validated by an expert radiologist (A.Ø.). The VOIs were resampled using Elastix23 to exactly match the resolution of the DWI images, and thus give direct voxel-to-voxel correspondence; see figure 1.

Left: example subtracted DCE-MRI image with drawn ROI. Middle: ROI exported to DWI images. Right: mean ROI signal curve

2.3 DWI analysis

Apparent transverse relaxation times T2 were calculated using an exponential decay model for images acquired with different echo times (53 and 77 ms) with both b0 and b50 diffusion weighting.

IVIM parameters were estimated by fitting Equation (1) to the measured data. Since the contribution to the signal from blood flow can be assumed to be negligible for large b values (above 200 s/mm2), the fitting was performed using a segmented approach24, 25: D was first calculated from the monoexponential decay for high b values (b ≥ 200 s/mm2); D∗ and ƒ were then estimated by nonlinear fitting of the biexponential IVIM equation to the data for all b values measured at TE = 77 ms, using the previously estimated values for D.

To estimate the parameters for the TE-corrected model, the diffusion coefficients, D and D∗, were calculated using the same procedure as described for standard IVIM, and given as input to the corrected IVIM model as fixed parameters. A biexponential fit was then done to estimate the signal fraction ƒcand the transverse relaxation times, T2t and T2p. Preliminary analysis showed that when both T2 parameters were allowed to vary freely this resulted in unstable parameter estimations of T2 for blood that sometimes took very high values, yielding non-physical solutions. Due to this, the transverse relaxation time in blood, T2p, could not be included as a free parameter, and fitting of the model was performed using different fixed values of T2p (100 ms and 150 ms).

SNR was estimated using a procedure similar to the one in reference.15 First, to estimate the noise; in each patient, a single slice where a homogeneous region of the breast muscle was visible was chosen from the b0 image. Then a region of interest (ROI) covering this homogeneous region was carefully drawn and the standard deviation of the measured intensity within that ROI was calculated. The average value of this standard deviation across all patients was then used as an estimate of the measurement noise. Second, the mean signal from each lesion VOI was divided by standard deviation to estimate the single-voxel measurement SNR.

Within the lesion, additional sources of intensity variations, such as tissue heterogeneity and physiological noise, may result in a larger standard deviation than expected from the pure measurement noise. Therefore, the standard deviation within each lesion VOI was also calculated and compared with the above standard deviation.

2.4 Statistical analysis

The Wilcoxon signed rank test was used to calculate p values for the estimation of transverse relaxation times at different b values and parameters estimated using the two different models. p values for the simulated data analysis were calculated using the double-sided Student dependent t test for paired samples. For the intratumor variability analysis, 30 bootstrap samples of the VOI measurement were generated from the original samples (VOIs) of measured voxel signals using random sampling with replacement. Parameter estimations were done as described above for standard IVIM and corrected IVIM (fixing the transverse relaxation time of blood T2p at 120 ms) to the signal mean value calculated in each bootstrap sample and for each patient. The parameter estimations and statistical analysis were done in Python 3.8 and MATLAB (MathWorks, Natick, MA).

3 RESULTS 3.1 Simulations

For the simulated data analysis, parameter estimations were calculated using T2p = 100 ms and T2p = 150 ms to investigate how this affected the estimated corrected signal fraction, ƒc, and what range of values ƒc can take depending on the transverse relaxation time in blood, T2p. The results are shown in Table 2. Both the standard and the corrected model are very sensitive to measurement noise, and have a high variance already at an SNR of 100. The corrected IVIM model gives large variations in the estimates of ƒc already at an SNR of 25, which is a typical single-voxel SNR for DWI protocols and reflects the SNR measured from the patient data in our current study. Using the longer transverse relaxation time for the parameter estimation ( T2p = 150 ms), which gives a more distinct separation between the signal contributions from the two compartments, decreased the variation in ƒc, but only to a limited degree. Already at an SNR of 100, ƒ calculated using a blood T2 of 100 ms is within the 95% confidence region of ƒc, while this happens at an SNR of 50 when a blood T2 of 150 ms is used.

TABLE 2. Results from signal simulation experiments ƒ ƒc T2p = 100 ms T2p = 150 ms SNR→∞ 0.127 0.101**** 0.0795**** SNR = 100 Median 0.127 0.100**** 0.0790**** P5% 0.084 0.068 0.0541 P95% 0.167 0.133 0.1039 SNR = 50 Median 0.130 0.099**** 0.0788**** P5% 0.054 0.45 0.0359 P95% 0.208 0.165 0.1284 SNR = 25 Median 0.153 0.114 0.0901* P5% 0.049 0.039 0.0295 P95% 0.296 0.249 0.1917 Parameter estimations of the pseudo-diffusion signal fraction, using the standard IVIM model f and the model correcting for echo time dependences, , at different compartmental T2 relaxation times. Pi%—the ith percentile; 3.2 In vivo data

As a first step to explore the echo time dependence of the measured signal in patients, the transverse relaxation times were calculated at b values b0 and b50 using the signal sampled at the two echo times. The results are shown in Figure 2. The median of the calculated transverse relaxation time at b50 was significantly lower (Wilcoxon signed rank test, p = 0.01) compared with b0 (68.7 ms and 74.5 ms, respectively), representing a decrease of about 8% and demonstrating a change in apparent T2 related to the removal of fast-diffusing (eg perfusion) signal contribution.

T2 transverse relaxation times calculated at diffusion weighting b = 0 s/mm2 and b = 50 s/mm2

The average value of measurement noise, as estimated from the standard deviation value within homogeneous muscle tissue ROIs in the b0 image, was 24, while the mean signal intensity across all patients within the lesion VOIs at b0 for TE = 77 ms was 610. This yields a nominal single-voxel SNR of 25. The within-lesion VOI standard deviation in the same b0 images gave an average value of 263 across all patients, which is more than 10 times larger than expected from pure measurement noise. This points towards an additional contributing factor within lesions that will influence the performance of associated IVIM analysis.

Descriptive statistics for the estimated parameters for the standard IVIM model are given in Table 3. Similarly as for the synthetic data analysis, setting all parameters to vary freely resulted in very high estimations of the transverse relaxation time for blood and concurrent large variations in ƒc (results not shown). Using the same approach as described for the analysis of the synthetic data and keeping the transverse relaxation time for blood fixed, more stable results for the corrected fraction ƒc were obtained, with significantly lower estimates for ƒc than for ƒ calculated using the uncorrected IVIM model (p < 0.01). The resulting estimated parameters with a transverse relaxation time for blood of 100 and 150 ms are reported in Table 3, again showing reduced variation associated with a higher fixed T2p.

TABLE 3. Descriptive statistics for the parameters ƒ, D, and D∗ estimated with the standard IVIM model, and the parameters (ƒc and T2t) estimated using the IVIM model corrected for echo time dependence by fixing the transverse relaxation time T2p to 100 ms and 150 ms Median P25% P75% Min. Max. D [10−3mm2/s] 1.13 1.00 1.26 0.86 1.69 D∗ [10−3mm2/s] 15.4 11.7 16.0 8.1 23.4 ƒ 0.109 0.102 0.117 0.067 0.155 T2p = 100 ms ƒc 0.094 0.080 0.106 0.056 0.126 T2t [ms] 73.0 66.2 87.1 57.9 104.8 T2p = 50 ms ƒc 0.074 0.062 0.084 0.044 0.100 T2t [ms] 71.5 64.9 85.4 57.0 101.6 Wilcoxon signed rank test gave p < 0.01 between ƒ and ƒc for both cases. Pi%—the ith percentile.

To obtain an understanding of how sensitive the models are to intratumor heterogeneity, we performed a bootstrap analysis to obtain estimates of the variance of the signal fractions. This gave an empirical bootstrap distribution of ƒ, which allowed estimation of the intratumor variability. The results show a smaller variation in the corrected ƒc compared with the ƒ calculated with the standard IVIM model, with a standard deviation of 0.021 compared with 0.017, as shown in Table 4.

TABLE 4. Mean value, μ, and standard deviation, σ, from bootstrap analysis, μ (σ) ƒ, μ (σ) Model Interpatient Intrapatient Standard IVIM 0.1094 (0.0206) — (0.0055) Corrected IVIM 0.0821 (0.0170) — (0.0041) 4 DISCUSSION

This study has shown that conventional IVIM modelling at 3 T in breast cancer systematically overestimates the pseudo-diffusion signal fraction ƒ, demonstrating dependence of the relative signal contributions on acquisition echo time. Simulated breast lesion signal curves showed that single-voxel calculation of the corrected pseudo-diffusion fraction ƒc was hindered by high noise sensitivity of the IVIM model, with ƒ expected within the 95% CI of ƒc for currently available clinical DWI protocols at 3 T. Clinical patient data acquired in this study clearly illustrated TE dependence for IVIM in the breast, although reliable estimations of blood transverse relaxation time T2p, necessary for freely fitting an extended IVIM model, was not possible due to noise sensitivity and clinical acquisition time limitations. Accurate determination of the pseudo-diffusion fraction, ƒc, could thus not be performed. Using illustrative values of 100 and 150 ms for blood T2p, our parameter estimations resulted in 15% to 46% overestimation.

Our median parameter values for tissue diffusion time D = 1.1x10−3 mm2/s and diffusion time for the perfusion compartment D* = 15x10−3 mm2/s, and mean value for standard IVIM signal fraction ƒ of 11%, were consistent with those of previous studies.19, 26 The mean value for estimations of tissue T2t was around 72 ms, comparable to the transverse relaxation observed in healthy fibroglandular tissue. 18, 27

Blood T2 strongly decreases with decreasing oxygenation level. Stanisz et al reported a transverse relaxation time for in vitro blood samples of 275 ± 50 ms at 3 T, with 95% blood oxygenation,28 while other studies report T2 around 150 ms or shorter, decreasing with oxygenation level.20, 29 For 70% oxygenated venous blood, Chen et al reported blood T2 values of 60–100 ms, depending on mixing time.20 Our study shows a systematic reduction in effective T2 for b50 compared with b0, strongly suggesting that blood T2 is longer than the tissue T2; finding an average tissue T2 of 72 ms, our results are consistent with an apparent oxygenation level of at least 70% within breast lesions.

Where T2 relaxation times differ between components, as demonstrated in our different T2 values between b0 and b50 measurements, longer acquisition echo times will result in more signal loss from the short T2 component, which introduces error to the apparent fraction ƒ estimated using the standard IVIM model.14, 16 Differences between tissue and blood T2 values are smaller at 3 T than at 1.5 T (references17, 27); although this will reduce differences between true and apparent ƒ at 3 T compared with equivalent studies at 1.5 T, the effect may still adversely affect interpretation of IVIM analysis.

Our analyses show that the pseudo-diffusion fraction of IVIM is indeed dependent on echo time for breast lesions at 3 T, mirroring observations by others in different tissues,14, 16 and provide motivation for accounting for this effect. However, analysis of simulated data using T2p as a free parameter in the fitting procedure resulted in unstable parameter estimations, suggesting that a full fitting of the extended IVIM model may not be feasible within clinical acquisitions. Acquisitions for the extended IVIM model require additional scanning across echo time,16 which is costly in terms of time, money, and discomfort for patients. Our study explored the possibility of introducing a clinically feasible acquisition scheme that allows for studying the echo time dependence of IVIM parameters; however, our results show that estimation of blood T2p requires a greater sampling of echo times in the low b-value range.

IVIM fitting is challenging, with model parameter estimates being sensitive to noise correction strategy and fitting algorithm choice.5, 30-32 An improved understanding of the echo time dependence of IVIM measurements would allow for optimization of the echo time to either enhance or reduce the relative contributions of the two model compartments. Although an extensive acquisition scheme designed for parameter estimation using the fully corrected IVIM model is unlikely to be feasible in the clinic, it may yield useful insight for improving the clinical relevance of IVIM in breast cancer, including any additional value of the corrected model to breast cancer diagnostics. Results from this study indicate that conventional IVIM analysis, neglecting T2, gives overestimated pseudo-diffusion fractions ƒ in malignant breast lesions. The degree of overestimation is linked to the underlying true ƒ, meaning that, while the returned values must be considered inaccurate, the relative value may be still be useful, and an artificial “amplification” of the parameter through choice of echo time may actually increase discriminatory characteristics.

Closely related to the challenge of fitting IVIM is understanding parameter repeatability. It is known that the mono-exponential ADC is repeatable in lesions across many organs and protocols,33 and is considered a useful clinical measure.34 A review of IVIM studies showed that while D is repeatable, with coefficients of variation (CoVs) generally smaller than 20%, pseudo-diffusion parameters commonly gave CoV in excess of 20%, in some cases reaching above 50% (ƒ) and up to 100% (D∗).35 Capturing perfusion characteristics is challenging, and appropriate acquisition protocol design is still an open topic.36-38 Since IVIM parameters are dependent on acquisition parameters, including echo time and b value but also diffusion time,39 added value from extended acquisition protocols will rely on standardization. Additionally, both scanner field strength and local blood oxygenation affect the T2 relaxation times, modulating the strength of the echo time dependence, and should also be considered when critically assessing the potential benefit of extended acquisition and models.

The conventional IVIM model implicitly assumes no or slow water exchange between the two model compartments. The signal fraction ƒ of the pseudo-diffusion compartment is expected to reflect vascular properties, but it cannot readily be interpreted as a directly measurable biophysical property. It may still provide a clinically useful biomarker, however, if it exhibits sufficient sensitivity, specificity, and repeatability; accurate quantification (of a biophysical property) is not necessarily as important as precision, or generalization of the signal model to account for all relevant acquisition parameters. In this sense, understanding how T2 relaxation effects influence estimations of ƒ may provide useful information.

A limitation common to all IVIM studies is the high noise sensitivity of the model. Clinically, the segmented fitting approach as used in this study is often preferred, although work on alternatives, in particular Bayesian methods, has been suggested for this issue. Another limitation to this study is the small number of patients, and the challenge of collecting additional scans during a clinical examination of patients with malignant tumors. A larger cohort including both benign and malignant tumors would allow for testing of the clinical value of the corrected IVIM model, with regards to lesion characterization and discrimination. A natural continuation of this study would collect an extensive acquisition scheme with more measurements in both the b value and echo time dimensions over an increased echo time range, to provide more information on the echo time dependence, and potentially allow for echo time optimization of IVIM measurements in breast cancer.