Susceptibility-induced distortion that varies due to motion: Correction in diffusion MR without acquiring additional data

Because of their low bandwidth in the phase-encode (PE) direction, the susceptibility-induced off-resonance field causes distortions in echo planar imaging (EPI) images. It is therefore crucial to correct for susceptibility-induced distortions when performing diffusion studies using EPI. The susceptibility-induced field is caused by the object (head) disrupting the field and it is typically assumed that it remains constant within a framework defined by the object, (i.e. it follows the object as it moves in the scanner). However, this is only approximately true. When a non-spherical object rotates around an axis other than that parallel with the magnetic flux (the z-axis) it changes the way it disrupts the field, leading to different distortions. Hence, if using a single field to correct for distortions there will be residual distortions in the volumes where the object orientation is substantially different to that when the field was measured. In this paper we present a post-processing method for estimating the field as it changes with motion during the course of an experiment. It only requires a single measured field and knowledge of the orientation of the subject when that field was acquired. The volume-to-volume changes of the field as a consequence of subject movement are estimated directly from the diffusion data without the need for any additional or special acquisitions. It uses a generative model that predicts how each volume would look predicated on field change and inverts that model to yield an estimate of the field changes. It has been validated on both simulations and experimental data. The results show that we are able to track the field with high accuracy and that we are able to correct the data for the adverse effects of the changing field.


Introduction
When an object, such as the human head, is placed in a homogeneous magnetic field it will disrupt that field and it will no longer be homogeneous. This disruption, the difference between the initial homogeneous field and the actual field, is known as a susceptibility-induced offresonance field.
The off-resonance field presents a particular problem for echo-planar imaging because of its low bandwidth in the phase-encode direction (Jezzard and Clare, 1999). It is not unusual to have a bandwidth per pixel of~10 Hz for EPI images acquired without in-plane acceleration (IPAT). Typical off-resonance values for the worst affected areas (near sinuses and ear canals) are~150 Hz for a 3T scanner, which leads to signal displacements of up to 15 pixels. This must be resolved and typical solutions aim at finding the offresonance field (field map) either through direct measurement (Jezzard and Balaban, 1995;Robson et al., 1997) or by acquiring at least two EPI images that are affected differently by the off-resonance field (Andersson et al., 2003;Morgan et al., 2004;Holland et al., 2010;Irfanoglu et al., 2015). Once the field is known it is easy to correct for the first order effects of the distortions (Jezzard and Balaban, 1995;Munger et al., 2000).
As a first approximation the off-resonance field is stationary in the object (the head) framework. That is to say, if the subject moves to the right inside the scanner, the off-resonance field will also move to the right. But this is not true of any rotations around an axis non-parallel to the magnetic flux (the z-direction). This means that EPI images acquired with the subject in different orientations will be distorted in slightly different ways, and even after a rigid-body alignment they will not line up perfectly (Andersson et al., 2001;Ward et al., 2002;Sulikowska, 2016). To correct such a set of images with a single field map will not address the issue of differentially distorted images.
A solution to the problem would entail finding a unique field map for each EPI volume. Andersson et al. (2001) attempts to calculate the fields directly from the observed shape differences remaining in an fMRI time-series after rigid-body motion correction. Others (Boegle et al., 2010) have suggested that if a map of magnetic susceptibilities for the object is known the field can be calculated for each observed orientation of the object. Another suggestion is to use the phase information in the complex images of an fMRI time-series to get an idea of the changing field (Ooi et al., 2013;Hutton et al., 2013;Dymerska et al., 2016). Finally, Alhamud et al. (2016) has suggested acquiring a dual echo-time volumetric navigator to get a low resolution fieldmap for each volume.
However, all of these suggestions have limitations, in particular if one wants to apply them to diffusion EPI images rather than fMRI. It is not trivial to measure or deduce magnetic susceptibilities for an object, diffusion sensitisation interferes with the phase information in the images and techniques requiring non-standard acquisitions are not available on most scanners. As a consequence there is presently no widely used method to correct diffusion data for a dynamically changing susceptibility-induced field.
In this paper we present an extension of the work in Andersson et al. (2001) and . It uses the same Taylor expansion of the susceptibility field as Andersson et al. (2001), but with a generative model for diffusion data Andersson and Sotiropoulos (2015). We have validated this new method using a combination of simulations  and experimental data.

Theory
How to represent the dynamically changing field Let us denote the susceptibility induced field in a framework that is fixed with respect to the object (the head) by ωðxÞ where x denotes any coordinate in the object framework. This field will depend on the orientation of the object with respect to the flux of the external field, so it is really ωðxÞ ¼ ωðx : rÞ (1) which we can think of as ωðxÞ given the object location r within the field, where r denotes the parameters of a rigid body transform. It implies that we would have to estimate a 3D field for any point r in a six-dimensional space, which at a first glance seems intractable. However, as we outlined in the introduction, the field is only expected to change as a consequence of out-of-plane rotations. Hence, we can reduce r to ½θ ϕ, i.e. any object rotation around the x-or y-axes. Furthermore, Sulikowska (2016) have shown that provided that θ and ϕ vary within a reasonably small range (small in this case being AE10 degrees) the field changes linearly with θ and ϕ. This means that we may be able to approximate it with a Taylor expansion around some point r 0 ¼ ½θ 0 ϕ 0 at which we know ωðxÞ.
where N is the order of the expansion and where ∂ω ∂θ is a spatial map of the same size as ωðxÞ and signifies the partial derivate of the field, at each point x, with respect to θ at the point ½θ 0 ϕ 0 (and similarly for ∂ω=∂ϕj r0 , ∂ 2 ω=∂θ 2 j r0 etc).
For the remainder of the theory section we will limit the derivation to a first order Taylor expansion of the field where it is assumed that only rotations around the x-and y-axes affects the field. In appendix B we show how this is extended to higher order and more movement parameters.
Equation (2) is demonstrated in an intuitive fashion in Fig. 1.

