MR imaging

Ten healthy subjects (age 21–34 years, seven males) were scanned using a 3 T Philips MR7700, software release 5.9.0, with a 32-channel head coil (Best, the Netherlands). The study was approved by the Swedish ethical review authority (Dnr 2020-00029) and all subjects signed informed consent. A subset of the data were included in a previous study evaluating methods for joint IVIM analysis of flow-compensated and non-flow-compensated dMRI data [10].

Three sets of dMRI scans based on a spin echo echo planar imaging sequence were acquired, designed for different types of IVIM analysis (upper panels in Fig. 1): (1) a three b-value scan with monopolar diffusion encoding gradients for a simplified IVIM analysis (sIVIM) [8, 11], (2) a ten b-value scan with monopolar diffusion encoding gradients for a conventional IVIM analysis assuming that the motion of blood is diffusion-like (pseudo-diffusion; referred to as the diffusive protocol) [8], and (3) two seven b-value scans with bipolar diffusion encoding gradients with one designed to be flow-compensated (FC) and the other one non-flow-compensated (NC) for joint analysis of FC and NC data assuming that the motion of blood is ballistic (referred to as the ballistic protocol) [12]. Timings of the diffusion encoding gradients are found in upper left and right panels of Fig. 1. Details about each scan are given below.

Fig. 1

Schematic illustration of the pulse sequences utilized in the study with actual timings displayed (upper left and right) along with the timeline of the full examination (lower panel). The scans with three and ten b-values used the monopolar gradient pulses while the scans with seven b-values used the bipolar gradient pulses to achieve a combination of flow-compensated (FC) and non-flow-compensated (NC) diffusion encodings. The width of boxes in the lower panel and their separation correspond to the actual scan duration and separation

All scans were acquired with a prototype pulse sequence enabling use of arbitrary gradient waveforms for diffusion encoding originally developed at Lund University [13] with the following common sequence parameters: TE = 80 ms, TR = 3700 ms, voxel size acquired 2 × 2 × 4 mm3 reconstructed to in-plane 1.875 × 1.875 mm2, 1 mm slice gap, SENSE = 1.9 (anterior–posterior direction), fat suppression by SPIR and gradient reversal. Dual readout was used for mitigation of Nyquist ghosting [14]. Six diffusion encoding directions (\(x,y,z, - x, - y, - z\)) were used for all diffusion-weighted scans (b > 0), while six repetitions were used for b = 0 to acquire the same number of scans regardless of b-value. The acquisition scheme was designed by looping through all combinations of b-values and diffusion encoding directions \(\left( ,d_ } \right)\), with \(d_ \, \in \,\left\ \right\}\), with simultaneous increments in both i and j to distribute acquisitions with the same b-value but different encoding direction over the scan time (see examples in lower panel of Fig. 1). All scans were acquired in a single session and each scan was repeated without moving the subject.

MR imaging: sIVIM protocol

The sIVIM protocol utilized monopolar diffusion encoding gradients using the three b-values 0, 200, and 800 s/mm2 with the number of repetitions of each b-value obtained from a Cramer-Rao Lower Bound-based optimization resulting in 2, 3 and 1 repetitions of the respective b-values for the given scan time set to be 4.5 min [15]. These repetitions were implemented as the b-value sequence 0, 200, 800, 200, 0, 200 s/mm2, repeated six times as described above to acquire data for all diffusion encoding directions.

MR imaging: diffusive protocol

The diffusive protocol utilized monopolar diffusion encoding gradients and a sequence of ten b-values similar to what is commonly seen in the literature: b = 0, 5, 10, 20, 30, 50, 100, 200, 500, and 800 s/mm2. The b-values were acquired in a “low–high” order i.e. b = 0, 800, 5, 500 etc. Total scan time was 7.5 min.

MR imaging: ballistic protocol

The ballistic protocol utilized two scans with bipolar diffusion encoding gradients (FC and NC) with seven b-values: 0, 5, 10, 20, 30, 100, and 200 s/mm2. The b-values were acquired in a “low–high” order i.e. b = 0, 200, 5, 100 etc. Total scan time was 10.5 min.

