Estimation and Removal of Residual Motion Artifact in Retrospectively Motion-Corrected fMRI Data: A Comparison of Intervolume and Intravolume Motion Using Gold Standard Simulated Motion Data

INTRODUCTION

Functional MRI (fMRI) and resting-state fMRI (rs-fMRI) play a major role in the assessment of functional activity and connectivity in the brain.^1,2 With these techniques, a time series of blood oxygenation level–dependent (BOLD) contrast images is collected throughout the brain over a 5-10 minute interval, and the temporal BOLD contrast change at each voxel is then compared to the external stimulus timing or other voxels. Correlation between the BOLD contrasts and its spatial similarity patterns are studied for their relation to neuronal activation or the formation of functional connectivity (FC) networks. Because a fast image acquisition technique is desirable to investigate whole brain activation and networks, multislice 2D echo planar imaging (EPI) scans are commonly used in rs-fMRI and fMRI studies to obtain temporal sampling typically of 1-2 s.^3–6

While rs-fMRI and fMRI scans are designed to detect neuronal signal fluctuations, in practice non-neuronal signals are also included in the data. The effective removal of these latter contributions is critical to producing unbiased results. Early research found that head motion artifact is one of the main sources of bias in fMRI studies.⁷ Additionally, physiologic signal fluctuation has also been found to affect MRI findings.^8–10

Head motion during rs-fMRI and fMRI acquisition can lead to alterations in the sensitivity of the radiofrequency (RF) transmitter/receiver, the spatial encoding process, and B0 field modulation; each of these contributes to various artifacts.^11–13 In addition to the local contrast modulation and image distortion from the B0 field alterations, other head motion artifacts also need to be addressed when a 2D EPI sequence is used. For instance, B0 field fluctuations from breathing patterns induced phase encoding (PE) direction image motion in 2D EPI acquisitions,¹⁴ and a different scale of PE direction image shift is reflected in each slice. Head motion also causes protons to shift between slices, altering the time between RF excitations. This causes the steady state of magnetization of each slice to be permuted, an effect known as the altered spin excitation history effect.^15,16

Previous fMRI studies have attempted to correct for head motion prospectively by altering the gradient after the motion is estimated through the use of an optical device,¹⁷ a navigator in the sequence,¹⁸ or real-time image volume registration.^19–22 However, by far the most common method for correcting head motion is retrospective 3D rigid volume motion correction, since motion-tracking equipment and sequences are not available at most sites. An algorithm for 3D rigid volume motion correction for a time series of 2D EPI scans has been proposed^23–25 and is commonly implemented in rs-fMRI and fMRI analysis using software such as AFNI,²⁶ BrainSuite,²⁷ FSL,²⁸ and SPM.²⁹

Despite these advances, residual motion artifact has been reported even after perfect motion correction.³⁰ The source of this artifact is known to be the partial volume (PV) effect of surrounding voxels due to the resampling of the target image, aligned to the reference resamples. To remove the residual motion artifact, volume motion corrected images are commonly “regressed out” using nuisance motion regressors, and volumes deemed to be excessively corrupted are censored (or “scrubbed”). Various motion nuisance regressors have been proposed. For instance, early studies described the use of position displacement to predict residual motion artifact.^31,32 Other researchers have proposed reducing residual motion artifact through the use of a voxel-wise combination of displacements and their previous time-shift and squared displacements³³ or through the use of nonlinear expansion of the realignment parameters, including temporal derivatives and corresponding squared regressors.³⁰ Although a number of studies have described the use of multiple motion nuisance regressors to reduce residual motion artifact,^34–42 it should be noted that increasing the number of nuisance regressors reduces the degrees of freedom (DOF) in the fMRI dataset. Ideally, the smallest possible number of nuisance regressors should be used.^43,44

To minimize the number of motion nuisance regressors required, the mechanism of the residual motion artifact must be clearly understood, and only these necessary motion nuisance regressor(s) should be employed. However, this is not feasible without motion-corrupted MRI data in which intravolume motion and the spin excitation history effect are emulated. Previous work has proposed that data-driven nuisance regressors can represent the unwanted non-neuronal signals including the residual motion artifact in fMRI data using various techniques, such as: principle component analysis,⁴⁵ independent component analysis,^46–51 the average signal in white matter and cerebrospinal fluid regions of interest (ROIs),^52,53 or local non-gray matter signal, as with ANATICOR.⁵⁴ However, most data-driven motion nuisance regressors are estimated volume-wise, in which it is assumed that head motion is restricted between volumes and the motion within the volume acquisition or slice-wise motion is negligible.

Intravolume motion correction, or the map-slice-to-volume approach, was first proposed for use in 2D EPI fMRI⁵⁵ and was later expanded to address the spin excitation history effect¹⁶ and geometric distortion.⁵⁶ A method to interpolate intervolume motion at each slice in the context of slice acquisition timing was then proposed.⁵⁷ This was followed by the proposal of an intravolume motion correction approach to measure in-plane and out-of-plane motion separately in each slice, which was validated in a cadaver study using the simulated prospective acquisition correction (SIMPACE) sequence.⁵⁸ A more recent study presented a different method for reducing intravolume motion artifact⁵⁹; with this technique, the intravolume motion is not directly measured or corrected, but the residual motion artifact reflected in each slice is instead removed in the context of the slice acquisition timing between volume acquisitions.