The model for correcting the images
We have previously described a model that we use to correct the observed diffusion images when the susceptibility-and eddy currentinduced off-resonance fields as well as the subject movement have been estimated where where b s i is the corrected volume i and f i is the original volume i. J x denotes the Jacobian determinant of the transformation at location x. R i is a rigid-body transformation for the ith volume and d x is a voxel displacement vector for location x. ω and eðβ i ; r i Þ are scalar off-resonance fields caused by susceptibility and eddy currents respectively. β i is a vector of parameters encoding the eddy currents for the ith volume and r i is a vector of movement parameters specifying the orientation of the object in that volume. a i are the acquisition parameters of the ith volume and serves to translate the scalar off-resonance field ω þ e into a vectorvalued voxel-displacement value d. For many protocols a i would have the same value for all i, but in for example the HCP (Sotiropoulos et al., 2013) protocol it would take one of two values and in the dHCP (Hutter et al., 2017) protocol one of four values. A more detailed version of the spatial model used inside eddy can be found in . For the purpose of this discourse it is useful to divide the transform into on the one hand rigid-body movement and eddy current-induced offresonance and on the other hand susceptibility-induced off-resonance. To that end we redefine x 0 as so that we can rewrite equation (4) where ω is the susceptibility field in Hz, t i denotes the total readout time for a slice of volume i and p x i denotes the x-component of the phaseencoding of volume i (and correspondingly for p y i ). Typically only one of p x i and p y i are non-zero for a given i and will have a value of 1 or -1.
Using the expression for ω from equation (2) we can rewrite equation (7) as where Δθ i and Δϕ i are the rotations around the x-and y-axes respectively of volume i compared to the reference r 0 orientation.

Adding Gaussian Process predictions
For the following we will assume the existence of an independent prediction s i of how the corrected image b s i from equation (8) should look. This is obtained using a Gaussian Process (GP) and is a linear combination of all the diffusion weighted images (DWIs) corrected for the current estimates of distortions and movements. It has been described in detail in Sotiropoulos (2015, 2016). Assuming that s i % b s i we can write where ε is the "error" or the deviation between s i and b s i . For the case where ∂ω=∂θj r0 , ∂ω=∂ϕj r0 etc are unknown (or not fully known, like in the early stages of the iterative estimation) that error signal will contain information about ∂ω=∂θj r0 , ∂ω=∂ϕj r0 etc.

Estimating the derivative fields
The formal derivation of the equations for estimating the field is in appendix A. Below we present an intuitive explanation of the method.

Intuitive explanation of our proposed method
Our explanation depends on two notions. The first is that there is a small number of "maps" that capture the variability of the off-resonance field over the course of an experiment. The second notion is that the residual variance, i.e. the difference between observations and predictions, contains enough information to allow us to estimate these "maps".
These maps represent one anatomical location per voxel, as defined by the reference (first) volume, and specifies the linear change with respect to the parameter in question (pitch or roll) at that location. Specifically is it the slope of a linear relationship between the strength of the off-resonance field and orientation parameter (pitch or roll) in units of Hz per degree. Formally these maps are maps of the partial derivative of the off-resonance field with respect to pitch and roll. Fig. 2 illustrates the forward model that is used. It shows how the observed differences are explained in terms of off-resonance changes caused by changes in pitch and changes in roll. Furthermore it demonstrates how these can be subdivided into changes in deflection of the signal and changes in stretching/compression of the signal, both along the PE-direction (in the figure the y-direction is the PE-direction).
To get an intuition for this consider, for example, the part that is given by Δθ 6 ð∂f 6 =∂yÞ ð∂ω=∂θÞ halfway down on Fig. 2. The map ∂ω=∂θ denotes the rate of change of the field with respect to pitch in units of Hz/ degree, so when multiplying this with Δθ 6 (degrees) we get a map with the change in field (in Hz) as a consequence of pitch change for volume # 6. In this forward model it is assumed that the derivative fields (∂ω=∂θ and ∂ω=∂ϕ) are known, but these are of course the maps we are interested in estimating. So, the next step is to describe a model to estimate the unknown entities ∂ω=∂θ and ∂ω=∂ϕ.   (2) and shows how the field ω i for any volume i is approximated by a linear combination of a measured field ω 0 and the derivative of that field with respect to θ (rotation around the x-axis) and ϕ (rotation around the y-axis). The weights for the derivative fields are given by the estimated movement parameters where Δθ i denotes the rotation around the x-axis of volume i relative the orientation that ω 0 was acquired in (and correspondingly for Δϕ i ). The maps in the figure are estimated from our simulations and the grey-scales are À40 to 100 Hz for the w i fields and À5 to 5 Hz/degree for the ∂ω=∂θ and ∂ω=∂ϕ fields. An intuitive description of the ∂ω=∂θ field is "How much the field changes if one nods forward (looks down) one degree". The corresponding description for ∂ω=∂ϕ would be "How much the field changes if one tilts one's head to the right (in the direction of the dark part of the field) one degree.
The matrix B contains a B-spline basis set that covers the entire 3D space of a data volume. The vector b consists of two "sub-vectors", b θ and b ϕ , that contain the weights for the B-spline basis such that the resulting weighted sums are the estimated fields.
In the formulation in Fig. 3 s À f is a very tall (Nn x n y n z Â 1) column vector, that is modeled as the product of a matrix X calculated from It demonstrates how the total difference can be described as the sum of the signal changes caused by a change in pitch and the signal changes caused by a change in roll. The signal changes can be further subdivided into changes in translation of the signal along the PE-direction (i.e. how the sampling point in the predicted image changes) and changes in Jacobian modulation (intensity changes caused by stretching/compression). The choice of volumes 1, 6 and 31 for demonstration is arbitrary and the model of course encompasses all volumes (both b ¼ 0 and DWI volumes). The symbol was used to denote Hadamard (or elementwise) product.
T contains the weights for the fields ∂ω=∂θ and ∂ω=∂ϕ such that B b b θ is an estimate of ∂ω=∂θ. B y denotes a matrix that is organised in the same way as B, but where the columns consists of splines that have been differentiated in the y-direction (PE-direction).
known entities and an unknown vector b. This is simply an over determined system of linear equations that can be solved in a least squares sense. However, the forward model described in Fig. 2 is a linearisation of a problem that is fundamentally non-linear with respect to the unknowns (b). Therefore solving the linear system in Fig. 3 represents the linearisation in a Gauss-Newton scheme and has to be repeated iteratively, each time updating b such that where X and f both depend on b ðkÞ . In appendix A we provide the formal derivation of the Gauss-Newton step in equation (10) and that is the proper reference for the proposed method. The calculation of X T X in equation (10) above is potentially very computationally expensive. The reason for that, and how to reduce some of that expense, is explained in appendix C.