Preprocessing of data

To minimize spatial misalignment due to motion and eddy currents, and to get all images to a common space, all images were registered to the first b = 0 image of the examination (first b = 0 image of the first acquisition of the sIVIM protocol was originally called mono3b, was the naming was changed during the review process. The instance was missed in the revision.) using affine registrations (eddy_correct script in FSL v. 6.0.4) [16]. Non-brain tissue was masked out by brain extraction (BET in FSL v.6.0.4) [17]. Three regions of interest (ROIs) were manually delineated in prefrontal white matter, centrum semiovale, and cerebellum for subsequent ROI-based analyses (examples in Fig. 2).

Fig. 2

b = 0 signal drift. The left column illustrates the investigated regions of interest (ROIs) for an example subject. The middle column shows how the average b = 0 signal, normalized to 100 for the first data point, fluctuated during the acquisition of an example scan from each protocol for the same example subject as in the left column. The colors correspond to different correction methods (see legend in right column) and each section corresponds to a given protocol. The right column shows group summaries where displayed data is the ROI median signal difference between each scan’s first and last b = 0 acquisition aggregated over all scans for a given protocol. Boxplots show min, 25th percentile, median, 75th percentile, and max.

Signal drift correction

Three methods for signal drift correction were utilized: (1) a global temporal correction as first suggested by Vos et al. [2], (2) a voxelwise temporal correction corresponding to this global method [6], and (3) a spatiotemporal correction as suggested by Hansen et al. [6].

The temporal corrections fitted a second order polynomial to the b = 0 signal as a function of time where the global method fitted a single polynomial to the median signal in the brain, while the voxelwise method fitted one polynomial for each voxel:

$$S\left( = 0} \right) = k_ + k_ n + k_ n^ ,$$

(1)

where \(S(n|b_ = 0)\) is the signal at acquisition \(n\) given that it is a b = 0 acquisition, and \(k_\) are the polynomial coefficients to be estimated. The polynomial coefficients were estimated by finding the least squares solution to data. The correction was then applied to all data, regardless of b-value, with corrected signal values calculated as:

$$S_}}} \left( n \right) = \frac_ + \hat_ n + \hat_ n^ }},$$

(2)

where \(S_}}} \left( n \right)\) is the corrected signal at acquisition \(n\) and \(\hat_\) are the polynomial coefficients estimated from the b = 0 signal, i.e. the subset \(}nb_ = 0\}\).

The spatiotemporal correction followed a similar procedure, but did instead fit a second order polynomial in both time and space:

$$S\left( = 0} \right) = p_ \left( \right) + p_ \left( \right) \cdot n + p_ \left( \right) \cdot n^ ,$$

(3)

where \(p_ \left( \right)\) is a polynomial containing all zeroth, first and second order combinations of the spatial coordinates x, y, and z, resulting in a total of 81 polynomial coefficients to be estimated. \(p_ \left( \right)\) can be found in Appendix A [6]. To avoid overfitting to outliers or data points with high leverage, the polynomial coefficients were estimated by bisquare regression as suggested by Hansen et al. [6, 18]. Before fitting this polynomial to data, a voxelwise normalization to the first b = 0 image was applied to remove the underlying anatomical signal distribution in the b = 0 images. As for the temporal corrections, the corrected data was calculated as:

$$S_ \left( \right) = \frac \right)}}_ \left( \right) + \hat_ \left( \right) \cdot n + \hat_ \left( \right) \cdot n^ }},$$

(4)

where \(\hat_\) are the polynomials estimated from the b = 0 signal, i.e. the subset \(}nb_ = 0\}\).

Simulation of acquisition order

