Towards motion insensitive EEG-fMRI: Correcting motion-induced voltages and gradient artefact instability in EEG using an fMRI prospective motion correction (PMC) system

The simultaneous acquisition of electroencephalography and functional magnetic resonance imaging (EEG-fMRI) is a multimodal technique extensively applied for mapping the human brain. However, the quality of EEG data obtained within the MRI environment is strongly affected by subject motion due to the induction of voltages in addition to artefacts caused by the scanning gradients and the heartbeat. This has limited its application in populations such as paediatric patients or to study epileptic seizure onset. Recent work has used a Moiré-phase grating and a MR-compatible camera to prospectively update image acquisition and improve fMRI quality (prospective motion correction: PMC). In this study, we use this technology to retrospectively reduce the spurious voltages induced by motion in the EEG data acquired inside the MRI scanner, with and without fMRI acquisitions. This was achieved by modelling induced voltages from the tracking system motion parameters; position and angles, their first derivative (velocities) and the velocity squared. This model was used to remove the voltages related to the detected motion via a linear regression. Since EEG quality during fMRI relies on a temporally stable gradient artefact (GA) template (calculated from averaging EEG epochs matched to scan volume or slice acquisition), this was evaluated in sessions both with and without motion contamination, and with and without PMC. We demonstrate that our approach is capable of significantly reducing motion-related artefact with a magnitude of up to 10mm of translation, 6° of rotation and velocities of 50mm/s, while preserving physiological information. We also demonstrate that the EEG-GA variance is not increased by the gradient direction changes associated with PMC. Provided a scan slice-based GA template is used (rather than a scan volume GA template) we demonstrate that EEG variance during motion can be supressed towards levels found when subjects are still. In summary, we show that PMC can be used to dramatically improve EEG quality during large amplitude movements, while benefiting from previously reported improvements in fMRI quality, and does not affect EEG data quality in the absence of large amplitude movements.


Introduction
Simultaneous electroencephalography and functional magnetic resonance imaging (EEG-fMRI) is a multimodal technique that aims to benefit from the temporal resolution of brain activity from the electroencephalography (EEG) with the spatial resolution of fMRI. EEG-fMRI has been widely applied to study healthy (Mulert et al., 2004;Laufs et al., 2003) and pathologic brain activity, such as for studying patients with epilepsy (Lemieux et al., 2001;Laufs et al., 2008;Centeno and Carmichael, 2014;Pittau et al., 2012;LeVan et al., 2010). More recently EEG-fMRI has been proven to be capable of mapping BOLD signal changes associated with epileptic seizures (Chaudhary et al., 2012a) using events detected in EEG. However this endeavour can be severely hindered by subject motion.
Subject motion during the acquisition of fMRI can lead to severe data quality degradation. The impact of motion on fMRI is well documented (Hajnal et al., 1994;Satterthwaite et al., 2012) and it causes large amplitude signal changes across consecutive fMRI volumes increasing temporal variance and increasing type 1and type 2 errors. A number of techniques ranging from purely data-based (post-processing) (Friston et al., 1996;NeuroImage 138 (2016) 13-27 Lemieux et al., 2007;Power et al., 2014;Tierney et al., 2015) to methods based on precise and independent measurement motion have been proposed (Eviatar et al., 1997;Eviatar et al., 1999). Recently, a prospective fMRI motion-related signal reduction system based on a camera-tracker system (PMC) and MR sequence acquisition update has been implemented and commercialised. Briefly, an MRI compatible camera is used to image a Moiré-phase tracker (MPT) attached to the head. The tracking information is converted into the head's position (three translations and three rotations). This information is then used to update the radio frequency pulses and the gradients, applied in the imaging process in real time with promising results Todd et al., 2015). EEG quality degradation due to motion is mostly the result of electromagnetic induction. Faraday's law states that magnetic flux changes through a conductor loop's area induce corresponding voltage fluctuations; in EEG circuits these voltages are superimposed onto the braingenerated signals, making the detection of activities of interest more difficult or even impossible. This is even more problematic for the application of EEG-fMRI in patients for whom the motion might be unavoidable, such as children or to study patients with epilepsy during seizures. Currently, most artefact correction methods are based on post-hoc EEG data processing techniques that rely on the identification of artefact waveforms and their subtraction with the aim of obtaining motion artefact-free signals. These processes do not take into account the measured motion (Lemieux et al., 2007;Chaudhary et al., 2012b).
Systems to detect motion induced voltages based on dedicated sensors, such as piezoelectric devices (Bonmassar et al., 2002) and carbon loop wires (Masterton et al., 2007), have shown promise. While these are able to deal with small amplitude movements (~1 mm and 1°), neither performed well in removing voltages induced by large motion events (Masterton et al., 2007). Furthermore, in the case of the wire loops, there is the requirement for additional technology related to the data acquisition and this can pose safety risks in some circumstances (Abbott et al., 2014). Finally, the additional requirement for fMRI data correction during subject motion is not addressed by these systems thereby limiting their impact. A PMC camera system has only previously been used for correcting the ballistocardiogram (BCG) artefact (LeVan et al., 2013) found in EEG data acquired inside MRI scanners. While promising this previous work did not address large scale movements or the effects of applying PMC itself on EEG quality during fMRI data acquisition which may suggest that correction of EEG and fMRI is problematic using this approach.
In this study, we aimed to focus on improving EEG data quality while suppressing fMRI motion artefacts using a commercially available PMCcamera system. We derived a model of voltage changes induced in the EEG from motion using the accurate measurement of head's motion recorded by the PMC system and used it to attenuate these artefactual signals. We tested this approach with and without large amplitude movements by modelling and removing motion induced voltages and assessed the EEG quality. Additionally, we determine the impact of PMC on the gradient artefact template temporal stability (i.e. the variance between EEG epochs for each fMRI volume or slice) which contains both motion and gradient artefact instability. We also verify our experimental findings by determining the effect of PMC updating the magnetic field gradients on the EEG gradient artefact (GA) in the presence of motion (i.e. the variability of the voltages induced exclusively by the magnetic field gradient switching with or without PMC).