This study extends previous work regarding the development of the slice-oriented motion correction method (SLOMOCO).⁵⁸ We collected ex vivo brain phantom data using the SIMPACE sequence with various types of intervolume/intravolume injected motion. Using SIMPACE data, we sought to investigate the cause of residual motion signal after motion correction. Based on these findings, we identify the most effective motion nuisance regressors and suggest how to improve the SLOMOCO pipeline using these regressors.

We note that the current SLOMOCO software package and its pipeline are now available via GitHub (https://github.com/wanyongshinccf/SLOMOCO). This offers improved version control features and update functionality. Additionally, users are able to raise issues directly through the website, to report bugs or unexpected behavior. Developers in the neuroimaging community are also able to submit pull requests.

METHODS

Ex vivo brain phantom

An ex vivo brain phantom was designed using a previously described technique.⁶⁰ The phantom was fixed with formalin and then soaked in Fomblin (FluidX, Salt Lake City, UT). Bubbles in the tissue and vascular structures were removed using a vibration machine and vacuum pump. Floating debris in the container was removed, and the brain phantom was then positioned in a 3D-printed holder inside an insulated container sized 17.78 cm (top) × 15.24 cm (bottom) × 26.67 cm (height) (Coleman Company, B003363V3Li, IL, USA). The phantom inside the container was then placed in a 20-channel head and neck coil (Nova Medical, Wilmington, MA) as shown in Fig. 1A.

Figure 1.A) Diagram of ex vivo brain phantom positioning in the scanner and 2D imaging plane. B) Injected motion in DICOM coordinates. A: anterior; P: posterior; L: left; R: right.

SIMPACE sequence

The SIMPACE sequence used in this study has been described previously.⁵⁸ Using the prospective acquisition correction (PACE) feature proposed by Thesen et al.,¹⁹ we altered the imaging plane of each slice from a time series of 2D EPI acquisitions with user-defined information. Using the SIMPACE sequence, we then injected intervolume (volume-wise) and/or intravolume (slice-wise) motion, generating realistic motion-corrupted EPI data (as shown in a movie file in Supplementary Fig. 1).

It should be noted that SIMPACE synthesizes the motion-corrupted MR data by altering the imaging plane before each slice and volume acquisition, resulting in the emulation of intervolume/intravolume motion, respectively. Since the position of the brain phantom is fixed, additional motion artifact from altered B0 and B1 inhomogeneity effects due to motion is not modeled. In addition, the use of ex vivo brains means that physiologic noise and its related motion on in vivo EPI images are not included in these data. However, they still remain highly useful test cases for understanding important features and consequences of motion in fMRI datasets.

MRI scanning

The brain phantom was scanned on a 3T scanner (Siemens Healthineers, Erlangen, Germany) using the SIMPACE sequence. Single-band SIMPACE data were acquired using the following parameters: repetition time (TR)/echo time (TE) = 2 s/30 ms, voxel size = 2 × 2 × 4 mm³, matrix size = 96 × 96, flip angle = 80°, ascending interleaved slice acquisition, and 21 slices. To replicate the positioning of an in vivo fMRI study, the EPI scan was conducted in the coronal imaging plane in the setting of a head-first supine and the coordinates in DICOM were manually changed to the axial plane after data collection. The injected and estimated motion was described in DICOM coordinates, and rotation was applied in a counterclockwise direction (see Fig. 1). The first 4 volumes of SIMPACE data were removed to obtain steady-state signals for analysis.

Injected motion patterns

Instantaneous intervolume motion

Intervolume motion was injected during a given TR, with a half voxel shift (2 mm) and a whole voxel shift (4 mm) in the z-directions at the 10^th volume. The 11^th volume was then moved back to the initial position. A total of 19 volumes were acquired (Fig. 2A). In the same manner, 4° of x-, y-, and z-rotation motion was injected during a given TR.

Linearly increasing shift intervolume motion

Linearly increasing shift and rotation motion was injected into the x-, y-, and z-axes over 30 volumes (60 s), with up to 4 mm of shift and 4° of rotation at the 10^th volume (Fig. 2B).

Figure 2.An example of injected shift motion after the removal of the first 4 volumes. (A) 4 mm of instant intervolume shift motion at the 6^th volume. (B) Linearly increasing intervolume motion up to 4 mm during 30 volumes. White dots indicate the acquisition timing of the first slices at each volume, with a gap of 2 s of TR.

Intervolume motion from in vivo datasets

Five individuals were scanned twice in a 3T scanner using the Young Adult Human Connectome Project fMRI protocol.⁶¹ This was approved by the local institutional review board, and all participants provided informed consent. The following parameters were used for these scans: TR/TE = 0.8 s/37 ms, voxel size = 2 × 2 × 2 mm³, multiband (MB) acceleration = 8, 72 slices, and 420 volumes. From 10 fMRI datasets, a time series of rigid volume motion was measured using AFNI’s 3dvolreg (Supplementary Fig. 2), and the estimated rigid volume motion parameters with a sampling rate of 1/0.8 Hz were temporally interpolated (down-sampled) with a half sampling rate to synthesize intervolume motion-corrupted SIMPACE data with TR = 2 s.

Intravolume motion from in vivo datasets

