The importance of correcting for signal drift in diffusion MRI

To investigate previously unreported effects of signal drift as a result of temporal scanner instability on diffusion MRI data analysis and to propose a method to correct this signal drift.


INTRODUCTION
Diffusion MRI is a powerful imaging technique that can reveal information about the neural fiber architecture in vivo (1,2). By sensitizing the magnetic resonance imaging (MRI) signal to self-diffusion of water molecules in the tissue, the measured signals can be interpreted as giving quantitative microstructural as well as orientational information. This has resulted in an explosive use in clinical research studies into diseases including Alzheimer's disease (3) and epilepsy (4), and a translation into clinical practice, for example, for improving neurosurgical outcome (5).
Results and conclusions in all these studies need highquality data to support the claims. In diffusion MRI, data quality is particularly relevant because of the large number of issues and artefacts arising from the data acquisition. Whereas subject motion is an artefact caused merely by the movement of the subject over the scan duraction, other artefacts are caused by the inherent need for acquisition speed in diffusion MRI. The single-shot echo-planar imaging (EPI) readout results in geometric distortions due to susceptibility differences as well as chemical shift artefacts, and the strong diffusion gradients cause eddy current induced artefacts. The combination of diffusion gradients and echo-planar imaging (EPI) readout put a high load on the system, causing heating effects that decrease the scanner's temporal stability (6,7). This instability has already been shown to decrease fat suppression efficiency (7).
In this work, we describe a previously unexplored detrimental effect of temporal instability of the scanner system on diffusion MRI data: a decrease in global signal intensity. We show evidence of signal drift in phantom data and in vivo data, and that this effect is present on scanners from multiple vendors (GE, Philips, and Siemens). We demonstrate the effects on commonly used DTI and DKI parameters as well as fiber tractography, and present an easy to use postprocessing method to correct for this artifactual signal decrease. The proposed signal drift correction method is implemented in MATLAB (MathWorks, Natick, MA, USA) and is made freely available on the author's website (http:// cmictig.cs.ucl.ac.uk/people/research-staff/50-sjoerd-vos) and in the MATLAB File Exchange website (http://www. mathworks.com/matlabcentral/fileexchange/54861) and will be incorporated into the new release of ExploreDTI 1 (www.exploredti.com) to facilitate rapid and widespread correction of signal drift in diffusion MRI data.

Data
In this work, we use four distinct experiments to characterize signal drift, to demonstrate its presence on all three investigated scanners from different vendors, and to quantify its effect on in vivo diffusion MRI analyses.

Signal Drift Characterization
First, the signal drift characteristics were investigated using acquisitions of a gel phantom and simulations of diffusion-weighted signals to determine signal drift patterns over time, reproducibility, possible causes, and the influence of acquisition protocols.

Phantom Acquisition
To determine the extent of signal drift and its characteristics in diffusion MRI acquisitions, a gel phantom with isotropic diffusion properties was scanned with a very high load on the gradient system. This acquisition consisted of 11 b ¼ 0-images and 25 gradient orientations on the half-sphere with b-values of 1000, 2000, 3000, and 9000 s/mm 2 for an acquisition time of 15 min. The experiments were performed on a routine clinical 3 T Philips Achieva system, without fat suppression. The ordering of the DWIs was done in two ways: (1) ordered by b-value from low to high; and (2) in randomized order. In both cases, the b ¼ 0-images were interspersed at regular intervals throughout the acquisition. Scans were acquired with and without updating the center frequency (f 0 ) after each acquired volume ("Dynamic Stabilization"), where the default is to only define this once at the beginning of the scan (so without updating). All scans were performed twice to check for consistency, and f 0 frequency was logged after each acquired volume.

Simulations
To further explore signal drift characteristics and the effect of drift on DTI and DKI metrics, simulations were performed to be similar to the acquired phantom data. The same b-values and gradient orientations were simulated with both ordered and randomized acquisition setups. The diffusion properties were based on those estimated from the phantom data, defining isotropic diffusion with an MD equal to that estimated in the phantom taking in account the effects of noise on the MD estimation. Therefore, the MD used as input for the simulations (MD ¼ 0.055 Â 10 À3 mm 2 /s) is lower than the actual MD measured in the phantoms (0.059 Â 10 À3 mm 2 /s). The diffusion tensor model, without kurtosis, was used to simulate the data with signal-to-noise ratio (SNR) equal to the estimated SNR of the phantom data, which was 44 6 13, by adding Rician noise to generate a dataset without drift (called the "unaffected" dataset). To compare this with a simulated dataset with drift, signal drift was applied based on the estimates of drift magnitude from the phantom data (called the "drift" dataset). The drift was defined using the quad-ratic function 100 À 0.0183 Â n À 2.25 Â 10 À4 Â n 2 , which results in %5% signal drift over the range of n ¼ 1-110 volumes. Adding signal drift was done before adding the Rician noise to the data, such that the actual SNR decreases for increased signal drift, as is the case for the acquired data. Similar simulations were performed for anisotropic systems using cylindrically symmetric diffusion tensors, without kurtosis, with simulated FA values of 0.41 and 0.81 and MD equaling 0.81 Â 10 À3 mm 2 /s. Values for SNR and drift magnitude in these simulations were equal to those used in the isotropic phantom simulations. All simulations consisted of 32,000 voxels (20 Â 40 Â 40) which was identical in size to the ROIs used to select data form the phantom measurements.