Data and methods
We acquired simultaneous EEG-fMRI data in three healthy subjects (two male, both 25 years old, and a female, 23 years old).

Task
During all recordings subjects were instructed to open and close their eyes every 60 (for sessions outside the scanner) or 30 s (for sessions inside the scanner). Two different movement tasks were performed: In the 'keeping still' sessions subjects were instructed to keep as still as possible. In 'motion' sessions subjects were instructed to move performing repeated repetitions of the following movements (see Fig. 1 and the Video 1 in Supplementary material): shaking their head side to side (right-left), nodding their head (back and forth) and rotating their head followed by a short period (~5 s) without movement. The subjects were instructed to start the first block with their eyes open allowing them to calibrate the motion's amplitude (via visual feedback of the marker position). In the second block the movements were repeated with eyes closed. Three repetitions of each block were made in each session.

EEG data acquisition
During all recordings subjects were instructed to alternate between eyes opened and eyes closed via a verbal queue. This was done in order to evaluate the practical contribution of EEG-motion artefact correction to measure physiological signal changes from the brain (the alpha rhythm in this case). Fig. 1 presents a schematic diagram summarising the acquisition setup.

MRI acquisition
All images were acquired at Great Ormond Street Hospital, London, United Kingdom, using a 1.5T Avanto (Siemens, Erlangen, Germany). The functional images were acquired using a gradient-echo echoplanar imaging (GE-EPI) sequence with the following parameters: 30 ascending slices (3 × 3 × 3 mm 3 ) covering the whole brain with a slice gap of 0.5 mm, slice/volume TR = 74/2220 ms, TE = 30 ms, flip angle = 75°. Each fMRI run was composed of 100 volumes, obtained in 222 s. The posterior half of a 32 channel head coil was used for signal reception to allow for video recording of the subject motion and visual feedback of the motion to be provided to the subject (see Section 2.6).

EEG recording
The EEG data were recorded using a MR-compatible Amplifier (Brain Products, Gilching, Germany). An MRI compatible cap (Easycap, Herrsching, Germany) with 31 EEG-electrodes placed as the 10-20 international system, and an ECG channel, all made of Ag/AgCl with internal safety resistors and referenced to FCz was used. The EEG data were acquired with a sampling rate of 5 kHz and band-pass filtered from 0.1 to 250 Hz. The EEG acquisition was synchronised with the MRscanner by the SyncBox device (Brain Products, Gilching, Germany), allowing the EEG acquisition to be time-locked with the fMRI acquisition gradient switching-related artefact, facilitating its correction (Mandelkow et al., 2006). Before analysing, all the EEG data were down sampled to 500 Hz.

Motion measurement and fMRI-prospective motion correction
A MR-compatible camera (Metria Innovation Inc., Milwaukee, USA) was used for tracking a Moiré Phase Tracking (MPT) marker  attached to a 'bite bar' specifically developed for each subject based on a dental retainer.
To make the bite bar a dental stone cast was produced from a dental alginate impression taken of each subject's maxillary dentition. A 'duallaminate' (DB Orthodontics, West Yorkshire, England), thermoplastic, 3 mm thick bite guard was then thermoformed to the maxillary dentition using a Biostar machine (Scheu-Dental, Iserlohn, Germany). The 'soft' side of the 'dual-laminate' ensured a snug fit around the teeth whilst facilitating easy insertion and withdrawal of the guard, whilst the 'hard' side ensured that the bite bar didn't move when the subject bit down. A 6 cm × 2 cm × 3 mm strip of 'dual-laminate' was then formed into a bar with a 90 degree angle using localised heat. The bar also incorporated housing for the MPT marker. Finally, the bar was joined to the incisor region of the bite guard using medical grade cyanoacrylate glue.
The MPT camera was used for recording the marker motion with six degrees of freedom, three translations (x, y and z axis, defined as rightleft, posterior-anterior and feet-head, respectively) and three rotations (about x, y and z axis, respectively) with a sampling rate of 85 Hz and recorded on a computer located outside the scanner room. The same computer was connected to the scanner to update the Radio Frequency and Gradients pulses used for MRI signal acquisition before every fMRI slice is acquired Todd et al., 2015). Motion tracking was used in the scanner to record the motion parameters throughout the experiments. Real-time updating of the MRI gradients was switched either on or off in the scanner's software in different sessions to test its effect (see Section 2.8).
We analysed the motion characteristics for each subject by calculating the root-mean-square (RMS) for each of the six velocities derived from motion measurements (three translations and three rotations) and detecting the minimum and maximum amplitude on the axis with the largest RMS. We also calculated the fast Fourier transform on the motion data from the same axis, to allow us to check the motion's spectral characteristics. For simplicity we reported the frequency of the highest motion-related peak (e.g. the highest peak in the spectra at a non-zero frequency).