Estimated rigid volume motion parameters from the 10 in vivo fMRI scans described in the previous section were temporally interpolated with a sampling rate of 10.5 Hz (= 21slices divided by a TR of 2 s). Note that temporally up-sampled motion is continuous and therefore may not represent realistic spontaneous motion within volume acquisitions. To address this, we amplified the intravolume motion by calculating the intervolume and intravolume motion separately, amplifying the intravolume motion by a factor of 2, and adding the amplified intravolume motion to the intervolume motion (Fig. 3). In summary, 10 intervolume/intravolume motion-injected SIMPACE datasets and 10 SIMPACE datasets with identical intervolume motion but with 2× amplified intravolume motion were generated.

Figure 3.An example of injected in vivo z-rotation motion in SIMPACE data across (A) the full time series duration and (B) a zoomed-in segment of the time series. The black and blue lines indicate the intervolume motion with a step-wise change every 2 s, and intervolume and intravolume motion with 21/2 Hz sampling rate. The red line represents 2× amplified intravolume motion.

Voxel-wise PV regressor

Because we observed a repeated residual motion artifact pattern in the SIMPACE data with injected linearly increasing intervolume motion (see Results), we proposed the use of a voxel-wise PV regressor. First, static images were generated with the reference image concatenated. The measured intervolume motion was then injected inversely onto each static image to simulate the motion-corrupted data, as previously described by Wilke⁶² and referred to as MOTSIM by Patriat et al.⁴⁵ Intervolume motion was then aligned back to the reference images with measured motion, synthesizing the time series of the residual artifact with the intervolume motion both forward and in reverse. This represented the resampling bias or PV effect due to the alignment (Fig. 4).

Figure 4.Generation of the voxel-wise partial volume (PV) regressor. (Step 1) Static time series images were generated with the reference image (ref) concatenated. (Step 2) The estimated volume motion was injected inversely (inv) onto each static image to simulate the motion-corrupted data (Motsim). (Step 3) Injected motion was then removed with previous estimated motion parameters, from which voxel-wise time series of residual motion artifact could be estimated by resampling both forward and in reverse.

Validation of PV regressor

We compared the residual time series after intervolume motion correction with a PV regressor on the 10 in vivo intervolume motion-injected SIMPACE datasets. The proposed PV regressor was tested as a motion nuisance regressor using the various motion correction pipelines. Six rigid volume motion parameters (6 Vol-mopa) and their derivatives (6 Deriv) were also tested as motion nuisance regressors.

In the 10 intervolume motion-injected SIMPACE datasets, the first 4 volumes were removed to obtain steady-state signals, and 3D rigid volume motion correction (VOLMOCO) was then performed based on the first volume using AFNI’s 3dvolreg. The residual time series processed using various combinations of motion regressors (in addition to linear detrending terms): 6 Vol-mopa (VOLMOCO + 6 Vol-mopa), PV regressor (VOLMOCO + PV), 6 Vol-mopa + PV regressor (VOLMOCO+ 6 Vol-mopa + PV), and 6 Vol-mopa and 6 Deriv (VOLMOCO + 6 Vol-mopa + 6 Deriv). The regression was conducted using AFNI’s 3dREMLfit. The voxel-wise residual sum of the squared (SoS) time series signal and standard deviation ($\sqrt{\frac{SOS}{vol}}$, where vol is the volume number) were then calculated. Averaged SD values in the gray matter (GM) were reported for each model. A GM mask was generated in native EPI space using SPM software²⁹ with a GM probability higher than 0.9. The top and bottom slices of the mask map were not included so as to exclude voxels with excessive out-of-plane motion.

Modified SLOMOCO (mSLOMOCO) method and its motion nuisance regressor model

We present a modified version of the original SLOMOCO (oSLOMOCO) method.⁵⁸ First, this modified SLOMOCO (mSLOMOCO) defines the individual reference for each target slice in an improved manner. After each volume was aligned to the initial reference volume, the affine transformation matrix was stored. A new reference volume was defined for each volume using the inverse affine transformation matrix. In-plane and out-of-plane motion calculations for each slice were identical to those used in oSLOMOCO.

The major change in the mSLOMOCO pipeline involved adding a PV motion nuisance regressor. While the oSLOMOCO pipeline employed 14 voxel-wise displacements and their derivative terms with time delay as motion nuisance regressors, the proposed mSLOMOCO pipeline included 6 volume-wise Vol-mopa, 6 slice-wise rigid slice motion nuisance parameters (Sli-mopa), and 1 voxel-wise PV regressor as motion nuisance regressors. Figure 5 displays a workflow of mSLOMOCO, with the modifications from previous work highlighted with yellow boxes. Volume-wise and slice-wise motion estimation and correction were conducted using AFNI commands, and motion nuisance components were regressed out using the lscov function in MATLAB.⁶³

Figure 5.Modified slice-oriented motion correction (mSLOMOCO) workflow. Yellow-colored steps are newly added in the mSLOMOCO pipeline. VOLMOCO: rigid volume motion correction; Ref. per vol: reference for each volume; moco: motion correction; Vox: voxel-wise; PV: partial volume; Vol: volume; Sli: slice; mopa: motion parameters.

Evaluation of mSLOMOCO and its motion nuisance regressor models