Multi-Vendor Comparison
To determine the magnitude of signal drift in realistic neuroscientific acquisition protocols across multiple scanners [e.g., (8,9)], a 15-min multi-shell acquisition protocol with b-values of 600, 1600, and 2500 s/mm 2 was created. Data were acquired of phantoms with isotropic diffusion properties on the same 3 T Philips Achieva, a 3 T GE MR750 scanner, and a 3 T Siemens Trio scanner to investigate signal drift behavior across scanners from different vendors. Across scanners there is a slight difference in TE and TRas a result of different maximum gradient amplitude and amplifier power-leading to a different number of images in the 15-min scan ranging from 140 to 150 DWIs, with b ¼ 0-images interspersed every 10 DWIs. Acquisitions were done with and without the respective dynamic f 0 updating options ("Dynamic Stabilization" on Philips, "Real Time Field Adjustment" on GE, and "Frequency Stabilization" on Siemens-part of a 2011 Works-in-Progress package for diffusion EPI).

In Vivo Acquisitions
Three individual datasets from the MASSIVE dataset [www.massive-data.org and (10)] were taken and examined for the impact of signal drift on diffusion MRI analyses. These datasets were selected based on having the lowest (dataset 1), median (dataset 2), and maximum (dataset 3) signal drift from the multi-shell datasets in MAS-SIVE, to illustrate the range of effects signal drift can have on diffusion MRI analyses. A detailed description of the data is given in (10), but these datasets consist of: 120 images acquired with b-values of 0, 500, 1000, 2000, and 3000 s/mm 2 , and a varying number of gradient directions for each b-value. Specifically, the three datasets included in this work consisted of: 11/16/31/31/31, 10/16/27/34/33, and 11/13/26/41/29 images (at b ¼ 0, 500, 1000, 2000, 3000 s/mm 2 , respectively). Importantly, the order in which these images were acquired was randomized. Acquisition details included: 3 T Philips Achieva, matrix 96 Â 96, FOV 240 Â 240 mm 2 , 2.5 mm slices for an isotropic voxel size of 2.5 mm 3 (11), SENSE acceleration ¼ 2.5, TE/TR: 100/7500 ms. Spectral Presaturation with Inversion Recovery fat suppression was used. Scan time per session was 15 min. All datasets were acquired on different days.

Signal Drift Correction
Mean magnitude intensity was obtained for each b ¼ 0image within a region of interest (either the entire phantom or the brain). Global signal decrease was then estimated with a linear and quadratic fit through all these mean signal intensities: and respectively, where n is the image index of the scanned images, Sðn j b ¼ 0Þ is the mean signal in image n given that it is a b ¼ 0-image; d, d 1 , and d 2 are signal drift coefficients; and s 0 is the signal offset at n ¼ 0. From the obtained fit of s 0 and d or s 0 , d 1 , and d 2 , a rescaling factor can then be calculated and applied for each image in the series to obtain a drift-corrected image: and S q ðnÞ ¼ SðnÞ 100 for the linear and quadratic correction, respectively, where S(n) is the uncorrected, or raw, image intensity of image n, and S l (n) and S q (n) are the corrected signal intensities in that image normalized to 100 for linear and quadratic corrections, respectively. This normalization facilitates comparison between different datasets and would allow for merging these datasets.

Data Analysis
Motion and distortion correction was performed on each acquired dataset with ExploreDTI (12) using a 12parameter affine registration in elastix (13). Each dataset was processed further with and without signal drift correction. Diffusion kurtosis tensors were estimated as in (14), wherein the diffusion tensor is also estimated.

Signal Drift Investigation
The quality of the linear and quadratic fits through Sðn j b ¼ 0Þ is compared to determine whether a quadratic term is needed to describe signal drift over time. The effects on signal stability over time with and without f 0 updating and with and without signal drift correction are investigated per b-value for the phantom datasets.