Mirror, video camera and screen
We were interested in large amplitude movements in the maximum range that can be reliably tracked with our motion tracking system within the range from − 10 to 10 mm and − 6 to 6°. Subjects were Fig. 1. Schematic overview diagram of the motion corrected (green arrows) and uncorrected (red arrows) EEG-fMRI acquisition. The MPT camera is mounted in the scanner above the subjects head (A) and is used to record head position (illustrated by a head model, B top row) via a MPT marker (B, second row). As the subject moves (as per the task in C) the marker movement is recorded by a camera (B, second row). fMRI will normally be affected by this motion (B, third row) however the marker movement is converted into the subjects head position which can be fed back to the scanner (A) in real time to update the RF and gradients pulses, stabilising the fMRI images (B, bottom row). The EEG data are also affected by the motion (D, top row). However, it is possible to remove the spurious voltages induced by this motion using a model derived from the same position information (3 translational and 3 rotational velocities are illustrated in the middle row of figure D), resulting in the corrected EEG (D, bottom row). provided with a visual display of the marker position in real time via a screen mounted outside the bore and a mirror to allow them to monitor their movements during the task. We also recorded video of the subject with the EEG in BrainVision Recorder (Brain Products, Gilching, Germany) via a second MRI compatible camera (NNL eye-tracker camera NordicNeuroLab, Bergen, Norway) in order to visually monitor movement and aid the identification of events in the EEG.
Moiré phase tracking-retrospective EEG motion artefact suppression (MPT-REEGMAS) The EEG time series data at any given electrode acquired inside the scanner in the presence of motion can be described as follows: where the voltage acquired by each electrode (EEG raw ) is the linear sum of the voltage generated by the neurons and physiological artefacts, such as eye-blinks and muscular activities (EEG physiological ) and a motion-related voltage(EEG motion ). The latter can be related to changes in magnetic flux by the Faraday's Law of induction: where dϕ B dt is the temporal derivative of the magnetic flux through the electrode + wire + head circuit area.
Here we assumed that dϕ B dt is purely a function of head motion relative to the laboratory (here scanner) coordinates. Motion can be due to the movement of the electrodes-head system and/or second order effects, such as those related to head acceleration and cap inertia. In both cases, we have a voltage induction in the electrodes proportional to the head motion in a fixed frame.
The following steps were all implemented in Matlab 2013a (The MathWorks, Natick, USA). The MPT-camera recorded motion data were low-pass filtered at 11 Hz using a zero-phase Butterworth filter of order eight resulting in smoother data in order to limit the variance of the derivatives. The 11 Hz cut-off was chosen after checking the motion parameters spectrogram and verifying that the main signal related to subject motion was below this frequency. We note that this rate is also the typical rate of 2D-EPI data acquisition, where there is rarely motion fast enough to corrupt within slice acquisition (i.e. b100 ms timescale). The 6 velocities were calculated as the temporal derivatives of the MPT-derived motion parameters. We also calculated the squared velocities to account for second order effects that might be present due to differences in movement between the skull and the EEG cap or non-linear movement of wires. The resulting 18 parameters (6 positions, 6 velocities, 6 velocities squared) were interpolated to 500 Hz to match the EEG sampling rate (after down sampling). We used a general linear model to estimate EEG(t) motion , for time points from t1 to tT, where T is the final time point based on the 18 motion parameters described above: (3) where M(t) p is the matrix of motion parameters at time t and β pi are the factors to be estimated. The Eq. 3 can be written as: Here, ε(t) represents the residual of the EEG temporal series with variance related to our motion model removed (i.e. EEG(t) physiological ).
β can be estimated in order to minimise the residual ε(t) Finally, the motion artefact corrected EEG(t) physiological can be recovered for each electrode at each time point by the following: We refer to the process of calculating the EEG(t) physiological using Eq. 6 and the tracking information as Retrospective EEG Motion Artefact Suppression (REEGMAS).

Experiments and processing Experiment 1: baseline EEG acquisition
We acquired baseline sessions of EEG data outside the scanner with subjects asked to keep still and to open and close their eyes every 60 s.