mSLOMOCO was evaluated with each of 10 in vivo intervolume/intravolume motion-injected SIMPACE datasets with 1× and 2× amplified intravolume motion. To investigate the effects and relative properties of the correction, the following combinations of terms were applied in parallel ordinary least squares models: without motion nuisance regressors (mSLOMOCO); with 6 Vol-mopa after SLOMOCO (mSLOMOCO + 6 Vol-mopa); and with 6 Vol-mopa, 6 Sli-mopa, and voxel-wise PV regressor (mSLOMOCO + 12 Vol-/Sli-mopa + PV). Note that the proposed mSLOMOCO pipeline employs 12 Vol/Sli-mopa + voxel-wise PV regressors as motion nuisance regressors. The different combinations of 6 Vol-mopa, 6 Sli-mopa, and PV regressor were tested with mSLOMOCO to investigate the efficacy of each regressor(s). To compare these results with the performance of VOLMOCO on intervolume/intravolume motion-injected SIMPACE data, 3 pipelines of VOLMOCO, VOLMOCO + 6 Vol-mopa, and VOLMOCO + 6 Vol-mopa + PV were also tested. To compare these results with those obtained using oSLOMOCO, we used the following SLOMOCO package: slomoco_afni_07_2014, from www.nitrc.org/projects/pestica. VOLMOCO was applied to SIMPACE data before oSLOMOCO. oSLOMOCO output (oSLOMOCO) and 14 voxel-wise motion nuisance regress-out outcomes (oSLOMOCO + 14 vox-reg) were presented.

Statistical analysis

Instantaneous intervolume motion

3D rigid volume motion was calculated and corrected using 3dvolreg in AFNI⁶⁴ based on the initial volume, and the estimated volume motion was compared to the injected motion. To investigate the source of the residual signal after motion correction, 15 time series of SIMPACE data with injected z-directional shifts of 2 and 4 mm at the 6^th volume after motion correction were fitted to the external stimulus at the 6^th ([0,0,0,0,0,1,0,…,0]) and 7^th ([0,0,0,0,0,0,1,0,…,0]) volumes, including the linear detrending term in the fitting model. Because the intervolume motion was injected and corrected at the 6^th volume, this volume represented the residual motion artifact directly related to motion and its correction, whereas the 7^th volume represented the time-delayed residual motion artifact or the spin history effect (if any).

Linearly increasing shift

Rigid volume motion was estimated and compared to the injected motion. To investigate the source of residual motion artifact after idealized motion correction, the SIMPACE datasets with injected intervolume motion were motion corrected, and the time series of the motion-corrected SIMPACE datasets or the residual motion artifact after motion correction was visually assessed.

In vivo intervolume motion-injected SIMPACE data

Estimated versus injected motion parameters

The estimated intervolume motion parameters were compared with the injected intervolume motion parameters. A Pearson correlation coefficient (CC) was calculated for the injected and estimated intervolume motion on the 10 in vivo intervolume motion-injected SIMPACE datasets.

We observed a noticeable discrepancy between estimated and injected motion on the 10 in vivo intervolume motion-injected SIMPACE datasets but only a marginal discrepancy between estimated and injected motion in instantaneous or slice motion-injected SIMPACE datasets (see Results and Supplementary Fig 3). We suspected that this was caused by a rounding error in the affine matrix calculation, with a corresponding altered gradient amplitude of imaging axis during the SIMPACE acquisition; when a small amount of motion was injected, the accumulated rounding error over 146 volumes would be noticeable. To test this hypothesis, the derivatives of injected x-, y-, and z-shifted motion were calculated for the 10 in vivo intervolume motion-injected SIMPACE datasets and the positive values of the derivatives were counted. Since the imaging axes are altered based on the previous coordinates during SIMPACE acquisition, the counted number of the positive derivatives in the injected motion would be proportional to the accumulated rounding error, if any. The number of the positive derivatives in x-, y-, and z-shifted motion were plotted against the discrepancy between the estimated and injected shift motion at the last volume on the 10 in vivo intervolume motion-injected SIMPACE datasets.

SD comparison in the GM after application of the motion nuisance models

F-tests for PV and 6 Deriv were conducted using a nested model of 6 Vol-mopa nuisance regressors. The F values for PV and 6 Deriv regressors were calculated by the variance ratio of model 2 (VOLMOCO + 6 Vol-mopa + PV or 6 Deriv) to model 1 (VOLMOCO + 6 Vol-mopa).

We hypothesized that the PV regressor would predict the residual motion artifact in the dataset with large motion drift or out-of-voxel dimension. To test this, we calculated the time series of volume displacement from the first reference volume using the "-maxdisp" option in 3dvolreg. Averaged and maximum time series of displacement in each SIMPACE dataset were calculated to represent the degree of out-of-voxel dimension motion.

The efficacy of the PV regressor was calculated from the reduction in temporal SD in the GM mask after VOLMOCO + 6 Vol-mopa with and without the voxel-wise PV nuisance regressor:

$$1 - \frac{SDafterVOLMOCO + PV}{SDafterVOLMOCO}$$

In vivo intervolume/intravolume motion-injected SIMPACE data

A Pearson CC was calculated for the injected and estimated intervolume/intravolume motion on each of 10 in vivo intervolume/intravolume motion-injected SIMPACE datasets with 1× and 2× amplified intravolume motion. The intervolume/intravolume motion was estimated from VOLMOCO and mSLOMOCO. Since VOLMOCO does not provide intravolume motion estimation, the estimated intravolume motion was set to zero for VOLMOCO estimated motion.