Combining the model with estimation of eddy currents
There is an extra caveat with the method outlined in the previous sections when there are also eddy current-induced distortions, as explained in this section.
Looking at the derivative fields in Fig. 1 it is clear that there are strong linear components in both ∂ω=∂θ and ∂ω=∂ϕ. This means that if one were to model the change of the field for those images as just a linear gradient (in the y-direction for ∂ω=∂θ and in the x-direction for ∂ω=∂ϕ) it would already explain a fair proportion of the variance seen in the data (this is the idea behind some attempts at real-time shimming (Ward et al., 2002;Alhamud et al., 2016)). And for diffusion weighted images we estimate a low order eddy current-induced field for each volume, hence this will include also the low (spatial) order parts of the susceptibility-by-movement interaction.
In itself it is not necessarily a problem that the EC-correction subsumes the low frequency component of the derivative field. One would expect the estimated susceptibility derivative field to be the high frequency component of the true derivative field, and that the superposition of the change of the susceptibility and the estimated EC field would still be correct. However, the b ¼ 0 volumes do not suffer from EC-induced distortions, so none is modeled. That means that the derivative field one needs for the b ¼ 0 volumes is the "full field", i.e. both the low and high spatial frequencies. That leads to problems if one wants to model a joint derivative field for the b ¼ 0 and the diffusion weighted images alike. The option to model separate fields for b ¼ 0 and diffusion volumes would be suboptimal for several reasons. It would increase computational effort, but above all it would negatively affect the precision of the estimated fields and potentially introduce shape differences between the corrected b ¼ 0 and diffusion volumes.
We have therefore opted for a scheme where we attempt to estimate the "full field" (low and high spatial frequencies) for b ¼ 0 and diffusion weighted volumes alike. It is done through a two step process where an initial set of EC parameters is estimated while limiting the estimates by a linear causal model of the diffusion weighting (referred to as a "linear second level model" in ). This limits the risk that the EC-model attempts to "explain" the shape differences that has been caused by the susceptibility field varying with subject orientation. However, this can potentially lead to suboptimal results because i) the linear second level model may be too strict and give slightly suboptimal results  and ii) the EC-model can still "explain" some susceptibility shape differences. For these two reasons, a second step is employed where iterations for estimating the susceptibility derivative fields (i.e. equation (A6)) are interleaved with iterations for estimating EC-parameters (without a second level model), subject movement and outliers Andersson et al., , 2017. A schematic of the algorithm is as follows:

Combining the model with estimation of intra-volume movement
We expect this method to be useful primarily for data/populations where there is a relatively high degree of movement. Therefore it is likely that the data will also be marred by intra-volume movement, i.e. there will be volumes where the subject has moved during the acquisition of the slices/MB-groups that constitute the volume. This results in a stack of slices that may be translated and/or rotated relative each other and may not even be parallel. It is important to ensure that the two corrections (intra-volume and susceptibility-by-movement) work seamlessly together.
As can be seen in Figs. 2 and 3 our model for estimating the derivative fields operates with a volume as the unit. When the algorithm (see algorithm 1 above) reaches the stage where it estimates the derivative fields, one set of intra-volume movement parameters have already been estimated and used to correct the data that goes into the estimation of the derivative fields. Hence, the problem is not primarily that volumes contaminated by intra-volume movement "confuse" the estimation of the fields. The problem is that volumes where the subject has performed a large movement (in particular if that includes a pitch and/or roll) have not been acquired under the influence of a well defined field.
We have solved that by examining the (slice-wise) movement parameters for each volume after the initial step (prior to the estimation of the derivative fields) to identify all volumes with an intra-volume standard deviation of Δθ and/or Δϕ greater than one degree. These volumes are considered "poorly defined" with respect to orientation, and they are excluded from the estimation of the derivative fields. However, it should be noted that they are included in the movement-eddy-current-outlier estimations that are interleaved with the derivative field estimation. That way these parameters can still be updated for those volumes in the light of the derivative fields. It should also be noted that when building the field for use in correcting a given volume, we use the intra-volume parameters (provided a slice-to-volume model was used in the estimation) so that the Δθ and Δϕ that are actually used in equation (2) vary from slice-to-slice/group-to-group.
In the simulations described below the method was tested with both a volumetric and a slice-to-volume movement model. For the experimental data with large subject movements it was inevitable that there would also be large intra-volume movement and only the slice-to-volume model was used to analyse those.