Experiment 2: EEG recording inside the scanner without scanning (no fMRI)
Two sessions of EEG data were recorded inside the scanner without MRI scanning: the first while keeping still (EEG 2S ) and the second while moving voluntarily (EEG 2M ).
The EEG data was down-sampled to 500 Hz and imported into Matlab. Two processed versions of each EEG dataset were obtained for further analysis: EEG 2S and EEG 2M without REEGMAS correction but with BCG artefact correction by Average Artefact Subtraction (AAS BCG ) (Allen et al., 1998) and EEG 2S-C and EEG 2M-C with REEGMAS correction followed by BCG artefact correction by AAS BCG as implemented in BrainVision Analyzer 2.0 (Brain Products, Gilching, Germany).
Experiment 3: simultaneous EEG-fMRI acquisitions Experiment 3a. EEG and fMRI data were acquired with the PMC system recording motion data but not updating the imaging gradients ('fMRI/ PMC-off'). We have acquired two successive sessions; with the subjects keeping still (EEG 3aS ) and while moving voluntarily (EEG 3aM ), respectively. Experiment 3b. EEG and fMRI data were acquired with the PMC system recording motion and updating the scanning gradients ('fMRI/PMCon'). Two successive sessions were acquired; with the subjects keeping still (EEG 3bS ) and while moving voluntarily (EEG 3bM ), respectively.
The EEG data from both Experiments 3a and 3b were down-sampled to 500 Hz, imported into Matlab and REEGMAS was applied as described in Section 2.7, resulting in motion-corrected EEG (EEG 3aS-C , EEG 3aM-C , EEG 3bs-C and EEG 3bM-C ). In order to correct the GA we have applied the Average Artefact Subtraction (AAS GA ) (Allen et al., 2000). Henceforward we refer to GA correction as AAS GA and to ballistocardiogram (BCG) artefact correction as AAS BCG . In this study we built the templates in two different formats, volume-volume (volume-wise) (Allen et al., 2000) and slice-slice (slice-wise) (Niazy et al., 2005). Firstly, a volume-wise GA template was formed by averaging volumes (n = 15, each 2220 ms, using 33.3 s of data in total) we chose a low number of volumes to reduce motion sensitivity. Secondly, a slice-wise GA template was formed by averaging slice epochs (n = 15, each 74 ms, using 1.1 s of data in total). This is possible because the gradients applied are identical for each slice with the only difference being the RF pulse frequency offset. Following this, AAS BCG was applied to both the REEGMAS corrected EEG (EEG 3aS-C , EEG 3aM-C , EEG 3bS-C , EEG 3bM-C ) and the original EEG without REEGMAS correction (EEG 3aS , EEG 3aM , EEG 3bS , EEG 3bM ). We define the GA template variability as the slicewise or volume-wise epoch-epoch differences that maybe caused by any factor (e.g. motion, scanner instability, GA). The GA is only the voltages resulting from the magnetic field gradients switching themselves.
After the processing described above, the EEG data from all experiments were then down-sampled to 100 Hz and band-pass filtered from 0.5-40 Hz in BrainVision Analyzer 2.0. This is a common filter width for clinical EEG visualisation. All further analysis was performed in Matlab 2013a.

EEG data quality assessment
For the data acquired during the moving sessions of Experiment 2, we evaluated the importance of each parameter to the motion regression model by calculating an F-test score using the function ftest.m.
The quality of the REEGMAS artefact correction was assessed by comparing the EEG obtained in the scanner (Experiments 2 and 3) to the EEG acquired outside the scanner room (Experiment 1) in the following ways: 1. We visually assessed the EEG data comparing the presence of physiologic signals/markers, such as eye-blink artefacts and presence of signal in the alpha-rhythm frequency band (8-12 Hz). To assess the impact of the motion correction across channels, the power topography over the scalp in the alpha rhythm frequency band (8-12 Hz in our study) was visualised. The alpha topography was the power in the frequency band 8-12 Hz averaged over an epoch of 5 s following the first period of eyes closed. We compared EEG data acquired during both 'keeping still' (EEG 2S ), 'motion' uncorrected (EEG 2M ) and corrected (EEG 2M-C ) sessions. 2. The EEG power spectral density (PSD) was calculated by applying the Welch method (pwelch.m function) using Hamming windows of 3 s with 1.5-second overlap during keeping still/motion sessions in eyes closed periods. The PSD was normalised at each frequency by the sampling rate and the number of samples inside each window (dB/Hz). The mean PSD was calculated for the baseline EEG (Experiment 1 during eyes closed periods). Additionally we calculated two standard deviations of this baseline EEG-PSD across time and assumed that points lying outside this range were likely to be due to artefacts. 3. The Root Mean Square Error (RMSE) defined as the difference between the uncorrected (or corrected) EEG-PSD and the baseline EEG mean PSD was calculated for each frequency from 0.5 to 40 Hz with a frequency resolution of 0.33 Hz. Then we calculated the average RMSE over the entire frequency range obtaining the Mean Root Mean Square Error (MRMSE). Finally, we applied a one-sample ttest in order to compare the MRMSE obtained for EEG data before and after REEGMAS correction.
For Experiment 3 data, we assessed GA template (slice-wise) stability for each EEG session by calculating the variance of Root Mean Square EEG amplitude (EEG-RMS) across slices based on the rationale that inter-slice variance is increased by noise and therefore summarises the stability of the GA template (slice-wise). The combined impact of motion and GA induced voltages can then be compared under the different conditions and corrections tested. The EEG-RMS was calculated from the raw/REEGMAS corrected down-sampled EEG, given by the following equation: where E i is the voltage measured by the electrode i and N is the number of electrodes. In this study each fMRI slice was acquired in 74 ms. The EEG was down-sampled to a rate of 500 Hz, resulting in 37 (k) points at which the GA template stability could be estimated (i.e. 37 samples were obtained in the 500 Hz EEG during the 74 ms slice acquisition period).
The GA template (slice-wise) variance was calculated at each of these 37 points across 3000 repetitions of the GA (j corresponding to the number of slices obtained in each session) by the following equation. where, We then compared the GA template stability (slice-wise) ρGA template for each session (motion versus keeping still) of Exp. 3a and 3b (fMRI/PMC-on and fMRI/PMC-off) to the fMRI/PMC-off, keeping still session using a one-sampled paired t-test. This required 37 tests, one for each time point in the template (k), therefore we applied a Bonferroni correction and considered significant differences at a p-value of b 0.05 corrected.