Frequency Drift Investigation
The phantom experiments in which f 0 was logged are used to determine the stability of the scanner in terms of center frequency drift, a known factor in EPI-based acquisitions [e.g., (6)].

Effect on Tensor Estimation
The performance of the kurtosis tensor fits were examined in the "Signal drift characterization" phantom data by looking at the residuals per acquired imaging volume after the voxel-wise fit, with the residuals defined as the absolute difference between the acquired signal and the signal reconstructed from the model fit. The median and interquartile ranges are shown for each dataset and DWI, which are the values and ranges over the entire phantom scan.

Effect on Quantitative Diffusion Measures
To investigate the effect of signal drift on the analysis of diffusion MRI data, comparisons between the dataset before and after signal drift correction were performed both quantitatively and qualitatively. The following quantitative measures were investigated for the phantom datasets and simulated data: all three eigenvalues (k 1 , k 2 , and k 3 ), fractional anisotropy (FA), mean diffusivity (MD), mean kurtosis (MK). For the in vivo datasets the FA, MD, and MK were also analyzed, as were the first eigenvectors (FEs).

Effect on Fiber Tractography
Qualitatively, the results of fiber tractography of the corticospinal tracts (CST) were examined for the in vivo data. Tractography was performed using the FE of the diffusion tensor as estimated during kurtosis tensor estimation (15). Tracts were seeded regionally from a single voxel within the bundle on an axial slice at the level of the decussation of superior cerebellar peduncles (16).

Signal Drift Characterization
The effects of signal drift as a function of scan progression are shown in Figure 1. All four datasets without frequency stabilization had pronounced signal loss of around 5-6% over the 15-min scan sessions. Signal loss occurred as a quadratic effect, as demonstrated by the higher accuracy of the quadratic fit over the linear fit, with increasing signal loss with scan duration (Fig. 1). Consequently, the signal drift correction based on the quadratic fit resulted in normalized signal intensities of the b ¼ 0-images that were more stable over time than those after linear correction. The signal drift in datasets with dynamic updating of f 0 showed more varying patterns, with moderate signal loss in two cases (0.5 and 1.5%) and slight signal increase in two other cases (0.2 and 0.5%), and was not correlated with acquisition order. For these datasets, the signal differences were smaller, and the quality of the linear and quadratic corrections appeared to be roughly equal. When plotting the signal intensities against the b-value of each acquired volume, ordered first by b-value and then by volume number, the effect of signal drift becomes extremely clear (Fig. 2). The 6% signal decrease as a result of drift causes some of the volumes with lower b-values to have lower signal intensity than volumes with more diffusion-weighting. The signal drift corrections used in this experiment were based on quadratic fits, and restored the signal intensities per bvalue to very stable levels. It is also clear that the randomized acquisition order shows larger variation in how strongly volumes within each b-value are affected compared to the ordered acquisition (Fig. 2a vs. 2b, respectively). In the ordered acquisition (b), the images at b ¼ 1000 s/mm 2 are only marginally affected by signal drift whereas the images scanned last, at b ¼ 9000 s/ mm 2 , have a considerably lower signal intensity before correction for all images. This analysis was also FIG. 1. Normalized mean signal intensities in the b ¼ 0images for the phantom datasets as a function of scan progression (i.e., how long after the scan started they were acquired). The left column shows the results from the randomized acquisition, the right column the ordered acquisition. The top two rows are the two repeats without frequency stabilization, the bottom two rows are the two repeats with stabilization. The black dots indicate the original signal intensities, the red line indicates the linear fit through these points, and the blue line indicates the quadratic fit through these points. The red and blue dots indicate the corrected mean intensities in the b ¼ 0-images from the linear and quadratic fits, respectively. The black dashed line shows the level of the normalized signal (100). The values in each panel indicate the mean 6 standard deviation over the b ¼ 0-intensities after linear (S0_l) and quadratic correction (S0_q).
performed on the simulated data, and the comparison of signal drift and drift correction between the acquired phantom data and the simulated datasets is shown in the Supporting Information Figure S1.
Even though SNR in the images decreases with bvalue, applying the correction based on signal drift estimates from the high-SNR b ¼ 0-images is a valid correction step for almost all SNR levels and b-values-as detailed in the Supporting Information Figures S2-S5.

Frequency Drift
The results of the phantom experiments where the center frequency was logged show approximately linear f 0 drift of around 50 Hz over the 15-min session (Fig. 3). There is no clear difference between f 0 drift in the experiments where f 0 was and was not dynamically updated (bottom vs. top panels in Fig. 3, respectively), nor between the randomized and ordered acquisitions (left vs. right columns in Fig. 3, respectively).