The averaged SD in the GM was calculated after VOLMOCO, oSLOMOCO, and mSLOMOCO. The motion nuisance regressor models were applied to each motion correction as follows: VOLMOCO + 6 Vol-mopa + PV, oSLOMOCO + 14 Vox-reg, and mSLOMOCO + 12 Vol-/Sli-mopa + PV. A paired Student t-test was conducted using calculated SD values in the 10 SIMPACE datasets.

RESULTS

Instantaneous intervolume shift motion-injected SIMPACE data

For the SIMPACE data with injected z-shifts of 2 and 4 mm at the 6^th volume, the estimated shifts were 1.975 and 4.002 mm, respectively. Most of the images demonstrated substantial residual artifact only in the 6^th volume with the half voxel shift and corrected SIMPACE data. No residual artifact was seen in the 7^th volume with either the half or single voxel shift, indicating no substantial time-delayed residual motion artifact after motion correction (Fig. 6). Single voxel shifted and corrected SIMPACE data also did not demonstrate substantial residual motion artifact, indicating that the residual motion artifact was mainly related to the PV averaging effect of the surrounding voxels due to alignment.

Figure 6.Student t-test result for the residual signals of SIMPACE data to the instant stimulus at the 6^th and 7^th volumes. Half (2 mm) and single (4 mm) voxel shifts of intervolume motion are injected in the z-direction at the 6^th volume and corrected. The color map is scaled from +20 (red) to –20 (blue), and thresholded at |t| >3.81 (degrees of freedom = 13, two-sided P < .001).

For SIMPACE data into which 4° of rotation had been injected, the estimated rotation values were 4.119°, 4.002°, and 4.121° in the x-, y-, and z-axes, respectively. Rotation-injected and corrected SIMPACE data demonstrated substantial residual motion in the 6^th volume but not in the 7^th volume (data not shown).

Linearly increasing shift intervolume motion-injected SIMPACE data

After linearly increasing intervolume motion correction, 2 periodic residual artifact patterns were observed in the x-directional voxel dimension of 2 mm (Fig. 7D). The residual motion artifact pattern was also found to be dependent on contrast in the adjacent voxels, leading to different patterns of residual motion artifact at each voxel (Fig. 7E1-8).

Figure 7.(A) SIMPACE data with linearly increasing 4-mm shift of intervolume motion in the x-direction; the green voxel highlights where the sample time series in the other panels are taken from. (B) Injected motion amplitude. Percentage time series signal at the green voxel location from (A), demonstrating results from (C) before and (D) after volumetric rigid motion correction. (E1-8) Percentage time series at 8 neighboring voxels of the green voxel after motion correction.

In vivo intervolume motion-injected SIMPACE data

Estimated versus injected motion parameters

With the exception of x-rotation parameters, Pearson CC values for the time series of injected and estimated motion parameters were typically higher than 0.96 (Table 1); the value for the x-rotation was 0.879. The relatively low correlation for x-rotation was a consequence of the small dynamic range of the injected motion parameters in SIMPACE datasets 7, 8, 9, and 10 (Supplementary Fig. 3 and Table 1).

In general, the offset between injected and estimated intervolume motion patterns was small (Supplementary Fig. 3). When VOLMOCO was applied, SoS values in the GM were 15% lower with estimated parameters than with injected parameters; when VOLMOCO + 6 Vol-mopa was used, SoS values were 16% lower with estimated parameters than with injected parameters (Supplementary Table 1). Supplementary Figure 4 presents the linear correlation between the number of positive derivatives in the injected shift motion and the discrepancy between the estimated and injected shift motion at the last volume (R² = 0.46, P < 10^–4). This finding indicates that the discrepancy between injected and applied motion was accumulated proportional to the frequency of positive or negative direction of the injected motion, supporting the hypothesis that applied motion in SIMPACE is systemically smaller than injected motion. For this reason, we used the estimated rigid volume motion parameters rather than the injected motion parameters for the following analyses.

SD comparison in the GM

The averaged SD value in the GM without motion injection or correction during SIMPACE acquisition was 4.92, which is a baseline SD for thermal noise during SIMPACE acquisition (Fig. 8). The averaged SD values in the GM across the 10 intervolume motion-injected SIMPACE datasets after intervolume motion correction were: 11.99 ± 3.94 with VOLMOCO only, 9.31 ± 2.42 with VOLMOCO + PV, 7.55 ± 1.90 with VOLMOCO + 6 Vol-mopa, 6.82 ± 1.40 with VOLMOCO + 6 Vol-mopa + PV, and 6.60 ± 1.32 with VOLMOCO + 6 Vol-mopa + 6 Deriv. Although VOLMOCO + 6 Vol-mopa + 6 Deriv provided the smallest averaged SD, the marginal difference in SD between VOLMOCO + 6 Vol-mopa + PV and VOLMOCO + 6 Vol-mopa + 6 Deriv was offset by a gain of 5 DOF. F-test results for VOLMOCO + 6 Vol-mopa + PV and VOLMOCO + 6 Vol-mopa + 6 Deriv based on VOLMOCO + 6 Vol-mopa were F(1,135) = 38.59 and F(6,132) = 8.90, respectively (all P < 10^–4).

Figure 8.Averaged temporal standard deviation (SD) in the gray matter (GM) after application of different motion nuisance regressors in 10 intervolume motion-injected and corrected SIMPACE datasets. The dashed line at 4.92 indicates average SD value in the GM when no motion was injected into the SIMPACE data. VOLMOCO: rigid volume motion correction; Vol-mopa: volume motion parameter nuisance regressor; PV: partial volume nuisance regressor; Deriv: Derivatives of Vol-mopa.

