Estimation and removal of movement-by-susceptibility induced variance in fMRI time series.
- Previous Work
- Assumptions behind Unwarp
- The "Realign & Unwarp" function
- Installation instructions
Even after realignment there is considerable variance in fMRI time series that covary with, and is most probably caused by, subject movements. It is also the case that this variance is typically large compared to experimentally induced variance. Anyone interested can include the estimated movement parameters as covariates in the design matrix, and take a look at an F-contrast encompassing those columns. It can be quite dramatic. The result is loss of sensitivity, and if movements are correlated to tas specificity. I.e. we may mistake movement induced variance for true activations. Because the movement induced variance is often very large compared to "true" activations false positives may ensue even from a relatively modest correlation between task and movement.
The problem is well known, and several solutions have been suggested. A quite pragmatic (and conservative) solution is to include the estimated movement parameters (and possibly squared) as covariates in the design matrix. Since we typically have loads of degrees of freedom in fMRI we can usually afford this. The problem occurs when movements are correlated with the task, since the strategy above will discard "good" and "bad" variance alike.
The "covariate" strategy described above was predicated on a model where variance was assumed to be caused by "spin history" effects, but will work pretty much equally good/bad regardless of what the true underlying cause is.
Others have assumed that the residual variance is caused mainly by errors caused by the interpolation kernel in the resampling step of the realignment. One has tried to solve this through higher order resampling (huge Sinc kernels, or k-space resampling).
Assumptions behind Unwarp
The "Unwarp" toolbox is based on a different hypothesis regarding the residual variance. EPI images are not particularly faithful reproductions of the object, and in particular there are severe geometric distortions in regions where there is an air-tissue interface (e.g. orbitofronal cortex and the anterior medial temporal lobes). In these areas in particular the observed image is a severely warped version of reality, much like a funny mirror at a fair ground. When one moves in front of such a mirror ones image will distort in different ways and ones head may change from very elongated to seriously flattened. If we were to take digital snapshots of the reflection at these different positions it is rather obvious that realignment will not suffice to bring them into a common space.
The situation is similar with EPI images, and an image collected for a given subject position will not be identical to that collected at another. Hence, even after a "successful" realignment there will be residual variance caused by the object having different shape at different time points. We call this effect susceptibility-by-movement interaction. "Unwarp" is predicated on the assumption that the susceptibility-by-movement interaction is responsible for a sizeable part of residual movem nt related variance.
The Realign & Unwarp function
Assume that we know how the deformations change when the subject changes position (i.e. we know the derivatives of the deformations with respect to subject position). That means that for a given time series and a given set of subject movements we should be able to predict the "shape changes" in the object and the ensuing variance in the time series. It also means that, in principle, we should be able to formulate the inverse problem, i.e. given the observed variance (after realignment) and known (esti ated) movements we should be able to estimate how deformations change with subject movement.
We have made an attempt at formulating such an inverse model, and at solving for the "derivative fields". A deformation field can be thought of as little vectors at each position in space showing how that particular location has been deflected. A "derivative field" is then the rate of change of those vectors with respect to subject movement. Given these "derivative fields" we should be able to remove the variance caused by the susceptibility-by-movement interaction. Since the underlying model is so re tricted we would also expect experimentally induced variance to be preserved. Our experiments have also shown this to be true. Indeed one particular experiment even indicated that in some cases the method will reintroduce experimental variance that had been obliterated by movement related variance.
In theory it should be possible to estimate also the "static" deformation field, yielding an unwarped (to some true geometry) version of the time series. In practice that doesn't really seem to work, hence the method deals only with residual movement related variance induced by the susceptibility-by-movement interaction.
Combined use with FieldMap
This means that the time-series will be undistorted to some "average distortion" state rather than to the true geometry. If one wants additionally to address the issue of anatomical fidelity one should combine Unwarp with a measured field-map. A field-map in the format that Unwarp expects can be created using the FieldMap toolbox.
The description above can be thought of in terms of a Taylor expansion of the field as a function of subject movement. Unwarp alone will estimate the first (and optionally second, see below) order terms of this expansion. It cannot estimate the zeroth order term (the distortions common to all scans in the time series) since that doesn't introduce (almost) any variance in the time series. The measured fieldmap takes the role of the zeroth order term. Refer to the FieldMap toolbox and the documents FieldMap.man and FieldMap_principles.man for a description of how to obtain fieldmaps in the format expected by Unwarp.
What derivatives do we need to model?
If we think of the field as a function of subject movement it should in principle be a function of six variables since rigid body movement has six degrees of freedom. However, the physics of the problem tells us that the field should not depend on translations nor on rotation in a plane perpendicular to the magnetic flux. Hence it should in principle be sufficient to model the field as a function of out-of-plane rotations (i.e. pitch and roll). One can object to this in terms of the effects of shimming (object no longer immersed in a homogenous field) that introduces a dependence on all movement parameters. In addition SPM/Unwarp cannot really tell if the transversal slices it is being passed are really perpendicular to the flux or not. In practice it turns out thought that it is never (at least we haven't seen any case) necessarry to include more than Pitch and Roll. This is probably because the individual movement parameters are typically highly correlated anyway, which in turn is probably because most heads that we scan are attached to a neck around which rotations occurr.
On the subject of Taylor expansion we should mention that there is the option to use a second-order expansion (through the defaults) interface. This implies estimating also the rate-of-change w.r.t. to some movement parameter of the rate-of-change of the field w.r.t. some movement parameter (colloquially known as a second derivative). It can be quite intresting to watch (and it is amazing that it is possible) but rarely helpful/necessarry.
Re-estimation of movement parameters
In addition to inducing residual (after realignment) movement-related variance, movement-by-susceptibility distortion changes may bias the estimation of the movement. For each iteration Unwarp gets a better idea of the true shape of each scan, and can potentially get a better estimate of movement. In Unwarp there is an option (default) to re-estimate (do a new realign) the movements between each iteration of estimating the fields. Our testing indicates that this is a good idea.
Jacobian intensity modulation
In the defaults there is also an option to include Jacobian intensity modulation when estimating the fields. "Jacobian intensity modulation" refers to the dilution/concentration of intensity that ensue as a consequence of the distortions. Think of a semi-transparent coloured rubber sheet that you hold against a white background. If you stretch a part of the sheet (induce distortions) you will see the colour fading in that particular area (you can also think of future appearance of a tattoo obtained when young and slim). In theory it is a brilliant idea to include also these effects when estimating the field (see e.g. Andersson et al, NeuroImage 20:870-888). In practice for this specific problem it is NOT a good idea.
Should I use it?
It should be noted that this is a method intended to correct data afflicted by a particular problem. If there is little movement in your data to begin with this method will do you no good. If on the other hand there is appreciable movement in your data (>1mm or >1deg) it will remove some of that unwanted variance. If, in addition, movements are task related it will do so without removing all your "true" activations.
The method attempts to minimise total (across the image volume) variance in the data set. It should be realised that while (for small movements) a rather limited portion of the total variance is removed, the susceptibility-by-movement interaction effects are quite localised to "problem" areas. Hence, for a subset of voxels in e.g. frontal-medial and orbitofronal cortices and parts of the temporal lobes the reduction can be quite dramatic (>90%).
It should also be noted that there are several reasons for residual movement related variance. Notably:
- Susceptibility-distortion-by-movement interaction: Is what we deal with here.
- Susceptibility-dropout-by-movement interaction: The susceptibility induced field inhomogenity will, in addition to distortions, cause signal loss due to through-plane dephasing (which will not be rephased by encoding gradients that are all in-plane). This signal modulation will also change with subject position.
- Spin-history effects: The signal will depend on how much of longitudinal magnetisation has recovered (through T1 relaxation) since it was last excited (short TR→low signal). Assume we have 42 slices, a TR of 4.2seconds and that there is a subject z-translation in the direction of increasing slice # between one excitation and the next. This means that for that one scan there will be an effective TR of 4.3seconds, which means that intensity will increase.
- Slice-to-vol effects: The rigid-body model that is used by most motion-correction (e.g. SPM) methods assume that the subject remains perfectly still for the duration of one scan (a few seconds) and that any movement will occurr in the few μs/ms while the scanner is preparing for next volume. Needless to say that is not true, and will lead to further apparent shape changes.
Unwarp will deal only with the first of these, which means that even after that correction the other components will remain. It is difficult to say which of these effetcs dominate. It will depend on ones scanner, and even on the specific data set.
Friston KJ, Williams SR, Howard R, Frackowiak RSJ and Turner R (1995) Movement-related effect in fMRI time-series. Magn Reson Med 35:346-355
Jezzard P and Balaban RS (1995) Correction for geometric distortions in echoplanar images from B0 field variations. Magn Reson Med 34:65-73
Wu DH, Lewin JS and Duerk JL (1997) Inadequacy of motion correction algorithms in functional MRI: Role of susceptibility-induced artefacts. J Magn Reson Imag 7:365-370
Hutton C, Bork A, Josephs O, Deichmann R, Ashburner J and Turner R (2002) Image distortion correction in fMRI: a quantitative evaluation. NeuroImage 16:217-240.
About method behind Unwarp Toolbox
Andersson JLR, Hutton C, Ashburner J, Turner R, Friston K (2001) Modelling geometric deformations in EPI time series. NeuroImage 13:90-919
The software was developed in part while JLRA was supported by STINT (Stiftelsen for internationalisering av h�gre utbildning och forskning.)
This software has been tested by the authors on data collected on two different make MR scanners. Undoubtedly there are numerous bugs that will surface when more people start using it. We are grateful for bug reports, and other input, and the more we receive the faster bugs will be weeded out.
Unwarp Toolbox Distribution
Unwarp Toolbox for SPM5 can be downloaded with the latest version of SPM5
Since Unwarp is a part of main SPM2 this distribution will (eventually, when bugs are fewer) be part of the general SPM updates. Meanwhile we recommend the following method for installing the update.
- Create a local directory, e.g. /usr/local/'whatever'/Unwarp_updates
- Download the Unwarp_ddmmyy.tar.gz file into that directory
- Unzip and untar, on a Linux machine
..]$ gunzip Unwarp_ddmmyy.tar.gz
..]$ tar xvf Unwarp_ddmmyy.tar.gz
- If you are on a Macintosh, or some other exotic brand, you will
need to compile spm_dilate_erode.c. It should work to simply cd into your Unwarp_updates
directory, start Matlab and at the Matlab prompt type
>> mex spm_dilate_erode.c
- Make sure that /usr/local/'whatever'/Unwarp_updates is included in your Matlab path ahead of your main SPM distribution.
- If (very optional) you want the help in the SPM graphical help (for the "Realign & Unwarp" button) to reflect Realign and Unwarp you should rename mod_spm.m spm.m and put it in your main SPM distribution directory. DONT rename it and leave it in the Unwarp_updates directory!
Jesper Andersson (jesper.andersson at ks.se)
Chloe Hutton (chutton at fil.ion.ucl.ac.uk)
John Ashburner (john at fil.ion.ucl.ac.uk)
Robert Turner (r.turner at fil.ion.ucl.ac.uk)
Karl Friston (karl at fil.ion.ucl.ac.uk)