Implementation
The method described in the present paper has been implemented in Cþþ as part of the eddy tool  and http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/EDDY) in FSL (see Smith et al. (2004) and http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/). Some of the steps of the algorithm, in particular the calculation of the derivatives and the spline-interpolation on an irregular grid, are very computationally intense. They have therefore been parallelised for NVidia GPUs using CUDA (Farber, 2011). All the results in the present paper have been generated using the GPU version of the software.
The current version of eddy will estimate parameters for eddy current-induced off-resonance fields and subject movement (see equations (4) and (5) and ). In addition it will also detect and replace movement induced signal dropout  and correct intra-volume movement . The work presented here has extended that to also include susceptibility-by-movement effects.
The new addition has increased the execution time for eddy. Not so much because of the actual estimation of the fields, but more for the additional iterations for estimation of eddy currents and movement that are interleaved with the field estimation. Hence, it effectively doubles the execution time compared to running eddy without the susceptibility-bymovement correction.

Simulations Data
The simulations that we use are similar to those described in Graham et al. (2017) where the authors characterise the problem that the method in the present paper is attempting to solve. In brief, the data was simulated using POSSUM, a highly realistic MRI simulator (Drobnjak et al., 2006) that has been extended to diffusion imaging . It contains an input object (head) that has been segmented into air and tissue and from which a susceptibility induced-field can be calculated for any orientation relative to the magnetic flux .
Using POSSUM, a diffusion data set was simulated with 60 b ¼ 1000 volumes and five interspersed b ¼ 0 volumes, where the diffusion directions were optimised on the whole sphere and were identical to those used for the FMRIB standard diffusion protocol. The simulated acquisitions were single-band with a 72 Â 96 Â 55 matrix, an interleaved slice ordering and an isotropic voxel size of 2.5 mm. The total readout time for a slice was 77.4 ms, corresponding to a bandwidth of 12.9 Hz/pixel in the PE-direction. The susceptibility field was calculated for a 3T static field and the eddy currents (when applicable) were commensurate with a Stejskal-Tanner diffusion encoding.
Various other options were varied according to a 2 Â 2 Â 2 Â 2 Â 2 factorial. The factors that were varied were Movement type: Volumetric or continuous movement. We expect this method to be of particular use for subjects that move more than average. Therefore we simulated data both with the movement injected "per volume" (a commonly used model) and with movement continuous in time (a more realistic model (Andersson et al., 2017)). The latter will result in a subset of "broken" volumes with movement between consecutive slices such that when stacked up they no longer represent anatomically faithful volumes.
Movement magnitude: Normal or large movement. The "large" movement trace was taken from a volunteer moving according to an agreed pattern including pitch and roll (rotation around the x-and ydirections). The "normal" movement was created from the same trace, but with all parameters divided by three.
Eddy current distortions: With or without eddy current-induced distortions.
Phase-encode direction: With phase-encode direction A→P or P→A. SNR: With an SNR of 20 or 40. The SNR is defined from the b ¼ 0 volumes as the average intra-cerebral signal divided by the noise.
Ten distinct simulations with unique noise were performed for each bin in this hyper-cube, resulting in a total of 320 data sets each consisting of 65 vol.
In addition to the simulations described above we also simulated data with exactly the same protocol but without any movement or offresonance fields. These "ground truth" data sets were simulated without noise, or with added noise to yield b ¼ 0 SNR of 20 and 40.

Analysis
If nothing else is specified the analyses below was performed using a first order Taylor expansion around x-and y-rotation, a spline knotspacing of 10 mm, 20 iterations and a regularisation λ of 10 (see equation (A13)). The static field was estimated from a pair of A→P and P→A b0-volumes (using the FSL-tool topup, Andersson et al. (2003)) from the simulations. These were the first volumes from their respective data sets, so before any movement was injected. Movement was estimated either using a volumetric model or with a slice-to-volume model with 17 degrees of freedom per volume.
Estimating the "true" derivative fields. The susceptibility fields used for the simulations were generated in a framework that was fixed with respect to the object (head). Therefore these were all reoriented to the space of the first volume (the reference space for eddy) using the known movement parameters, after which the phase was fitted to a general linear model consisting of a constant and the (known) rotations around the x-and yaxes (Δθ and Δϕ respectively). The resulting parameter maps for Δθ and Δϕ are considered as the "true" ∂ω=∂θ and ∂ω=∂ϕ fields.
Checking the assumptions of the Taylor expansion. One of the assumptions of the model is that a first order Taylor expansion is sufficient to track the changes of the field with head orientation. This was tested on the simulated data by generating R 2 -maps for a first order model. First the "true" first order derivative fields were calculated, as described above, after which the fraction of variance explained by these was calculated.
Comparing estimated and true derivative fields. The "true" derivative fields were estimated as described above and compared to the ∂ω=∂θ and ∂ω=∂ϕ maps estimated using the method described in the present paper.
Comparing estimated and true off-resonance fields. For each volume we had full knowledge of the actual off-resonance field that had affected that acquisition. We were able to compare that, on a volume-by-volume basis, to the estimated field. The estimated field was calculated in two ways, as The sum of the static susceptibility field estimated with topup and the eddy current-induced field estimated by eddy. The sum of the static susceptibility field (topup), the eddy currentinduced field and the change in the susceptibility induced-field caused by the change of orientation of the head. The two latter estimated by eddy.
This analysis was limited to the data simulated with volumetric movement because the slice-to-volume case does not have fields that are defined for entire volumes.
Comparing corrected and true images. After correcting the simulated data for movement and distortions it was compared to the "ground truth" data simulated without any movement or off-resonance effects. This was done by calculating the correlation between paired (with the same diffusion encoding) volumes. We compared the results when excluding or including susceptibility-by-movement effects in the correction.
The effects of the diffusion encoding direction (being fixed in the scanner framework) changing in the subject framework as the subject rotates its head (Leemans and Jones, 2009) were not included in any of our simulations. Hence, the volume pairs should have identical contrast.