Figure 9 shows the linear correlation between the reduction in averaged SD in the GM after application of the PV nuisance regressor based on linearly detrended residual signals and the maximum displacement output of 3dvolreg across the 10 intervolume motion-injected SIMPACE datasets. A larger SD reduction was observed with larger maximum and mean values of time series maximum volume displacement for brain voxels. The largest reduction was seen in the SIMPACE 6 dataset (Supplementary Fig. 2), which demonstrated a drift greater than the 2-mm shift in the z-direction. This finding agreed with the results seen for linearly increasing intervolume motion-injected SIMPACE data. PV was found to predict the repeated residual artifact pattern out of voxel dimension (Fig. 7D). For this reason, the final mSLOMOCO pipeline employed the PV regressor rather than derivative terms as nuisance regressors.

Figure 9.Plot of the reduction in temporal standard deviation (SD) in the gray matter (GM) with voxel-wise partial volume (PV) nuisance regressor based on the linear detrending and maximum displacement (mm) output of 3dvolreg in 10 intervolume motion-injected SIMPACE datasets. Black dots represent maximum values and white dots represent mean values of time series maximum displacement with linear fit lines for each.

In vivo intervolume/intravolume motion-injected SIMPACE data

Accuracy of the estimated motion parameters

Compared to VOLMOCO, mSLOMOCO demonstrated similar or higher Pearson CC values for estimated versus injected motion parameters, especially for the y-shift motion parameter (Fig. 10). The accuracy values for x-rotation and y-shift estimation were relatively lower than those for other motion parameters with both VOLMOCO and mSLOMOCO. When intravolume motion was increased, the accuracy of the motion parameter estimation decreased; this degraded accuracy was mainly observed when VOLMOCO was used (Fig. 10). The representative CC values for y-shift motion were as follows: 0.79 ± 0.09 versus 0.50 ± 0.16 from VOLMOCO, and 0.87 ± 0.10 versus 0.79 ± 0.14 from mSLOMOCO with 1× versus 2× amplified intravolume motion, respectively. An example of estimated motion using VOLMOCO and mSLOMOCO is shown in Figure 11.

Figure 10.Comparison of motion estimation using rigid volume motion correction method (VOLMOCO, blue) and modified slice-oriented motion correction method (mSLOMOCO, red). Pearson correlation coefficient (CC) was calculated for the time series of injected and estimated motion parameters. For VOLMOCO, the estimated intervolume motion parameters were temporally concatenated for intervolume/intravolume motion during volume acquisition. rot: rotation.

Figure 11.(A-C) Examples of injected and estimated in-plane motion in intervolume/intravolume motion-injected SIMPACE data. Black, blue, and red lines indicate the injected intervolume/intravolume motion, estimated intervolume motion using 3dvolreg, and estimated intervolume/intravolume motion using modified slice-oriented motion correction (mSLOMOCO), respectively. (D) Zoomed-in portion of the image shown in (C).

SD comparison in the GM after motion correction

We found that a) oSLOMOCO and mSLOMOCO perform similarly in reducing the residual motion artifact, and b) mSLOMOCO pipeline with 12 Vol-/Sli-mopa and PV regressors outperforms oSLOMOCO pipeline with 14 voxel-wise regressors, and VOLMOCO pipeline with 6 Vol-mopa and PV regressors. Paired Student t-tests for SD values in the GM demonstrated significant differences between VOLMOCO and oSLOMOCO (P < 10^–3) and between VOLMOCO and mSLOMOCO (P < 10^–3), but not between oSLOMOCO and mSLOMOCO (P = 0.49). mSLOMOCO with 12 Vol-/Sli-mopa and voxel-wise PV regressors demonstrated a significantly lower SD in the GM than VOLMOCO + 6 Vol-mopa + PV and oSLOMOCO + 14 Vox-reg with both 1× and 2× amplified intravolume motion injection (paired t-test, P < 10^–3).

The averaged SD in the GM was observed to be higher when the intravolume motion was amplified, and the difference in averaged SD between 1× and 2× amplified intravolume motion was lower with the mSLOMOCO pipeline (5.47± 0.86 versus 5.76 ± 1.09) than with the VOLMOCO (7.69 ± 1.77 versus 10.47 ± 3.00) or oSLOMOCO pipeline (7.57 ± 2.04 versus 8.33 ± 2.62). Figure 12 and Supplementary Table 3 show the SD values in the GM with the different motion correction methods and their motion nuisance regressor models.

Figure 12.Averaged temporal standard deviation (SD) comparison for each of ten in vivo intervolume/intravolume motion-injected SIMPACE datasets with 1× (white bars) and 2× (gray bars) amplified intravolume motion using VOLMOCO, oSLOMOCO, or mSLOMOCO and their motion nuisance regressor models. Dashed line indicates 4.92 of SD from SIMPACE data without motion injection. VOLMOCO: rigid volume motion correction; Vol-mopa: volume motion parameter nuisance regressor; PV: partial volume nuisance regressor; SLOMOCO: slice oriented motion correction; oSLOMOCO: original SLOMOCO; mSLOMOCO: modified SLOMOCO; Vox-reg: voxel-wise nuisance regressor; Sli-mopa: slice-wise rigid volume motion parameter nuisance regressor.