Results
Experiments 1 and 2: EEG quality assessment inside the scanner with REEGMAS without scanning vs baseline EEG In all subjects, in the moving session (EEG 2M ), the maximum RMS velocity was for translations along the X axis: Vx RMS = 22.5, 23 and 20 mm/s for subjects #1, #2 and #3, respectively (Range: [− 42.9 to 28.6 mm/s], [− 32.2 to 38.6 mm/s] and [− 41.2 to 19.4 mm/s]). The peak frequency of motion was 0.43, 0.71 and 0.55 Hz, for subjects #1, #2 and #3, respectively. All motion-related model regressors explained a significant amount of variance. The motion-related quantity that explained the most variance in the EEG was the velocity of motion as expected followed by position and squared velocity (see Supplementary material Table S1).
EEG 2S data (recorded during the keeping still sessions of Experiment 2, in scanner without fMRI acquisition) were of high quality with a clearly visible alpha rhythm during epochs of eyes closed (orange regions, Fig. 2A) and the presence of eye-blink artefacts (blue region, Fig. 2A). The REEGMAS correction improved the visual appearance EEG 2S attenuating BCG artefacts (red region, Fig. 2A), even prior to applying AAS BCG on EEG 2S-C . REEGMAS correction did not change the visual appearance of alpha rhythm nor the eye-blink artefacts. In the moving session (EEG 2M ), large amplitude voltages were clearly visible during subject movement that can be seen in the MPT tracking derived velocity measurements (Fig 2B, middle panel). Following REEGMAS correction a substantial qualitative improvement in EEG 2M was observed in EEG 2M-C with strong attenuation of the large amplitude motion-related voltages (Fig. 2B). Application of REEGMAS improved the visual appearance of physiological signals such as the alpha rhythm (orange region in Fig. 2B) in the keeping still ( Fig. 2A) and moving (Fig. 2B) sessions. Similar results were obtained for the other two subjects (see Supplementary material Fig. S1 and S2).
In the keeping still session of Experiment 2, the PSDs were similar for the EEG 2S (dashed black curve) and EEG 2S-C (dashed red curve) for all the subjects (Fig. 3A, top row). In general the quantitative analysis (MRMSE) did not show statistic differences (p N 0.05) between EEG 2S and EEG 2S-C (Fig. 4A, see also Supplementary material Table 2a).
We observed an increase in EEG power, predominantly in the frequency range of 0.5-7 Hz during moving sessions (EEG 2M ) for all subjects (Fig. A, bottom row). REEGMAS decreased the power in this frequency range to baseline levels for all subjects. Similar results were obtained in other electrode locations (data from one temporal, parietal, central and frontal electrode are provided in Supplementary material Figures S3-S6). The motion correction decreased the MRMSE significantly (p b 0.05, Supplementary material Table S2), for all subjects (Fig. 4A, right hand side), when compared to the uncorrected EEG data. We also note that residual gradient artefact can be observed in some subjects at the frequency of the slice acquisition and harmonics (13.5 & 27 Hz). In some subjects/conditions the power at these frequencies decreased towards to the levels found before without scanning suggesting AAS GA has provided a better correction of the GA following REEGMAS.
As expected, the electrodes with the largest voltage values in the alpha band (8-12 Hz in our study) are the parietal/occipital electrodes Representative EEG data from subject #2 obtained in Experiment 2 (inside the MRI scanner, no fMRI acquisition). A ten second epoch following the first eyes closed period is shown for the subject keeping still (A) and subject moving (B) sessions. In both A and B the EEG is displayed before (EEG 2S and EEG 2M ) (top) and after (EEG 2S-C and EEG 2M-C ) (bottom) REEGMAS correction. The middle section shows the translational (green cm/s) and rotational (red degrees/s) velocities measured by the MPT motion tracking system. The blue ellipses highlight the eye-blink artefacts; red ellipses highlight the cardiac pulsation-related artefacts and the orange highlight the alpha-rhythm detection. The black arrows point to residuals remaining after motion correction. To visualise the effect of REEGMAS the EEG data presented in A-B were not AAS BCG corrected. when subjects were still (Fig. 5, top row). The alpha rhythm (Fig. 5, middle row) was contaminated by motion induced voltages, mainly in frontal electrodes (subjects 1 and 2) and in temporo-occipital electrodes (subject 3). After REEGMAS, the alpha power was distributed over the occipital channels and more closely corresponded to the topography seen in the still session for all three subjects (Fig. 5