Effect on Tensor Estimation
Signal drift, being a modification of signal intensity as a function of image index, will also appear as a covariate in the residuals of the model fit. Figure 4 shows this effect on the absolute residuals of the diffusion kurtosis model fitting for each volume from the phantom data. After signal drift correction the residuals in each shell are more homogeneous, indicating that the fit is more equally influenced by all DWIs.

Effect on Quantitative Diffusion Metrics
Estimated DTI metrics from the phantom data are shown in Figure 5, which shows that in the ordered acquisition the signal drift caused erroneous estimation of most of these DTI metrics whereas a randomized acquisition resulted in smaller errors in the estimated DTI metrics. The MK measure, being the metric based on the higherorder terms of the signal equation, gives the most specific indication of the risk of signal drift in ordered acquisitions. In an ordered acquisition signal drift affects the higher b-value data more severely, which leads to artifactually high signal loss at high b-values and an underestimation of the MK (as shown in Fig. 5). In DKI analyses, where positive kurtosis values typically indicate restricted diffusion, this would lead to an underestimation of restriction.
The simulated data show similar results (Fig. 5, right panel). Here, in contrast with the phantom data, the true DTI metrics are known. Although a fully isotropic system was simulated (FA ¼ 0), due to noise and eigenvalue sorting the FA is estimated to be around 0.15 in the simulation where no signal drift was present. Signal drift affected the FA, causing an underestimation in ordered and an overestimation in randomized acquisition. After signal drift correction both sets returned to the value estimated in the driftfree reference data. The MD was affected more strongly by drift, being overestimated by roughly 12% in the ordered acquisition and 5% in the randomized acquisition. It is important to note that the FA and MD errors were caused by overestimations of all eigenvalues in the ordered data and overestimations of the first and second eigenvalues in randomized data. The kurtosis for the ordered acquisition shows a similar underestimation as in the phantom data. Here, however, signal drift correction leads to a near perfect correction of MK. The MK from the randomized acquisition is not affected by signal drift, as could be expected since drift affected all shells equally (a multiplication in signal magnitude is simply an addition in the log of the signal). For the anisotropic simulations, the results are very similar: scalar diffusion metrics are prone to misestimation in the presence of signal drift (Supporting Information Fig. S6). The misestimation is most pronounced for ordered acquisitions, resulting in more than 5% errors in MD and FA, and considerable erroneous estimates of MK. For randomized FIG. 3. Plot of center frequency (f 0 ) drift versus time since scan start for the phantom datasets from the "Signal drift characterization" experiment. This is the change in f 0 at the beginning of each volume with respect to f 0 at the beginning of the scan. The left column shows the results from the randomized acquisition, the right column the ordered acquisition. The top row shows the experiments without f 0 updating with identical intensity ranges for all four graphs; the bottom row with f 0 updating. Each graph shows the two repeats of this experiment (one indicated by red dots and the other by black stars). For the acquisitions with f 0 updating the measured frequency drift was corrected before acquiring the next imaging volume.
acquisitions the effects are most subtle, mostly affecting the DTI metrics.
Comparing the phantom and simulated data, the FA values in the phantom data are slightly higher, 0.18 versus 0.15. In both cases, the values are higher than the true FA of 0, but it is also clear that signal drift-correction gives similar values as the data that was not affected by drift (in case of simulated data). In general, the phantom data exhibited a small effect of drift correction on diffusion tensor metrics compared to the simulations, whereas the kurtosis metric was affected more in the phantom data than the simulated data.  4. Absolute residuals (median and interquartile range) for uncorrected (left column) and corrected (right columns) phantom datasets for the randomized acquisition without (a) and with (b) f 0 updating and the ordered acquisition without (c) and with (d) f 0 updating. One can clearly appreciate that the residuals are more homogeneous for each b-value after signal drift correction. This is most pronounced for the dataset without f 0 updating (a and c) but homogeneity in residuals is also improved for the dataset with f 0 updating (b and d).
The residuals are smaller after signal drift correction in the randomized than the ordered datasets.