Human data
A healthy and experienced volunteer was scanned three times in the same session, twice with the FMRIB standard diffusion protocol and once with an augmented protocol. The standard protocol consists of 60 unique directions with a b-value of 1500 and five interspersed b ¼ 0 volumes. The augmented protocol had the same 60 directions with a b-value of 1500, but instead of five b ¼ 0 volumes it had a total of 61 b0-volumes such that the full sequence was "b0, b0, DWI, b0, DWI, b0 … b0, DWI". Each of these scans was preceded by the acquisition of two b ¼ 0 volumes with identical scanning parameters except for the phase-encode direction being reversed. For the first scan the subject was instructed to lie as still as possible. For the second scan he was instructed to make a set of movements (rotations around the x-and y-axis), where after each movement the new position was held for~10 volumes. For the third scan the instruction was to try and mimic the movement from the second scan, but this time hold each new position for~20 volumes. Note that, similarly to the simulations with intra-volume movement injected, this means that the human data contained volumes contaminated by intra-volume movement.
The data was acquired on a Siemens Magnetom Prisma system with a 32-channel receive head coil and monopolar Stejskal-Tanner gradients were used for the diffusion encoding. The resolution was 1.5 mm isotropic with 84 slices of a 150 Â 150 matrix. A multi-band factor of 4 was used, no in-plane acceleration, phase-encoding A→P and a 6/8 partial k-space. The TR was 3.35 s, the TE was 88.8 ms and the total readout-time (time between centre of first and last echo) was 107 ms, yielding a bandwidth per voxel of 9.32 Hz.

Analysis
We performed several analyses of the data to explore different aspects of the method. If nothing else is stated all estimations/corrections below were performed with a slice-to-volume movement model with 17 degrees of freedom per volume, eddy current correction and detection and replacement of outliers. The susceptibility-by-movement estimation was performed using a knot-spacing of 10 mm, a regularisation λ of 10 and 20 iterations in 5 interleaves with movement, EC and outlier estimation as outlined in algorithm 1.
Test-retest reproducibility of the estimated derivative fields. By excising relevant b0-volumes from the third scan we mimicked the data from the second scan. Thus we had two corresponding, but independent, data sets with similar levels of movement from which to estimate the ∂ω=∂θ and the ∂ω=∂ϕ fields. In order to ensure that we performed the Taylor expansion around the same point in orientation space we used the static field estimated from scan two for both scans. We also replaced the first b ¼ 0 in the third scan by the first b ¼ 0 from the second scan, thus ensuring the same reference space for both analyses.
Ability to estimate derivative fields from diffusion data. It is potentially more difficult to estimate the derivative fields from diffusion data than from data with inherently similar contrast. The differing contrast/information content in the different diffusion weighted volumes is one complication, and there is also the difficulty of estimating it in the presence of volume specific eddy current-induced distortions. To investigate that, the third scan was sub-sampled into two data sets, one where the additional b ¼ 0 volumes had been stripped away to yield a data set consisting of 60 diffusion weighted volumes and 5 interspersed b ¼ 0 volumes as in the standard protocol. The other sample consisted only of the 61 b ¼ 0 volumes. For this analysis the static field was estimated from the first b ¼ 0 volume from the third scan and the immediately preceding b ¼ 0 scan with reversed PE-direction. The two samples have the same initial volume, so the reference space was the same. The estimation from b ¼ 0 data was performed in the same way as for the diffusion weighted data, except that eddy currents were not estimated.
Ability to correct the data for susceptibility-by-movement effects. For this analysis the images acquired in the first scan were considered as "ground truth". Both visual inspection and the estimated movement parameters confirmed that there were very little movement and no obvious movement related dropout. It was corrected using a volume-to-volume movement model, eddy current correction and detection and replacement of outliers. Data set two and three (the latter after stripping of extraneous b ¼ 0 volumes) were used as "test" data sets. They were corrected as described above, once with and once without susceptibilityby-movement estimation. The fields estimated from the reversed PEdirection data were used as the static fields. After within-scan correction a rigid-body motion transform was calculated for scans 2 and 3 relative to scan 1, and the former data sets were resampled into the space of scan 1 using spline interpolation. This entails an additional interpolation, but that should not matter for the comparison between results with and without susceptibility-by-movement estimation.
Unlike the simulations, these data will be affected by the diffusion encoding direction rotating in the subject framework when the subject rotates its head. However, that should only lead to a slightly lower correlation for the volumes with the greatest rotation and will not bias the comparison between results with and without susceptibility-bymovement.