To study the effect the order in which data are acquired has on derived model parameters, data with acquisition protocol ordered by b-value (ordered protocol) were generated, from which parameter maps were derived and compared with those obtained from the actually acquired data (mixed protocol). Data corresponding to the ordered protocol were generated by reordering data corrected with the spatiotemporal polynomial correction based on b-value, such that an increasing b-value order was obtained with all encoding directions acquired in direct succession. This is in contrast to the mixed protocol used for acquisition of data where b-values are mixed and the different encoding directions for a given b-value were distributed over the scan time. Following the reordering, the inverse spatiotemporal polynomial correction was applied to simulate the signal drift that previously had been corrected for. This enabled comparison of uncorrected data from the original mixed protocol and uncorrected data from a simulated ordered protocol as well as their corresponding parameter maps.

IVIM parameter estimation

Before parameter estimation, all scans were directionally averaged with a geometric average to mitigate potential effects of background gradients [19, 20]. To remove potential scaling differences between the FC and NC scans from the ballistic protocol, the two scans were each normalized by dividing with the median b = 0 signals within the brain mask before being combined into a single data set. Specifically tailored analysis approaches were then used to obtain a set of IVIM parameters from each protocol, as described below.

IVIM parameter estimation: sIVIM

Data from the sIVIM protocol were used to estimate the IVIM parameters D and f via a segmented algorithm [8, 21]. First, a monoexponential signal representation was fitted to data corresponding to the two non-zero b-values:

$$S\left( b \right) = Ae^ .$$

(5)

Since only two b-values were used, analytical solutions could be used to first calculate the diffusion coefficient D as:

$$D = \frac - \ln S_ }} - b_ }},$$

(6)

where Si is the signal at b-value bi. The intercept term A was then calculated as:

In a second step, the perfusion fraction f was calculated as:

$$f = 1 - \frac }},$$

(8)

where S0 is the signal at b = 0.

IVIM parameter estimation: diffusive regime

Data from the diffusive protocol were used to estimate the IVIM parameters D, f, and D* of the diffusive regime by a least squares fit to data with the following biexponential signal representation:

$$S\left( b \right) = S_ \left( \right)e^ + fe^ }} } \right),$$

(9)

where D* is the pseudo-diffusion coefficient [8].

IVIM parameter estimation: ballistic regime

Data from the ballistic protocol were used to estimate the IVIM parameters D, f, and vd of the ballistic regime by a least-squares fit to data with the following biexponential signal representation:

$$S\left( \right) = S_ \left( \right)e^ + fe^ }} e^ v_^ }} } \right),$$

(10)

where c is the flow encoding factor, and Db and vd are the diffusion coefficient and velocity dispersion of blood, respectively [8, 12]. To stabilize the fit, Db was fixed as suggested by Wetscherek at el. [22] with a value of 1.75 µm2/ms in accordance with Ahlgren et al. [12]. It is worth noting that Eq. 10 could have been used for processing data from the diffusive protocol as well, but Eq. 9 was used for comparability to the literature.

Software

All processing of data, unless otherwise stated, was implemented in Python 3.11 with direct calls to packages numpy 1.24, matplotlib 3.7.1, scipy 1.10.1, statsmodels 0.13.5, and nibabel 5.0.1. Project specific functions, including b-value optimization, signal drift corrections, and IVIM parameter estimation methods, were gathered in a publicly available python package (https://github.com/oscarjalnefjord/ivim, SHA-1 hash 4f063a1).

Statistical analysis

For analysis purposes, signal drift was approximated as signal difference between each scan’s first and last b = 0 acquisition normalized to the scan duration. Difference in signal drift among brain regions, and difference in IVIM parameter values (both average and difference of repeated scans) between correction methods or b-value order were statistically tested using the Friedman test with the exact method for calculating the test statistic (Python package permutation-stats [23]) [24]. If a Friedman test returned a p-value < 0.05, pairwise signed-rank tests were employed where p < 0.05 was considered statistically significant.

Identification of example data

A data set suitable for displaying typical results was identified based on the signal drift in all scans. Data from all subjects, repetitions, regions of interest (ROIs) and scans were ranked based on their distance from the median signal drift of the particular ROI and scan. The subject and repetition with the overall lowest rank was defined as a typical example and used throughout the results section for displaying example data on a subject level.