Multi-Vendor Comparison
In the 15-min neuroscientific acquisition protocol signal drift was observed on all three investigated scanners, with signal drift magnitude varying up to 2% across scanners (Fig. 6). On the Philips scanner (Fig. 6, top row), signal drift of around 1% was ameliorated by frequency stabilization which reduced signal drift to around 0.5%. On the GE scanner, signal drift was around 0.5-1% (Fig. 6, middle row). A second repeat of these figures resulted in drift of around 0.6% irrespective of field updating, indicating that dynamic stabilization on this scanner does not influence signal drift magnitude. On the Siemens scanner signal drift magnitude seems the largest, around 1.5-2%, when the scanner had not been heavily used before the diffusion MRI acquisition (Fig. 6, bottom row). After a period with a heavy gradient load, signal drift was around 0.5-1%. These patterns were irrespective of whether frequency stabilization was turned on or off. On the Philips and GE scanners, no differences were observed between scans with or without prior heavy gradient use. These experiments were all performed twice for reproducibility, and showed strongly reproducible results. The multi-vendor comparison reiterates that a quadratic drift correction provides better results than linear correction.
Signal drift magnitude was not the only property of drift that varied between scanners. As can be observed from Figure 6, drift on the GE scanner seems to vary linearly with elapsed time during a scan, whereas a quad-ratic term is present in drift in the data from the Philips and Siemens scanners. Specifically, signal drift in the data from the Philips scanner is concave downwardsdrift increasing more rapidly with scan progressionwhereas on the Siemens scanner it is concave upwardsdrift reducing, or stabilizing, as the scan progresses.

In Vivo Data Analysis
Significant global signal decreases were also observed in all investigated in vivo datasets. Similar to the phantom datasets, the quadratic fit through the b ¼ 0-intensities provided the best results. Based on these quadratic fits the signal decrease from first to last image in each of the 15-min sessions was 5.0, 8.3, and 16.7%, corresponding to a signal loss of 0.04, 0.07, and 0.14% per image (Fig. 7, Table 1).

Effect on Quantitative Diffusion Metrics
Differences in diffusion tensor and kurtosis tensor metrics between corrected and uncorrected dataset are shown in Figure 8 for all three datasets. Pronounced differences are present for all three datasets and throughout the white matter. The pattern of over-or underestimation in the uncorrected datasets is very variable between the datasets. For all datasets, the histograms of the FA difference are skewed towards positive difference, that is, where the uncorrected datasets underestimate the "correct" value. Such patterns are not visible for the MD FIG. 5. Calculated DTI/DKI metrics for the acquired real phantom data (left) and simulated data (right). Different metrics are shown in each panel: the three eigenvalues of the diffusion tensor (k 1 , k 2 , k 3 ); the mean diffusivity (MD); fractional anisotropy (FA); and mean kurtosis (MK). Each panel shows data from the ordered and randomized acquisition. For the phantom data, the frequency-stabilized, driftaffected, and drift-corrected data are shown in blue, red, and green, respectively. For the simulated data, the unaffected, drift-affected, and drift-corrected data are shown in blue, red, and green, respectively. The colored boxplots indicate the 25th-75th percentiles, the white line within the box the median, and the error bars the 5th-95th percentile range. For the simulated data, the solid gray line plotted behind the boxplots indicates the value for the unaffected data and the dashed line indicates the value set in the simulations. and MK, which are more centered around zero mean deviation but with considerable variance.
The orientation of the diffusion tensor's FE is also affected by signal drift, giving a considerable angular deviation between the FEs from the corrected and uncorrected datasets (DFE, in Fig. 8 and Table 1). For voxels with an FA larger than 0.4, the modes of the angular deviation histograms are 1.0, 1.3, and 3.0 . Large portions of the angular deviations are above these modes-66, 67, and 59% of the histogram areas, respectively-and the angular deviation can greatly exceed a few degrees. Table  1 also states these percentages for all voxels with FA larger than 0.2 and 0.7.

Effect on Fiber Tractography
The effects of erroneous FEs is shown most pronounced by the effects on tractography. An example of the CST is shown in Figure 9. The tracts end up in different gyri when signal drift correction is not performed, caused by The black dots indicate the original signal intensities, the red line indicates the linear fit through these points, and the blue line indicates the quadratic fit through these points. The red and blue dots indicate the corrected mean intensities in the b ¼ 0-images from the linear and quadratic fits, respectively. The black dashed line shows the level of the normalized signal (100). The values in each panel indicate the mean 6 standard deviation over the b ¼ 0-intensities after linear (S0_l) and quadratic correction (S0_q). Phantom data from the Philips and GE scanners are scanned without (left) and with (right) dynamic field updating. Dynamic frequency updating had no effect on the Siemens scanner (data not shown), but there was a strong effect of signal drift reduction with a "cold" scanner (left) or a scanner that had been through 15 min of scanning (right). (17)].