Simulations
Looking at movies of the simulated data ( Fig. S1 in the supplementary material) after correcting for eddy currents, subject movement and a static susceptibility field it was clear that in the areas with a strong offresonance field from susceptibility it was not sufficient to correct with a single field. The susceptibility-by-movement effects were obvious as apparent shape changes between volumes. It was also clear that these were very well corrected by the susceptibility-by-movement model. Fig. 4 shows the R 2 -maps for two different levels of the brain. It can be seen that the maps are close to 1 for most of the intra-cerebral voxels. In fact they are greater than 0.99 for most of the brain at these levels. The regions that are less than 0.99 coincide with areas where it was problematic to estimate the "true" field due to phase wrapping. Together with the observation that there are fewer voxels with an R 2 less than 0.99 for the large movement case makes us conclude that a first order model is sufficient for the simulations. If those areas had truly reflected a need for a higher order expansion one would have expected more voxels with R 2 less than 0.99 for the large compared to the normal movement case.

Checking the assumptions of the Taylor expansion
Comparing estimated and true derivative fields Fig. 5 shows estimated versus "true" derivative fields for some of the simulations. The particular simulations used for Fig. 5 used a volumetric movement model, phase-encoding in the A→P direction and an SNR of 40. It can be seen that the estimated derivative fields capture the features of the true fields very well for both "normal" and "large" movements and regardless of the presence of eddy current-induced distortions.
The correlations between "true" and estimated derivative fields are tabulated for all the different simulations in Tables S1 and S2 in the supplementary material S. Fom the results in Fig. 5 and Tables S1 and S2 we conclude that Fig. 4. This figure shows the R 2 -maps that demonstrate the proportion of true variance in the simulated off-resonance fields explained by a first order model. The left column shows b ¼ 0 images for anatomical guidance. The middle and rightmost columns show the R 2 -maps for the "normal" and "large" movement simulations respectively. The greyscale goes from 0 (black) to 1 (white). rows respectively) when using simulations with volumetric movement, phase-encoding A→P and an SNR of 40. The first column shows the true b ¼ 0 images for anatomical guidance. Columns 2-4 show the results for the "normal" movement case. Column 2 shows the "truth", column 3 the estimated fields for the "No eddy currents" case and column 4 the estimated fields when they were jointly estimated with eddy current-induced fields. Correspondingly columns 5-7 show the results for the "large" movement case, where column 5 is the "truth" and columns 6 and 7 shows the estimated fields in the absence and presence of eddy currentinduced fields respectively. The units of the colorbars are Hz/degree.
We are able to accurately estimate the rate of change of the field from the diffusion data. That ability is largely independent of PE-direction, the absence or presence of eddy-currents, the absence or presence of intra-volume movement artefacts and of SNR in the range 20-40. The one things that appears to have an appreciable effect on that ability is the magnitude of movement. We believe that is partly explained by a poorer "estimate of truth" (can be seen as ringing artefacts in Fig. 5) for the normal movement case. But also that it reflects less information about the susceptibility-by-movement in data with less movement.

Comparing estimated and true off-resonance fields
The results from the comparison of estimated and true off-resonance fields are shown in Fig. 6. It shows the simulations for phase-encode direction A→P and an SNR of 40. The results for phase-encode direction P→A and SNR 20 are sufficiently similar that we have chosen to not include those results. As expected it shows that the advantage of our method is greater for "large" than for "normal" movement, though there is a clear advantage also for "normal" movement. The loss of accuracy of the estimated fields when not taking susceptibility-by-movement into account (grey lines) is strongly related to a "rotation proxy" (dashed black line) that we calculated as a composite measure of the x-and yrotations (see legend to Fig. 6). It can also be seen that the accuracy of the estimated fields were poorer when the simulated data contained eddy current-induced distortions. However, the relative increased accuracy offered by the method in the present paper remained the same.

Comparing corrected and true images
Figs. 7 and 8 show the correlation between true and corrected images for a volumetric and slice-to-volume movement model respectively. The correlation was calculated for a single slice corresponding to the most basal slice shown in Figs. 4 and 5, i.e. a slice that is appreciably affected by susceptibility-by-movement effects. These (Figs. 7 and 8) correspond to the simulations with phase-encoding in the A→P direction and an SNR of 40. Similar figures for both both SNR levels (20 and 40) and for slices at different levels (in the inferior-superior directions) of the brain are shown in Figs. S3-S34 in the supplementary material. The appearance of the plots in Figs. 7 and 8 are similar to those in Fig. 6, which simply shows that there is a close relationship between getting the fields right and getting the corrections right.

Human data
Similarly to the simulated data, the human data also exhibited clear residual distortions after correcting for eddy currents, subject movement and a static susceptibility field. This manifested itself as changes in distortions between volumes, and can be seen as a movie in S2 in the supplementary material. It can also be seen how these apparent shape changes largely disappear when correcting also for susceptibility-by- Fig. 6. This figure shows the correlation between true and estimated off-resonance fields for all volumes of the simulated data. The data used for this figure was simulated such that the phase encode direction was A→P and the SNR was 40. The solid black and solid grey lines represent the correlation between true and estimated off-resonance fields for the method in the present paper and the method assuming a constant susceptibility field respectively. The scale for both these curves is found on the left-hand of the graphs. The dashed black line shows a "proxy" for the rotation relevant to susceptibility field. It was calculated as where R x and R y denote rotation around the x-and y-axes respectively. The scale for that curve is found on the right hand side of the graphs. The graphs in the left column show the results for "normal" subject movement and the right column for "large" movement. The top row shows the situation when no eddy current distortions were simulated and no attempt was made to estimate eddy currents. The bottom row shows the situation when eddy currents were included in the simulations and eddy currents and susceptibility-by-movement were jointly estimated by eddy. movement effects.
Test-retest reproducibility of the estimated derivative fields The outcome of estimating the derivative fields from two different runs is shown in Fig. 9. It can be seen that the reproducibility is good and the correlation between the two ∂ω=∂θ fields is 0.945 and for the ∂ω=∂ϕ fields it is 0.880. It can also be seen that there is a general agreement with the appearance of the derivative fields estimated from the simulated data (see Fig. 5). When considering Fig. 9 it should be remembered that the fields are estimated from completely distinct data sets with distinct subject movements and with a slightly different starting position of the subject.

