Highly accelerated fMRI of awake behaving non-human primates via model-based dynamic off-resonance correction
Mo Shahdloo1, Daniel Papp2, Urs Schüffelgen1, Karla L. Miller2, Matthew Rushworth1, and Mark Chiew2
1Wellcome Centre for Integrative Neuroimaging, Department of Experimental Psychology, University of Oxford, Oxford, United Kingdom, 2Wellcome Centre for Integrative Neuroimaging, FMRIB, Nuffield Department of Clinical Neurosciences, University of Oxford, Oxford, United Kingdom


Dynamic B0 field inhomogeneities due to vigorous body motion - after mechanical head stabilisation - introduce severe artefacts in accelerated fMRI of awake behaving non-human primates (NHPs) by invalidating the calibration data used for unaliasing reconstruction. Here, we propose a method to estimate dynamic field perturbations via model-based frame-to-frame comparison of EPI reference navigators. These estimates can be used to improve reconstruction quality by matching each data frame to the calibration data, and simultaneously correcting the geometric distortions. The proposed method successfully estimates field perturbations and improves reconstruction quality in accelerated NHP fMRI, without the need for sequence modification or extra acquisitions.


Accurate unaliasing in accelerated image reconstruction often relies on consistency of the undersampled acquisition with some reference calibration data1. However, this assumption can be violated in a number of different ways, including by bulk motion, significant contrast changes, or fluctuating B0 field inhomogeneities. The last factor is an especially important source of artefacts in fMRI of awake non-human primates (NHPs), where significant unpredictable field inhomogeneities are caused by body movements that are likely to occur even if the head is fixed. Hence, estimating dynamic field perturbations and incorporating them in the reconstruction framework can reduce ghosting and residual aliasing artefacts that are not addressable in post processing.

Recently, multi-channel free induction decay navigators (FIDnavs) have been used to successfully estimate low‐spatial‐order dynamic field changes in fMRI2. However, utilising FIDnavs requires sequence modification, and generating an accurate forward model that depends on extra reference images that should be acquired with contrast parameters matched to the navigator, which might be difficult to achieve in some scenarios.

To mitigate these issues, we used reference EPI navigators to estimate dynamic field perturbations, without the need for contrast-matched reference images. These short navigators that sample the central line of k-space three times at each TR are already present in typical EPI sequences for correcting Nyquist ghost artefacts. Here, we used the spatial encoding provided by multi-channel navigators to estimate first order field perturbations in every TR. Field estimation quality is demonstrated in phantom and in vivo experiments, and improved reconstruction performance through significant reduction of residual aliasing is demonstrated using simulations of accelerated fMRI in NHPs.


As described in (3), $$$N_c\times N_c$$$ GRAPPA operators can be trained on calibration data to predict each k-space point using one adjacent sample on the Cartesian grid from $$$N_c$$$ channels. Letting ($$$G_x,G_y$$$) be the operators and ($$$\Delta_x,\Delta_y$$$) be grid spacing across x and y directions, arbitrary shifts of k-space data $$$S$$$ can be realised using fractional powers of the operators3:


Letting $$$s_j^n(t)$$$ be the channel $$$j$$$ navigator signal at frame $$$n$$$, signal at the reference frame ($$$n=0$$$) is expressed as:
where $$$\rho_j$$$ is the coil-sensitivity-encoded object magnetisation, and $$$\Delta\omega$$$ is the off-resonance. Field perturbations in NHP fMRI are mainly caused by motion in body parts that are distant from the imaging volume, and have low spatial frequency4. Hence, off-resonance at frame $$$n$$$ can be modelled as a linear spatial perturbation superimposed on the off-resonance at the reference frame:


yielding the navigator signal as:


Linear field perturbation is thus manifested as linear shifts in k-space:


where $$$S_j^n$$$ is the k-space data at frame $$$n$$$.

We modelled the shifted k-space data from the multi-channel three-line navigators using GRAPPA operators (illustrated in Fig.1):


and solved for $$$\alpha=b_x/\Delta_x$$$ and $$$\beta=b_y/\Delta_y$$$ using data from $$$N_c$$$ channels. Nelder-Mead simplex gradient-free nonlinear solver was used5. To improve numerical stability, only the five central k-space samples were used.

To examine the perturbation estimation accuracy, a bottle phantom was scanned on a 3T scanner using a 15-channel NHP receive coil and an EPI sequence with parameters: TE/TR= 30/2000ms,FA=90,FOV=192mm,1.5mm isotropic resolution. To mimic dynamic field perturbations, separate scans were performed where first-order shim terms (X,Y) were manually adjusted up to ±20μT/m in 5μT/m increments across scans.

To test the method in vivo, we acquired 50 frames of brain data from an awake macaque monkey using scan parameters matched to the phantom scan. We estimated field perturbation at each TR and distorted the fully-sampled reference using conjugate-phase reconstruction (CP)6. We then compared the distorted and acquired frames.