Impact of motion on gradient artefact stability with and without REEGMAS
The AAS GA correction of EEG data (EEG 3aS and EEG 3bS ) acquired in the keeping still sessions (Fig. 6A, top row) resulted in EEG data of good quality for both volume (Fig. 6A, middle row) and slice wise (Fig. 6A, bottom row) templates independent of the fMRI/PMC being off (Fig. 6 left hand side) or on (Fig. 6 right hand side); the alpharhythm was clearly visible (Fig. 6A orange regions) and its frequency's power distribution was comparable to that recorded in Exp. 1 ( Fig. 2B and C, top rows) for all three subjects. We verified that during the sessions. The mean PSD for each subject is presented for Experiment 2 (in-scanner, no fMRI) (A), Experiment 3a (fMRI/PMC-off) (B) and 3b (fMRI/PMC-on) (C). The shaded grey area represents two standard deviations from the mean baseline (Exp. 1) spectra obtained outside the MRI scanner. The black dashed line refers to the mean power spectra of the EEG acquired in-scanner (Experiments 2 and 3) but not corrected by REEGMAS and the red dashed line to the same data corrected by REEGMAS. Electrode O2 is shown due to its sensitivity to detect the alpha-rhythm. keeping still sessions (EEG 3aS and EEG 3bS ) the GA template (slice-wise) had small variability (ρGA template ), meaning that it was temporally stable in fMRI/PMC-off (Fig. 7a, bottom, left hand side) and fMRI/PMC-on (Fig. 7a, bottom row, right hand side) acquisitions. As expected, the AAS GA correction of EEG data (EEG 3aM and EEG 3bM ) acquired during the moving sessions (Fig. 6B, top row) resulted in decreased EEG data quality compared to the still sessions when No REEGMAS was applied and standard volume-wise AAS GA correction was used (Fig. 6B second  row). However, the strong residuals in the EEG data (3aM and 3bM) corrected by volume-wise based AAS GA (Fig. 6b, second row) were attenuated (EEG 3aM and EEG 3bM ) when corrected by slice-wise AAS GA (Fig. 6b, third row). Following REEGMAS, AAS GA and AAS BCG , the EEG data (EEG 3aM-C and EEG 3bM-C ) quality was sufficiently improved to show physiological electrical activity, such as the alpha rhythm, in both volume (Fig. 6B, third row) and slice-wise (Fig. 6B, fourth row) based AAS GA , the latter was less contaminated by GA. Furthermore eye-blink artefacts (blue circles in Fig. 6B, second row) were clearly distinguishable from motion events following motion correction (blue circles in Fig. 6B, fourth and fifth rows). The visual and quantitative improvement was observed for all subjects.
The GA templates (volume-wise and slice-wise) variability was dramatically increased during motion (Fig. 7A, first and second rows). The GA template variability was substantially reduced by slice-wise AAS GA compared to volume-wise AAS GA (Fig. 7A, third and second rows respectively). The GA template stability (slice-wise) was then further improved by REEGMAS where visually the stability approaches that seen when the subject was still (c.f. Fig. 7A third row and fourth row respectively). The PSD analyses of the data from Experiment 3 with fMRI  shown results comparable to the EEG acquired in Experiment 2 without fMRI acquisition. For all Subjects there was a significant reduction (p b 0.05) of the MRMSE after applying REEGMAS prior to GA correction of the EEG data (EEG 3aM and EEG 3bM ) acquired in moving sessions of Experiment 3, a and b, (Fig. 4B-C).