Ability to estimate derivative fields from diffusion data
The fields estimated from a diffusion data set (60 DWI þ 5 b ¼ 0) compared to a field estimated from only b ¼ 0 volumes (61 vol) are shown in Fig. 10. The agreement is very good with a correlation of 0.966 for the ∂ω=∂θ fields and 0.975 for the ∂ω=∂ϕ fields.
Ability to correct the data for susceptibility-by-movement effects Fig. 11 shows a similar pattern to the simulated data with respect to our ability to correct the data. The greater our "rotation proxy" the greater the difference between the correction with and without susceptibility-by-movement correction. The overall correlations after correction are lower than for the simulated data. This would be expected since the rotation of the diffusion gradients in a subject framework leads to increasing differences in contrast with increasing subject rotation. We also find that the measured fields tend to be "sharper" (i.e. contain higher spatial frequencies) than the simulated ones, making their correction a little more difficult.

Discussion
In this paper we have demonstrated our ability to track the movement-induced changes in the susceptibility-induced off-resonance field. It is very important to be clear that this is done without any need to actually measure any field, either by traditional field measurement or by reversed gradients. The method achieves this purely from the deviation of the observed images from the expected/predicted images. This is encouraging since it will enable workers to correct their data for this artefact without the need to alter protocols, acquire additional data or invest in costly hardware. It also means that the method can be applied to "historical" data that suffers from high levels of subject movement.
It should be noted that if a single, static, field is measured (for example by a volume pair with reversed gradients at the start of the scan) the expansion in equation (2) is performed with that field as the constant. In that case the corrected data will be internally consistent (corrected for susceptibility-by-movement effects) and anatomically faithful (corrected for susceptibility). If no susceptibility field is provided the constant term in equation (2) will be zero. The data will still be corrected for Fig. 7. This figure shows the correlation between true and corrected images for all volumes of the simulated data. The data used for this figure was simulated such that the phase encode direction was A→P, the SNR was 40 and all movement was "intervolume", i.e. any movement was assumed to occur between acquisition of consecutive volumes. Correspondingly the analysis used a volumetric movement model. The solid black and solid grey lines represent the correlation between true and corrected images for the method in the present paper and the method assuming a constant susceptibility field respectively. See the legend for Fig. 6 for more details.
susceptibility-by-movement effects, and hence be internally consistent, but will not be anatomically faithful.
The results from our simulations and from our human experimental data are in good agreement, and indicate that the errors due to susceptibility-by-movement effects are determined by the amount of pitch and/or roll relative to the orientation of the subject when the field was measured. They also indicate that the method in the present paper can to a large degree eliminate this error. The rate-of-change maps are also in good agreement with those reported by Sulikowska (2016), both in magnitude and general appearance. Note that Sulikowska (2016) appears to have chosen an opposite convention with regards to the sign of the rotation (our fields are predicated on positive Δθ and Δϕ corresponding to nodding forward and tilting one's head to the right respectively), and also differ in terminology such that what we denote "roll" is denoted "yaw" in Sulikowska (2016).

When do we need this?
The effects of the susceptibility field changes are relatively small in most of the brain and depend strongly on how much the subject moves. It is therefore of greatest importance for subjects who are uncooperative, such as babies, children or patients that find it hard to remain still. In this respect it is similar to corrections for intra-volume movement , though the way its effects manifest themselves is very different. Intra-volume movement introduces variance in the data that is, as a first approximation, proportional to the temporal derivative of the movement. The changing susceptibility, in contrast, causes variance that is proportional to the movement, specifically to rotation around any axes non-parallel with the magnetic flux.
Our assumption is therefore that the present method will be important, together with the intra-volume motion correction , for all categories of subjects that move a lot.
In addition we predict that it will be useful for very long examinations (such as for example the diffusion acquisitions in the HCP project, Sotiropoulos et al. (2013)), where even slow changes in position can add up to sizeable movements over the course of an hour or more. Since the susceptibility-by-movement effects are proportional to relative orientation it doesn't matter how slow the movement is and it can potentially have a significant impact on such studies.

Higher b-values
In the present paper we have tested the method on simulated data with b-values up to 1000 and on human data with b-values up to 1500. It has previously been demonstrated  that eddy can correct for eddy current-induced distortions and subject movement for high b-value data. That ability hinges on its internal predictive model for the diffusion weighted images. We therefore expect this method to similarly work well on data up to at least the b-values of 5000 that was used in . Fig. 8. This figure shows the correlation between true and corrected images for all volumes of the simulated data. The data used for this figure was simulated such that the phase encode direction was A→P, the SNR was 40 and movement was continuous, i.e. volumes occurring during periods of rapid movement were corrupted by intravolume movement. Correspondingly the analysis used a slice-to-volume movement model. The solid black and solid grey lines represent the correlation between true and corrected images for the method in the present paper and the method assuming a constant susceptibility field respectively. See the legend for Fig. 6 for more details.

Possible extensions
Another source of changing field is the breathing of the subject. Even though the lungs are a relatively long distance away from the object we image (the head) its effects on the field are experienced as a very smoothly varying off-resonance field at the location of the head. Such a smoothly varying field is in principle easy to compensate for with shimming, but the problem is that it varies during the respiratory cycle. Hence, it is another source of a field varying over time.
The framework we have suggested in the present paper can be extended to take this source into consideration. One would then replace (or add to) Δθ and Δϕ with a low-order basis expansion of the point in the respiratory cycle (for example measured using a bellows belt). This could be for example p, cosðπpÞ and cosð2πpÞ where p is a number between 0 and 1 defining the point in the respiratory cycle. It would also need to be changed so as to allow for different values of p, cosðπpÞ and cosð2πpÞ for different slices/multi-band-groups in Fig. 3 and equation (A8).
Our method could in principle also be combined with prospective motion correction (Maclaren et al., 2013) where subject movement is estimated by some external device (Maclaren et al., 2012) or by a 3D navigator (Alhamud et al., 2012) and used to inform the prescription of subsequent acquisitions. In that case the external recordings of pitch and roll would be used in lieu of the retrospectively estimated ones in equations (2), (8) and (A8). It would also mean that the algorithm as described in 1 would be simplified in that no movement parameters would need to be estimated.
Is it a problem that the derivative fields are poorly estimated when subjects move little?
No, it is not. It is the variance that is introduced by the susceptibilityby-movement that drives the estimation of the derivative fields. In the absence of susceptibility-by-movement-induced variance there is nothing to drive the estimation and it will be dominated by the regularisation of the fields, resulting in smooth estimated fields with values close to zero. When used for correction of the data (equation (2)) these "close to zero" fields will be multiplied by very small values for Δθ and Δϕ, and effectively not contribute at all. Which is the desired outcome when there is no susceptibility-by-movement-induced variance.

Alternative methods
As far as we know the only viable alternatives to the method suggested in the present paper are the dynamic shimming methods suggested by Ward et al. (2002) and by Alhamud et al. (2016). The former of those methods was implemented for gradient echo EPI and has as far as Fig. 9. This figure shows the test-retest agreement of derivative fields estimated from two different data sets in the same subject. The first column is shown for anatomical reference and the two rightmost columns show the estimated derivative fields. The top two rows show the rate-of-change of the field with respect to pitch, and the bottom two rows with respect to roll. The units of the fields are Hz/degree rotation.
we know not been applied to diffusion data. Furthermore, those methods are only capable of correcting for field changes that are constant or linear in space. As can be seen from Figs. 5, 9 and 10 that is not sufficient to fully capture the actual field changes due to subject movement. Our tests (data not shown) indicate that a linear combination of linear fields in the principal directions can capture roughly half of the observed field changes. We don't know if either of those methods are freely available or implemented by any vendor.
The method suggested by Boegle et al. (2010) hinges on having a map of susceptibilities in exact register with the data one wishes to correct. Neither the acquisition or registration of such a map is trivial.
Methods based on tracking the phase of each voxel in the time-series (Ooi et al., 2013;Hutton et al., 2013;Dymerska et al., 2016) can not be used with spin-echo EPI data.

Relation to earlier work
The ideas and principles behind the suggested method are similar to those described in Andersson et al. (2001), but there are many important differences.
The inclusion of a predictive model for diffusion weighted images has allowed us to extend the method to diffusion data.
The co-existence of volume-specific eddy current-induced fields with a changing susceptibility-induced field presents additional challenges that have been solved. It is implemented within a framework that models and corrects for intra-volume movement. The forward model has been complemented with an additional term (described as "Signal change caused by change in stretch/compression of signal" in Fig. 2) compared to the model that was used in Andersson et al. (2001). The use of a spline basis set (as opposed to the DCT set used in Andersson et al. (2001)) and the re-ordering of summations described in appendix C have facilitated estimating the derivative fields at higher spatial resolution with realistic computational and memory demands.
Where to next?
In a series of papers we have described the development of a comprehensive tool to correct diffusion data for some of the most important sources of artefacts. The current version (as described in the present paper) corrects for eddy current-induced distortions, susceptibility-distortions (including changes with subject orientation) and subject movement (inter-and intra-volume bulk motion as well as Fig. 10. This figure shows the agreement between derivative fields estimated from b ¼ 0 volumes only (middle column) and from a diffusion data set (right column). The first column contains a corresponding b ¼ 0 volume after correction for susceptibility distortions, and is shown for anatomical reference. The top two rows show the rate-ofchange of the field with respect to pitch, and the bottom two rows with respect to roll. The units of the fields are Hz/ degree rotation. motion induced signal drop-out). Additional sources that we hope will be addressed in future work include Spin history: The signal in a given voxel depends on the time that has passed since it was previously excited, with shorter time leading to less signal. In the absence of movement this time is identical for all voxels in all volumes. When there is movement, that time can vary considerably and cause unwanted variance over volumes. With knowledge of the history of excitation for each voxel (which requires an intra-volume movement model) it should be possible to correct.
Receive bias field: There is, especially for coil arrays, an inhomogeneous intensity bias field caused by differential sensitivity to the RFsignal from different parts of the FOV. This is, to a first approximation, fixed in the scanner framework. When the subject moves within that, stationary, field it will result in unwanted variance across the volumes.
Long time-constant eddy currents: We have seen some evidence that the first one or two slices/MB-groups for a given volume (diffusion gradient) are affected by different eddy currents compared to the remaining slices/MB-groups. This violates the assumption within eddy of a single, constant, eddy current-induced field per volume.

Conclusion
We have augmented our framework for simultaneous correction of susceptibility-and eddy current-induced distortions and subject movement effects with a model that estimates how the susceptibility-induced field changes with subject orientation. fields) copies of the Jacobian volume J i concatenated horizontally.