DISCUSSION
The aim in this work was to demonstrate the presence of signal drift in diffusion MRI data, and to present a simple correction that could be applied as part of the image processing pipeline after acquisition. Diffusion-weighting causes signal attenuation, and any uncorrected signal drift over the scan duration will thus be interpreted as overestimated diffusion. Even for the 15-min acquisitions included in this study, pronounced signal drift was observed. Although our results show effects on diffusion tensor and kurtosis tensor metrics, these effects are generalizable to any diffusion MRI method.
In functional MRI, which is typically also acquired using EPI, signal drift is a well-known effect (18,19). In the time-course analysis of functional MRI data, this signal drift is commonly included as a confounding factor. Although an identical approach would not work for diffusion MRI data, we have proposed a similar, quadratic, correction to compensate for signal drift. Figures 1 and 6 show that signal drift is an issue common throughout scanners from all three investigated vendors, and that the proposed quadratic correction restores temporal stability in the acquired signals. In vivo, the quadratic fit also provides better results than the linear fit, restoring the normalized signal intensities of all b ¼ 0-images more closely to equal intensities.

Effect on Quantitative Diffusion Metrics
Signal drift correction, overall, resulted in more homogeneous residuals (Fig. 4), indicating that all information was integrated into the fit more equally. Additionally, Figure 5b and Supporting Information Figure S6 show that signal drift correction causes the range of all investigated DTI and DKI metrics to correspond better to the reference data in both isotropic and anisotropic simulations. Here, the drift-corrected data had the same median values and range as the data that was simulated without drift, where the data corrupted by signal drift deviated significantly.
The estimated FA in the phantom data was higher than the estimated FA from the simulations, 0.18 versus 0.15, which could be caused by imperfections resulting from acquisition and processing. The signal decay remained Gaussian, as was expected in a medium without restrictions, as evidenced by the median MK being zero.
Although results of the simulations are very similar to the results obtained from phantoms they are not a perfect match. In the simulations, we assume that the drift is perfectly quadratic and that the only effect of drift is signal loss. In reality, the drift only approximates a quadratic function.
In vivo, the specific effects of signal drift on DTI and DKI metrics are strongly dependent on when in the session each diffusion-weighting gradient orientation is acquired. If, either by design or by chance, all gradient directions at the end of a scan are primarily left-right diffusion-weighted, signal drift will cause a systematic bias along that left-right axis. With the randomly ordered gradients in our design, we did not see such a bias, but   rather observed nontrivial alterations in the estimated diffusion tensor and kurtosis tensor and the resulting FA, MD, MK, and FE (Fig. 8). Any erroneous estimation of these quantitative measures will affect results obtained from these datasets. The most likely effect is an increase in variance within or between groups, leading to a reduction in statistical power to detect genuine effects.

Effect of Fiber Tractography
The effect of angular errors of the FE can lead to significant differences in fiber tractography results. Even for relatively low angular errors per voxel, on the order of 1 , the effects on tractography can be pronounced due to error accumulation along the tract (17). This is most clearly shown in Figure 9, with tracts ending up in different gyri with and without signal drift correction. For group studies of connectivity, any signal drift would lead to a reduction in statistical power, similar to the effects of quantitative diffusion metrics. For singlesubject cases, as for instance in neurosurgical planning [e.g., (5,20-22)], any angular errors caused by signal drift could significantly alter any planning approaches or clinical decisions based on tractography.

Causes of Signal Drift
Although the main purpose of this work was to characterize signal drift and to provide a simple correction strategy, we could speculate about the origin of the drift. Scanner instability as a result of gradient heating can cause signal drift in various ways. First, heating can cause B 0 field drift and thus f 0 drift as previously described in (7) and also shown in Figure 3. It should be noted that although the acquisition protocols in the multi-vendor comparison were near-identical, the load on the gradient systems between scanners is not necessarily equal. Differences in maximum gradient amplitude and amplifier power could cause a different amount of heating, with heating also depending on the internal construction of the scanner and cooling systems. Depending on the chosen excitation and fat suppression techniques, f 0 drift could lead to reduced efficiency of the excitation pulse and/or a shift of the fat suppression pulse toward the water peak, both causing a decrease in excited magnetization. Another possible cause is the effect of magnetic field drift on Nyquist ghost correction. Typically this is calibrated once at the beginning of the scan and f 0 drift may cause the correction to be less effective, increasing the amount of signal attributed to the ghost. Finally, transmit coils heat up as a result of the high duty-cycle and this heating potentially results in an altered transmission energy and thus a flip angle other than the prescribed 90 . Signal drift was observed on all three included scanners, but magnitude and temporal pattern of the drift varied considerably. Dynamic frequency updating reduced signal drift magnitude on the Philips scanner, and resulted in more stable b ¼ 0-intensities over time than without dynamic updating (Figs. 1 and 6). On the GE and Siemens scanners, dynamic field or frequency stabilization had no effect. Signal drift patterns varied as well, with Philips data having strong quadratic downwards drift (i.e., more signal loss with each acquired image); GE exhibiting a linear drift; and Siemens a quadratic upwards drift (i.e., less signal loss with each acquired image). These patterns give a strong indication that signal drift is caused by multiple different phenomena on different scanners. On the Philips scanner the quadratic drift term was absent when frequency stabilization was on, indicating frequency drift might be one of the causes of signal drift. This can be supported theoretically: when f 0 is initially set correctly to yield maximum signal, the constant term and quadratic term with a negative coefficient are the leading term in a power series expansion of the signal as a function of frequency offset.
With the observed linear f 0 drift (Fig. 3), the signal is expected to decrease quadratically with time since scan start. In the GE data, however, stabilization did not reduce signal drift, indicating no frequency drift was present. Coupled with the absence of a quadratic terms, this further strengthens the rationale that frequency drift causes a quadratic signal drift. The Siemens data from a "cold" scanner shows a quadratic upward signal drift pattern, with signal drift reducing with each acquired volume, could indicate stabilization or a warming-up pattern that matches with the results from showing that data from a warmed-up scanner had less signal drift. This, however, is not linked by frequency drift given that frequency stabilization did not alter results.