Impact of PMC on gradient artefact template stability and resultant EEG data
For the fMRI/PMC-off (Fig. 7B, top row) and fMRI/PMC-on (Fig. 7B, bottom row) scans, the RMS variance across the 3000 artefact epochs was higher for the EEG data acquired during moving sessions (EEG 3aM and EEG 3bM ) than the EEG data acquired during still sessions (EEG 3aS and EEG 3bS ). When the REEGMAS was applied (EEG 3aM-C and EEG 3bM-C ) the variance was strongly attenuated for both acquisitions; fMRI/ PMC-off (Fig. 7B, top row) and fMRI/PMC-on (Fig. 7B, bottom row).
For all subjects and scan conditions (fMRI/PMC-off and fMRI/PMCon), there was a significant reduction (p b 0.05 corrected for multiple comparisons) in ρGA template due to REEGMAS motion correction ( Table 1). The GA template (slice-wise) variance reduction after applying REEGMAS varied between 62.6% (Subject #2, fMRI/PMC-on scan) and 81.89% (Subject #3, fMRI/PMC-on scan). The variance of motioncorrected EEG (EEG 3aM-C and EEG 3bM-C ) and still EEG data (EEG 3aS and EEG 3bS ) was not statistically different in all the subjects for fMRI/PMC-on acquisitions and for subjects #1 and #2 for fMRI/PMC-off acquisitions, in subject 3 there was a significant difference.

Discussion
In summary, we have proposed a method for reducing motionrelated artefacts in EEG data acquired inside the MR-scanner. This method is based on recording subject's head position using a MPT system and using this information to model and supress the motioninduced voltages in the EEG data. Our main finding is that the proposed method is able to considerably reduce motion-induced voltages both during small and large amplitude movements in the absence (Fig. 2 and Supplementary material S1-S2) and presence of fMRI data acquisition (Fig. 6). The motion-induced voltage correction was shown to increase the visual quality of the EEG data and the alpha-rhythm power topography (Fig. 5). This was confirmed by quantitative analysis in the frequency domain where the MRMSE (Fig. 4) showed a consistent improvement in EEG data quality in 3 subjects when they were moving. The PSD analyses (Fig. 3 and Supplementary material Fig. S3-S6) showed that the motion correction is beneficial in frequencies between approximately 0.5-8 Hz (Fig. 3). In this study, motion information at frequencies above 11 Hz were filtered based on the assumption that the camera provides more accurate tracking at lower frequencies and that head motion is predominantly in this range. While the motion task was predominantly low frequency, higher frequency events were Fig. 6. EEG-fMRI data after AAS GA based on volume or slice-wise templates and AAS BCG correction before (top) and after (bottom) REEGMAS. All figures show a 10 s epoch of the first eyes closed period, for Subect#1. The fMRI/PMC-off acquisition is presented on the left hand side and fMRI/PMC-on the right hand side. A) Data acquired in keeping still sessions. Top panel illustrates the RMS of translational (transl) (blue) and rotational (rot) (red) velocities, middle panels illustrate the EEG data after volume-wise AAS GA correction and the bottom panels the EEG data after slice-wise AAS GA correction. B) Top panel illustrate the RMS of translational (blue) and rotational (red) velocities. Middle panels illustrate the EEG data after volume-wise or slice-wise AAS GA correction but without motion-artefact correction (Before REEGMAS). Bottom panels illustrate the EEG data after volume-wise or slice-wise AAS GA correction and motion-artefact correction (After REEGMAS). The blue regions identify the eye-blink artefacts and the orange highlight the alpha-rhythm detection. also corrected such as those relating to cardiac pulsation-related artefacts. Previous studies (Bonmassar et al., 2002;Masterton et al., 2007;Moosmann et al., 2009;LeVan et al., 2013) have shown the potential for correcting motion-induced voltages induced by small movements (in the order of 1 mm) by linearly modelling motion-induced voltages based on head motion detection. Here we demonstrated for the first time the potential correction of motion-induced voltages in the order of 10 mm, 6 degrees and velocities up to 54 mm/s (See Video 1 in  (Bonmassar et al., 2002;Masterton et al., 2007;Chowdhury et al., 2014) presented approaches that may be able to capture non-linear effects (e.g. scalp pulsation) that can contribute to pulsation-related artefact in EEG data. However, the dominant contribution is due to rigid body rotation (Yan et al., 2009) which REEGMAS corrects and the non-linear terms (velocities squared) in the REEGMAS model also explained significant variance. EEG in the MRI scanner is degraded by large amplitude voltages induced by magnetic field gradients. These can be effectively corrected when they are temporally stable using template subtraction methods (Allen et al., 1998;Allen et al., 2000). The GA template stability can be degraded by subject motion, scanner instabilities (Goldman et al., 2000) or magnetic field gradient instabilities. However on most modern MRI systems the latter contribution is small. Subject motion affects the GA template used for the AAS GA method in two ways. Firstly, motioninduced voltages are added to those induced by the switching gradients, therefore modifying the GA template. The voltage induced by motion is of comparable magnitude to GA voltages (480 μV was recorded in this study for a motion with velocity of 20 mm/s and a rotation of 6°around the x-axis at electrode T8, subject #1 while the GA voltage in this case was 50 to 9000 μV). Secondly, if the subject is in a different spatial location the magnetic field gradient experienced within the EEG circuit is different and correspondingly so is the voltage induced. From previous theoretical work (Yan et al., 2009) the expected alteration in GA is spatially varying but will be 75 μV for 1mm translation on z-axis and 370 μV for a 5mm translation at electrode T8, assuming azimuthal angle of 6°b etween the head and B0 (z-axis). Our experimental data that showed a maximum volume-volume GA template difference in voltage of 382 μV for a head translation on x-axis of 2.44 mm (Exp. 3a, Fig. 6B, middle panel, left hand side). For a slice-wise template both head displacement and correspondingly maximum slice-slice GA template voltage change were reduced to 120 μV and 0.89 mm respectively. For these experiments GA template epoch differences are due to the summation of motion induced voltages and GA artefact differences as explained above. After applying REEGMAS, where motion induced voltages are mainly removed the residual variability is predominantly due to GA changes associated with shifts in spatial position. In this case the maximum volume-volume GA difference for the same conditions (Exp. 3a) was 107.7 μV and slice-slice GA of 1.06 μV, broadly comparable to the expected values (Yan et al., 2009). Therefore the error in the GA template found experimentally due to GA variability related to changes in head position are modest provided a slice-wise template is used.

Supplementary material). Previous studies
For fMRI with the PMC 'on', the magnetic field gradients are updated based on head motion. In this case when the subject moves, the axes of the magnetic field gradients are maintained in the same plane relative to the EEG circuit. In some special cases the circuit should experience the same magnetic field gradients from epoch to epoch as long as the EEG circuit and head move in perfect harmony. In this scenario the PMC updates to the acquisition should contribute to increasing the GA stability. However, translations along each axis due to subject motion will cause a different magnetic field gradient to be experienced from epoch to epoch in either case of the fMRI/PMC being On or Off. Overall this suggests there should not be a significant penalty in GA variability from PMC gradient updates and even be an advantage. We therefore conducted simulations based on Yan et al., 2009 (see Supplementary material-simulation for details). The simulations showed that the variance ρGA template is slightly increased for some of the 37 points of a slice-wise epoch acquisition for fMRI/PMC-on during motion sessions. However, the mean variance considering the whole GA template epoch is not significantly increased for fMRI/PMC-on acquisitions when compared to fMRI/PMC-off acquisitions. This confirmed our experimental findings (Fig. 7B and Table 1) that despite possible confounding differences in voluntary movement between sessions suggested PMC-on and off sessions had similar GA stability.
In practice not all of the EEG equipment will move with the head i.e. there will be increased GA instability from the cables connecting the amplifier to the cap electrodes. This necessitates the use of an alternative approach for minimising GA variance during movement; we have demonstrated that a slice based GA template approach is an effective solution requiring only equally temporally spaced slice acquisition during the volume acquisition period. Our data does confirm that there was not a significant penalty in EEG quality for using PMC with the large potential benefits of increased fMRI image stability. This is crucial because it suggests that PMC can be used for EEG-fMRI data acquisition to improve fMRI and also EEG data quality during motion. Moreover, the amplitudes of the motion-related EEG artefacts reduced by our approach are on the same order (~10mm and 6°) as the range of motion-related artefacts reported to be corrected by the PMC camera-system when applied to fMRI acquisitions (Todd et al., 2015). Therefore, our results represent the first successful implementation of motion-related artefact correction for both EEG (offline correction) and fMRI (prospective) by using the same MPT system to monitor subject motion.
For motion correction methods a key requirement is that data quality is not reduced in the most compliant subjects. In our group when the subjects were keeping still, our method did not change the alpharhythm wave shape (orange region, Fig. 2a) or modify the eye-blink artefact (blue region, Fig. 2a). However, the greatest data quality benefit was obtained during large amplitude motion, where the motionrelated artefact correction improved the representation of the alpharhythm signal (orange region, Fig. 2b), and preserved or recovered the eye-blinks (blue regions, Fig. 2b). This improvement is clearly illustrated by the topographic distribution of the alpha-rhythm, which after applying motion correction (Fig. 5c), is comparable to that expected (Fig. 5a). More advanced analysis such as source localisation of in-scanner EEG (Grova et al., 2008;Vulliemoz et al., 2009;Maziero et al., 2015) could be severely affected by these motion induced changes in topography. At the same time, the motion parameters (red and green tracers, Fig. 2b) provide a useful aid for visual analysis of EEG data that is frequently used for studies of epilepsy (Salek-Haddadi et al., 2006;Gotman et al., 2006;Thornton et al., 2010), where the motion parameters displayed alongside the EEG can be used to determine the likely origin of EEG changes and thereby facilitate interpretation. The ability to scan and obtain EEG of sufficient quality during ictal events, which are clinically important but typically suffer from motion (Chaudhary et al., 2012b), should improve the yield and quality of these studies e.g. by allowing the ictal phases to be recorded with sufficient reliability while supressing motion related variance in the fMRI time series.
EEG acquired in the scanner suffers from a characteristic cardiacrelated artefact (often termed the BCG artefact). The cardiac-related artefact was attenuated by our approach while subjects were keeping still. A recent study (LeVan et al., 2013), used similar methodology to ours (camera-tracker system) but applied it exclusively to correct the cardiac-related artefact. In our results, we verified that the motioncorrection attenuated the cardiac-related artefact during 'keeping still' sessions (red regions, Fig. 2a and Supplementary Figures S1-S2). Our method cannot correct for non-rigid components of the BCG artefact, however realistic modelling of the artefact suggests that bulk head motion is its dominant cause (Yan et al., 2010). We have additionally shown how the system can be used during fMRI motion correction both in terms of effective motion-related voltage suppression and the maintenance of GA temporal stability. In contrast to our approach, LeVan et al. (2013) attached the marker on subject's forehead. The forehead is a region susceptible to involuntary movements such as eye blinks, which would be recorded as head motion, increasing tracking noise. We developed a marker fixation device based on a dental retainer thereby limiting the marker movement to movement of the skull. This may explain why in our study we were able to correct for a wide range of motion amplitudes and frequencies.
Marker fixation has been extensively discussed in the literature (for more details see Maclaren et al., 2013) due to the importance of accurate tracking information for PMC to be beneficial to MRI, and the teeth have been suggested as the best place to attach the marker to. Previous motion tracking studies (Masterton et al., 2007) involved the use of additional circuitry, such as loops and wires in close proximity to the subject. In our acquisitions the camera is spatially distant and electrically isolated from the subject (being mounted on the top of the bore) and is therefore unlikely to pose an increased safety risk (i.e. of localised subject heating (Lemieux, 1997)).
In our data following motion correction, residual artefacts were visible in the corrected EEG data (e.g. black arrows, Fig. 2), during very high amplitude movements. This residual artefact could be due to nonlinear effects when the cap and head do not move together, for example. The residual artefacts appear to occur across channels and so additional spatial (or cross channel) constraint on the motion related artefact removal may further improve the motion artefact correction. There are also limitations to the tracking especially for faster and larger amplitude motion  which might be reduced by using multiple markers and outlier detection and removal in the tracked data. It must be emphasised that these remarks relate to extreme motion, which might be considered by some sufficient grounds to discard entire datasets. It is possible that these residual artefacts are due to GA instability although motion related voltages were in general the source of greater magnitude voltage changes in our data. In the context of the study of difficult-to-scan subjects (due to disease or mental age or other factors making it difficult for them to control their movement), our method offers the possibility of offering a more robust scanning protocol, thereby optimising the use of patient and scanner time.
The ability to record motion-insensitive EEG-fMRI data should enable the improved study of epilepsy patients during seizures and more generally populations both healthy and with pathology where combined recordings of EEG-fMRI can provide important information about brain function, especially considering that motion differences between populations (Satterthwaite et al., 2012) can lead to bias. Additionally, studies related to different neuroscience areas, such as speech production and motor tasks should also benefit from our approach. ) and corrected EEG during moving (EEG 3aM-C and EEG 3bM-C ) and also for the keeping still session for both fMRI/PMC-off and fMRI/PMC-on scans. The brackets indicate which comparisons were statistically significant different (p(bonferroni) b 0.05), where the '*' indicates the higher value.