Figure 13.Example of the standard deviation (SD) map for intervolume/intravolume motion-injected SIMPACE data after application of (A) VOLMOCO + 6 Vol-mopa + PV and (B) SLOMOCO + 12 Vol-/Sli-mopa + PV, scaled up to 14 and 10, respectively. (C) SD maps from 10 SIMPACE datasets with intervolume/intravolume motion injection after SLOMOCO + 12 Vol-/Sli-mopa + PV.

DISCUSSION

In this study, we generated motion-corrupted MR data by altering the imaging plane coordinates before each volume and slice acquisition from an ex vivo brain phantom. In this manner, the intervolume/intravolume motion-injected SIMPACE data became a gold standard motion-simulated dataset to emulate motion effects during MR data acquisition. Testing with these SIMPACE datasets, we validated using the described PV motion nuisance regressor and mSLOMOCO method as a means to reduce the effects of subject motion in fMRI data analysis. In the tests performed here, the proposed mSLOMOCO pipeline (including a PV regressor) effectively reduced the residual motion artifact and outperformed each of the standard VOLMOCO and oSLOMOCO pipelines.

In vivo motion characteristics

We observed slow drift and spontaneous estimated motion in the x-shift/z-shift and y-rotation/z-rotation and periodic estimated motion in the y-shift and x-rotation (Supplementary Fig. 3). In real data, y-shift (= PE directional) motion is typically related to artificial motion caused by B0 field fluctuation due to breathing.¹⁴ Periodic x-rotation motion may also be related to actual breathing-related motion and/or to estimated motion due to PE directional artificial head motion. Empirically, this y-directional motion continuously fluctuates from 0.2 to 0.4 Hz, indicating that intravolume motion always exists across slices, even in fast fMRI acquisitions, when 2D EPI is used. For this reason, intravolume y-shift and x-rotation motion patterns are the most difficult to estimate and are believed to be one of the main sources of bias in intervolume motion assumptions; this is reflected in the relatively lower CC between injected and estimated y-shift motion with VOLMOCO than with mSLOMOCO (Fig 10). Figure 11D shows the volume-wise estimated y-shift motion with VOLMOCO (blue line) versus the continuous y-shift motion that occurs during intravolume acquisition (black line).

Injected and estimated motion on in vivo SIMPACE datasets