Ordered versus Randomized Acquisition
For both the acquired phantom data and the simulations, the ordered acquisition gave the largest error in DTI and DKI metrics (Fig. 5). In the simulations, signal drift correction performed a near-perfect correction, with driftcorrected data yielding nearly identical DTI/DKI values as the data that was not corrupted by signal drift for both ordered and randomized acquisitions. In the phantom acquisitions, the uncorrected randomized data yielded results closer to signal drift-corrected data than the uncorrected ordered data did. This was most pronounced for the MK and FA. Further, the residuals of the DKI fit (Fig. 4) are lower after signal drift correction in the randomized than the ordered datasets.

Improvements to Signal Drift Correction
To keep the proposed method applicable to as many studies as possible, we investigated only linear and quadratic correction, which could be reliably obtained from the low number of b ¼ 0-images in this study (around 10) or even less. The various causes of signal drift complicate an accurate a priori estimate of signal decrease patterns. Further, given the fact that the x-, y-, and z-gradients are irregularly switched on and off to provide diffusion sensitization, heating effects will likely be nontrivial. Despite those multiple unknowns, Figures  1, 6, and 7 show a consistently higher accuracy of corrections based on quadratic rather than linear fits. For longer scans or other manifestations of signal drift, the proposed quadratic correction may be suboptimal. A simple and data-driven way to improve on this would be to do a piece-wise linear correction, or a smoothed version thereof. A regular spacing of the b ¼ 0-images would likely be beneficial here (e.g., every tenth image), to ensure all parts of the signal drift are sampled equally.

Other Implications
Dynamic updating can help reduce signal drift, and it is, therefore, suggested to use this option, although this option might not be available on all clinical scanners or might increase scan time. Even with dynamic updating there is residual signal drift (Figs. 1-6). Our result show that (i) the simulations show signal drift-correction restores DTI/DKI measures to their "correct" values ( Fig.  5 and Supporting Information Fig. S6), and (ii) the signal drift-correction yields more stable global b ¼ 0-intensities than frequency-stabilization alone. Therefore, we believe that the proposed signal drift correction method is complementary to frequency stabilization and thus should be used even if frequency stabilization is on, and that the method is a legitimate alternative if frequency stabilization is not available or feasible.
In our data, we have roughly 110 DWIs over multiple shells to model 21 parameters (DTI þ DKI). Given perfectly randomly ordered gradient directions, both on the shells and within the scan duration, such an overdetermined system is estimated to be affected by the average signal drift over the session: the first, least-affected, DWIs could partially compensate for the last, most-affected, DWIs. For cases with more parameters and/or fewer measurements, such a compensation will not occur and the effects of signal drift will be more drastic.
Correction is especially relevant for multi-shell acquisitions on standard clinical systems. By default, these acquire data with b-values ordered from low to high. Not correcting for signal drift would mean that for higher b-values the signal attenuation caused by diffusion would be overestimated, which would correspond to an artefactual reduction in restricted diffusion. The expected effects on multi-shell diffusion methods such as CHARMED or NODDI (8,23) would in such instances mean that the restricted volumes would be underestimated. Moreover, because these models generally require more data to acquire resulting in longer scan times, higher b-values and thus a higher duty-cycle, and have a high number of parameters, the effects of signal drift will affect the fit more strongly.
The four phantom datasets without f 0 updating in Figure 1 (top two rows) show some variation in the magnitude of signal drift, even for the same acquisition schemes acquired on the same day (top vs. second row). The variation in signal drift is also present in vivo, ranging from 5.0 to 16.7% (Table 1), which was the maximum range of the multi-shell datasets in MASSIVE (10). This larger variation over dozens of sessions in the multisession project could possibly be due to a different load on the gradient system as a result of prescribed gradient orientations and b-value, a higher complexity of the imaged object (head vs. homogeneous phantom), or because of a larger inherent variation over different days.
The multi-shell acquisition protocols used in this study are representative of current research protocols, and as such provide a good indication of the magnitude of effects. For protocols that are less demanding of the scanner, for instance by having a lower b-value, effects will not necessarily be lower given the variety of causes for signal drift.
Given the multiple causes of signal drift, the variation between vendors, scanners, and scanner environment, and factoring in the large range of diffusion MRI scan protocols used across studies and sites, it is important to investigate the extent of signal drift and the impact on analyses for each site or scanner.
To demonstrate the widespread presence of signal drift in diffusion MRI data, and thus the importance of signal drift correction, we have examined datasets from the Human Connectome Project (HCP) [www.humanconnectome.org (24,25)]. This is regarded as one of the highest-quality neuroimaging databases to date, with data acquisition and processing being optimized by a collaboration between multiple institutions. Even in such high-quality data signal drift is present, as demonstrated in Supporting Information Figure  S7. This indicates that optimal use of HCP diffusion MRI data would also benefit from signal drift correction.