To demonstrate the effect of calibration inconsistency correction on reducing residual ghosting and unaliasing artefacts, data and navigators at the reference frame were distorted via CP, assuming a linear field perturbation (0.00±9.62Hz, mean±std across the FOV). Perturbation was estimated and data were corrected via inverse operators. Inconsistent and corrected frames were retrospectively undersampled at R=2,4 and GRAPPA-reconstructed. Multiband data were simulated with MB=2,R=3 by collapsing the data across slices and reconstructed using split-slice-GRAPPA7.


Figure 2 shows the shim terms estimated on the phantom. Modulated shim terms are accurately estimated, with absolute error of 0.56±0.10μT/m (mean±sem across frames).

Figure 3 shows the estimated distortions in vivo for two representative frames. The method accurately estimates the distortions in vivo.

Figure 4 shows the effect of field perturbation leading to calibration inconsistency, and correcting for it, on the accelerated reconstructions. Considerable ghosting artefacts arise from calibration inconsistency that cannot be corrected in post processing. The proposed method corrects the inconsistency, reduces artefact power, and corrects for geometric distortion in a single step.

Discussion and Conclusion

Calibration inconsistency due to dynamic field perturbations leads to significant aliasing artefacts when scanning awake behaving NHPs. To address this, we propose a method to estimate field perturbations using GRAPPA operators that does not require sequence modification or extra acquisitions. The method accurately predicts field perturbations in phantom and in vivo NHP data. We also demonstrate significant artefact reduction in synthesised accelerated EPI acquisitions by estimating and correcting the perturbations. The proposed improvements make highly accelerated awake behaving NHP imaging more robust and reliable, reducing the gap between what is possible with NHP protocols and state-of-the-art human imaging.


M.C. is supported by the Royal Academy of Engineering (RF201617\16\23). The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust (203139/Z/16/Z).


  1. Polimeni JR, Bhat H, Witzel T, Benner T, Feiweier T, Inati SJ, et al. Reducing sensitivity losses due to respiration and motion in accelerated echo planar imaging by reordering the autocalibration data acquisition. Magn Reson Med. 2016;75(2):665–79.
  2. Wallace TE, Polimeni JR, Stockmann JP, Hoge WS, Kober T, Warfield SK, et al. Dynamic distortion correction for functional MRI using FID navigators. Magn Reson Med. 2021;85(3):1294–307.
  3. Griswold MA, Blaimer M, Breuer F, Heidemann RM, Mueller M, Jakob PM. Parallel magnetic resonance imaging using the GRAPPA operator formalism. Magn Reson Med. 2005;54(6):1553–6.
  4. Pfeuffer J, Shmuel A, Keliris GA, Steudel T, Merkle H, Logothetis NK. Functional MR imaging in the awake monkey: effects of motion on dynamic off-resonance and processing strategies. Magn Reson Imaging. 2007;25(6):869–82.
  5. Nelder JA, Mead R. A Simplex Method for Function Minimization. Comput J. 1965;7(4):308–13.
  6. Man LC, Pauly JM, Macovski A. Multifrequency interpolation for fast off-resonance correction. Magn Reson Med. 1997;37(5):785–92.
  7. Cauley SF, Polimeni JR, Bhat H, Wald LL, Setsompop K. Interslice leakage artifact reduction technique for simultaneous multislice acquisitions. Magn Reson Med. 2014 Jan 1;72(1):93–102.


Figure 1. Linear magnetic field perturbation is cast as a linear shift of the k-space data for each navigator line at the reference frame (green dots), yielding shifted navigator lines (red dots). Separate GRAPPA operators (Gx, Gy) trained on calibration data to map each k-space point to their adjacent point across the orthogonal Cartesian grid can shift points on the grid (blue arrows). Partial shifts (red arrows) can be cast as fractional powers of the operators. Field perturbations were estimated by finding the shifts between navigators at each frame and those at the reference frame.

Figure 2. Dynamic field perturbations were introduced on the water phantom by manually modulating the linear shim terms across acquisitions. Perturbations were then estimated using the proposed method. Empty markers show the applied shim changes and filled markers show the estimates. The proposed method can accurately estimate the applied linear shim changes across (a) X, and (b) Y directions.

Figure 3. Field perturbations were estimated at each frame of the in vivo data. Distorted frames were then synthesised using the reference frame and the estimates. One representative frame exhibits distortion in the posterior parts (frame 20), while the other suffers in the anterior parts (frame 41). Actual distorted frames (left) and synthesised frames (middle) are shown. Red contour shows the object outline at the reference frame. Right column shows the error between the synthesised and actual frames. The proposed method can accurately estimate the field perturbations in vivo.

Figure 4. (a) Distortion due to a simulated linear field perturbation was applied on the fully-sampled reference frame from in vivo data. (b) In-plane accelerated data were synthesised by retrospectively undersampling at R=2, 4. Multiband data were synthesised by adding the data across two slices (R=3x2). Field perturbation leads to calibration inconsistency, manifested as strong unaliasing artefacts. The perturbation was estimated and the data were corrected, yielding a reconstructions that have significant reduced unaliasing artefacts and corrected geometric distortion.

Proc. Intl. Soc. Mag. Reson. Med. 29 (2021)