Although there was only marginal difference between estimated and injected motion on instantaneous motion-injected SIMPACE datasets (ie, 1.975 and 4.002 mm for estimated motion vs 2 and 4 mm for z-directional injected motion), the discrepancy between estimated and injected motion on the in vivo intervolume motion-injected SIMPACE datasets was noticeable (Fig. 11 and Supplementary Fig. 3). Our data indicated that the discrepancy between estimated and injected motion was increased at the last volume when the injected motion at each volume was applied in a positively or negatively skewed manner over total volume acquisition (Supplementary Fig. 4). As the differential motion was injected based on the previous volume in SIMPACE acquisitions, a rounding error in the affine matrix and corresponding gradient amplitude in the re-defined imaging axis could accumulate when motion was injected in a single direction and could be cancelled out when motion was positively and negatively counterbalanced. The specific injected motion information and SIMPACE data will be publicly available (to be announced in https://github.com/wanyongshinccf/SLOMOCO).

Voxel-wise PV regressor

The PV regressor generated for this study was effective in reducing the residual motion artifact, especially in cases with large motion drift. We generated this regressor by aligning the static images to each time point of motion⁴⁵ and then re-aligning them back to the reference. This process included 2× amplified resampling artifact in the PV regressor, which may have introduced bias in the ability to predict residual motion artifact, especially in regions surrounded by different types of tissue contrast (eg, the subcortical region near ventricles). Higher resolution EPI images than those used in these fMRI datasets could be useful in generating these data and minimizing this bias; however, it would be challenging to generate similar tissue contrasts and levels of EPI distortion/artifact with these higher resolution EPI images in an fMRI dataset. Despite these limitations, we recommend including a PV regressor as a motion nuisance regressor even in cases of rigid intervolume motion assumption or VOLMOCO, as adding the PV nuisance regressor reduces the DOF by only one, while appearing to notably reduce subject motion effects.

We note that the effectiveness of the PV regressor was also dependent on the interpolation algorithm used. We tested the different interpolation options in 3dvolreg (-linear, -cubic, -quintic, and -wsinc5) to generate PV regressors and found that the -cubic option generated the smallest SD in GM when the PV nuisance regressor was added to VOLMOCO + 6 motion parameters (Supplementary Table 3).

VOLMOCO vs mSLOMOCO

With intervolume/intravolume motion-injected SIMPACE data, mSLOMOCO outperformed VOLMOCO. When 1× and 2× amplified intravolume injection was used, the average SD values in the GM were 29% and 45% smaller with the mSLOMOCO + 12 Vol-/Sli-mopa + PV pipeline than with the VOLMOCO + 6 Vol-mopa + PV pipeline, respectively. Intravolume motion and its amplitude did degrade intervolume motion correction with both VOLMOCO and mSLOMOCO, indicating that spontaneous motion could be one of the main residual artifact sources even after motion correction, as shown in previous studies.^36,38,65 However, considering that the SD value in the GM in ex vivo brain EPI scans without motion injection was 4.92, the mSLOMOCO pipeline led to increases of 11% and 17% for SD in the GM across the 10 SIMPACE datasets with 1× and 2× amplified intravolume motion injection, respectively, after intervolume/intravolume motion correction, showing an additional 6% increase due to the imperfect motion correction on 2 × amplified intravolume motion. It should be noted that VOLMOCO produces a 57% increase in SD value when intravolume motion is amplified by two. This could be explained by the fact that the amplitude of intravolume motion, induced by spontaneous head motion, reduces the efficacy of VOLMOCO. These results suggest that mSLOMOCO minimizes bias from residual motion artifact that is caused by intravolume motion.

oSLOMOCO vs mSLOMOCO

We found that oSLOMOCO and mSLOMOCO performed similarly in reducing residual motion artifact (Fig. 12 and Supplementary Table 3). However, the different motion nuisance regressor models of oSLOMOCO and mSLOMOCO produced significant differences with both 1× and 2× amplified intravolume injections (all P < 10^–3). It should be noted that mSLOMOCO + 12 Vol/Sli-mopa + PV produced the smallest variation of SD in the 10 different motion patterns of SIMPACE datasets with 1× and 2× amplified intravolume motion, indicating a robust performance in reducing residual artifacts with minimal dependency on the different motion patterns and amplitudes. The improvement seen with the mSLOMOCO pipeline is mainly induced by the new motion nuisance regressor model including the voxel-wise PV regressor.

Spatial variation of residual motion artifact after motion correction

Figure 13A and 13B show an example of SD maps for SIMPACE 6 dataset after application of the VOLMOCO and mSLOMOCO pipelines, displaying the spatial variation in SD. The spatial variation of the residual motion artifact was also observed even without the motion nuisance regressor models (result not shown). Figure 13C shows the same location of SD maps for 10 SIMPACE datasets after the application of mSLOMOCO with 12 Vol-/Sli-mopa + PV, demonstrating the different spatial patterns of the residual motion artifact according to 10 injected motion patterns. These findings suggest that it is important to characterize head motion in fMRI data to effectively predict the degree of residual motion artifact and to consequently censor bad-quality fMRI data or the corresponding volume (discussed in the next section).

Limitations and future directions

In this study, we generated and tested SIMPACE data using fixed MR scan parameters (eg, TR of 2 s and voxel size of 2 × 2 × 4 mm³ for the in vivo dataset). The effects of subject motion depend on tissue contrast, whether VOLMOCO or SLOMOCO is applied. It is expected that using a short TR with in-plane acceleration or simultaneous multislice (SMS) excitation⁶⁶ will lead to relatively small intravolume motion within a TR, improving VOLMOCO performance. However, the relatively short TR and small flip angle with SMS acceleration reduce the signal-to-noise ratio (SNR) and contrast-to noise-ratio (CNR) in brain tissues, in addition to the slice cross-talk and acceleration penalty.⁶⁷ The effect of reduced SNR and CNR on the efficacy of the motion estimation with VOLMOCO and SLOMOCO was not tested in this study. Additional studies employing variations in MRI parameters are therefore needed to further generalize the results observed here.

In this study, we did not present the head motion index to be utilized for volume censoring or the exclusion criteria for fMRI datasets due to excessive head motion. Motion-corrupted volumes are often censored and motion-corrupted fMRI data are commonly excluded in fMRI data analysis.^36,65,68 However, it remains a challenge determine which motion is bad motion for an fMRI study and how bad it is. SIMPACE data with a priori motion parameters and without physiologic noise could aid in the determination of which volumes should be censored and what could be the exclusion criteria for fMRI datasets. This was beyond the scope of the current study but will be assessed in future research. Validation of the proposed SLOMOCO technique with a motion nuisance regressor will be conducted in a large number of in vivo fMRI studies in the future.

CONCLUSION

This study assessed the use of intervolume and intervolume/intravolume motion-injected SIMPACE data from an ex vivo brain phantom and investigated residual motion artifact signal on intervolume motion-corrected SIMPACE data. Based on our findings, we propose the use of a voxel-wise PV nuisance regressor to remove residual motion artifact. We found that the proposed modified SLOMOCO pipeline with 12 Vol-/Sli-mopa and PV regressors reduced residual motion artifact in intervolume/intravolume motion-injected SIMPACE data more effectively than standard VOLMOCO and oSLOMOCO pipelines. Moreover, the SLOMOCO pipeline reduced residual motion artifact variations across the different patterns of head motion. The modified SLOMOCO script can be found at https://github.com/wanyongshinccf/SLOMOCO and the tested intervolume and intervolume/intravolume motion-injected SIMPACE data will be made publicly available (to be linked via the GitHub address, above).

Data/Code Availability

SLOMOCO software is available now in https://github.com/wanyongshinccf/SLOMOCO

SIMPACE data will be shared in a public (will be announced in the link above).

Author Contributions

Wanyong Shin: Conceptualization, Methodology, Investigation, Writing - Original Draft, formal analysis, Software development.

Paul Taylor: Software development, Reviewing and editing the manuscript.

Mark J. Lowe: Writing, Reviewing and Editing. Supervision. Project administration.

Acknowledgements

Authors appreciate M.D Sanghoon Kim, Ph.Ds Jagjit Singh and Jacquline Chen for ex-vivo brain phantom. Authors also appreciate Megan Griffiths for editing the manuscript. The authors have no funding sources to declare.

Conflicts of Interest

The authors declare no competing interests.