CONCLUSION
In this work, we have shown the detrimental effects of signal drift on diffusion MRI analyses that were present on all three investigated scanner vendors. The presence of signal drift can be easily examined when multiple b ¼ 0-images are acquired throughout the scan. From the same data, one can determine simple correction factors to mitigate the effects of signal loss on the subsequent analyses. Furthermore, the effect can be minimized by randomization of b-values throughout the acquisition. Whereas this is not a replacement for preventing the effects of signal drift to occur, it is essential to correct for the effects if they occur.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article. Figure S1. Mean signal of each volume of the acquired (left) and simulated (right) datasets. Figure S2 True simulated noise free signal (red) and the drift-affected noisy data (black). The left panel shows the data ordered as one would scan a randomized diffusion MRI acquisition; the right panel reordered by the signal (where the noise level was kept equal). Figure S3. Quadric fit of signal of known drift for different SNR levels. Figure S4. The left graph shows the difference in signal between the driftcorrected data and the simulated noisy data for each SNR level. The right graph shows the difference in signal between the drift-corrected data and the simulated noise-free data. The colors are indicated on the far right, with the colored area in each boxplot in the graphs indicating the 25 th -75 th percentile range, the white bar within the boxplot the median, and the error bars the 5 th -95 th percentile range. Figure S5. Phantom data acquired with randomized b-values (b 5 0, 1000, 2000, 3000, 9000 s/mm 2 ). The columns represent two repeats of the data with (left two columns) and without (right two columns) dynamic stabilization. The data was corrected for drift using quadratic fits based on data from each of the different b-values. Signal drift corrections estimated from the b 5 0, b 5 1000, b 5 2000, and b 5 3000 s/mm 2 data was identical, and only for the correction based on the b 5 9000 s/mm 2 there is a small (1-2%) overcorrection of the signal drift. Figure S6. Calculated DTI metrics for the simulation with FA 5 0.41 (left) and FA 5 0.81 (right). Different metrics are shown in each panel: the three eigenvalues of the diffusion tensor (k 1 , k 2 , k 3 ); the mean diffusivity (MD); fractional anisotropy (FA); and mean kurtosis (MK). Each panel shows data from the ordered and randomized acquisition, showing the data from the unaffected, drift-affected, and drift-corrected data in blue, red, and green, respectively. The colored boxplots indicate the 25 th -75 th percentiles, the white line within the box the median, and the error bars the 5 th -95 th percentile range. The solid grey line plotted behind the boxplots indicates the value for the data without signal drift; the dashed line indicates the value set in the simulations. Figure S7. Signal drift in ten example subjects from the HCP database.
The rows indicate ten separate subjects. The first column is the combined diffusion MRI dataset for this subject, as is regularly available for download as the "processed" dataset. This is the combination of the individual sessions of the second through last column, where the second and third column are the two sets of identical gradient directions but with opposite phase encoding direction. The fourth and fifth, and the sixth and seventh columns represent the second and third sets of gradient directions, respectively, each with opposite phase encoding direction. In each figure, the horizontal axis represents the number of the diffusion-weighted volume; the yaxis the normalized signal intensity). The red dots indicate the drift-affected signal, the black line the quadratic fit, and the blue dots the drift-corrected signal intensities. The dashed line at signal intensity 100 is shown as a visual reference for interpretation of signal drift.