NeuroImage 280 (2023) 120353

Contents lists available at ScienceDirect

NeuroImage
journal homepage: www.elsevier.com/locate/ynimg

NeuXus open-source tool for real-time artifact reduction in simultaneous

EEG-fMRI
Gustavo Caetano , Inês Esteves , Athanasios Vourvopoulos , Mathis Fleury , Patrícia Figueiredo *
ISR-Lisboa/LARSyS and Department of Bioengineering, Instituto Superior Técnico – Universidade de Lisboa, Lisbon, Portugal

A R T I C L E I N F O A B S T R A C T

Keywords: The simultaneous acquisition of electroencephalography and functional magnetic resonance imaging (EEG-fMRI)
EEG-FMRi allows the complementary study of the brain’s electrophysiology and hemodynamics with high temporal and
Real-time spatial resolution. One application with great potential is neurofeedback training of targeted brain activity, based
Artifact reduction
on the real-time analysis of the EEG and/or fMRI signals. This depends on the ability to reduce in real time the
LSTM
Toolbox
severe artifacts affecting the EEG signal acquired with fMRI, mainly the gradient and pulse artifacts. A few
methods have been proposed for this purpose, but they are either slow, hardware-dependent, publicly unavai­
lable, or proprietary software. Here, we present a fully open-source and publicly available tool for real-time EEG
artifact reduction in simultaneous EEG-fMRI recordings that is fast and applicable to any hardware. Our tool is
integrated in the Python toolbox NeuXus for real-time EEG processing and adapts to a real-time scenario well-
established artifact average subtraction methods combined with a long short-term memory network for R
peak detection. We benchmarked NeuXus on three different datasets, in terms of artifact power reduction and
background signal preservation in resting state, alpha-band power reactivity to eyes closure, and event-related
desynchronization during motor imagery. We showed that NeuXus performed at least as well as the only
available real-time tool for conventional hardware setups (BrainVision’s RecView) and a well-established offline
tool (EEGLAB’s FMRIB plugin). We also demonstrated NeuXus’ real-time ability by reporting execution times
under 250 ms. In conclusion, we present and validate the first fully open-source and hardware-independent
solution for real-time artifact reduction in simultaneous EEG-fMRI studies.

1. Introduction movement due to the heart pumping the blood towards the head,
expansion of the scalp due to arterial vessel pulsation, and the Hall effect
Since 1999, the simultaneous acquisition of electroencephalography creating voltages between the blood ions in the presence of the magnetic
and functional magnetic resonance imaging (EEG-fMRI) has been used field (Debener et al., 2008; Yan et al., 2010; Mullinger et al., 2013;
to study brain function (Bonmassar et al., 1999). This multimodal Abreu et al., 2022). This artifact is already present inside the static
technique provides two distinct views on neuronal processes by magnetic field, even when no MRI is being acquired. In contrast to GA,
combining measures of electric activity with fine temporal resolution PA exhibits a more irregular pattern, following the cardiac cycle. Besides
(EEG) with measures of the brains’ hemodynamics with fine spatial GA and PA, the EEG recorded inside the MRI scanner is also affected by
resolution (fMRI) (Jorge et al., 2014; Abreu et al., 2018; Ritter and spontaneous head movements (also changing the shape of the GA and
Villringer, 2006). One of the greatest challenges are the artifacts induced PA), as well as the environment, including the powerline, lights, venti­
on the EEG data in the MR environment. The largest in amplitude is the lation and the helium pump (Bullock et al., 2021; Mullinger et al., 2008;
Gradient Artifact (GA), which is induced by the time-varying magnetic Mulert and Lemieux, 2023; Nierhaus et al., 2013; Neuner et al., 2014;
field associated with the image acquisition on the electric circuit formed Rothlübbers et al., 2015).
by the electrodes and the subject’s head (Allen et al., 2000; Niazy et al., Several methods have been developed to reduce these EEG artifacts
2005). Using typical fMRI sequences, this variation is exactly repeated in post-processing (offline). Most of them are based on artifact average
every repetition time (TR). In turn, the Pulse Artifact (PA) is induced by template subtraction (AAS) techniques (Allen et al., 2000; Moosmann
changes occurring through the cardiac cycle, including bulk head et al., 2009; Sun et al., 2009; Zhang et al., 2019; Wan et al., 2006; Freyer

* Corresponding author.
E-mail address: patricia.figueiredo@tecnico.ulisboa.pt (P. Figueiredo).

https://doi.org/10.1016/j.neuroimage.2023.120353
Received 6 April 2023; Received in revised form 21 August 2023; Accepted 28 August 2023
Available online 29 August 2023
1053-8119/© 2023 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-
nc-nd/4.0/).
G. Caetano et al. NeuroImage 280 (2023) 120353

et al., 2009; Sartori et al., 2009; Huang et al., 2012; Ferreira et al., command line: neuxus script.py) (Legeay et al., 2022) (https://github.
2014), which work by averaging a set of segments of a signal with a com/LaSEEB/NeuXus). Each class instance (called node) has, at most,
periodic artifact to obtain an artifact template and subtract it from the an input and an output variable (called ports) and a method (called
signal. Other methods include blind source separation techniques such update) that transforms the data from the input to the output port in a
as principal component analysis (PCA) or independent component particular way (e.g. removing the Gradient Artifact of each channel).
analysis (ICA) (Niazy et al., 2005; Ryali et al., 2009; Acharjee et al., The data inside the ports takes the form of tables (called chunks), with
2014; Negishi et al., 2004; Liu et al., 2012; Li et al., 2017; Jing and rows representing time and columns channels. The data chunks may
Sanei, 2006), filtering (Frigo and Narduzzi, 2014; Hoffmann et al., 2000) contain signals (one per channel) or markers (in which case there is only
or dictionary learning approaches (Xu et al., 2007). Alternatively, one channel and time, and the data value is the name of the marker). By
methods have also been proposed that are based on the measurement of creating nodes and setting the output ports of some as the input ports of
artifact signals using additional sensors (Abreu et al., 2022). This can be others, the user can design a pipeline, and execute it, to make each node
achieved directly by using a complete layer of reference electrodes update sequentially, in a loop, and process data in real time.
duplicating the EEG electrodes but insulated from the brain (Luo et al.,
2014; Chowdhury et al., 2014). In a relatively more practical approach, 2.2. Real-time artifact reduction
a limited number of sensors may be placed in specific locations (or
simply obtained by insulating a few EEG electrodes from the scalp) and The EEG is acquired by the acquisition software and streamed to
the motion artifacts in each electrode location may then be estimated as NeuXus. In the specific case the acquisition software is BrainVision
linear combinations of the sensor signal time-courses (Jorge et al., 2015; Recorder (used here), the data is streamed by RDA in blocks with a fixed
Xia et al., 2014; Masterton et al., 2007; Bonmassar et al., 2002). number of channels but a variable number of time points (which de­
Although the PA and other motion artifacts may be more efficiently pends on the acquisition itself, and was experimentally observed to vary
reduced using sensor-based techniques, and despite the fact that carbon between 50 and 250 points at a 5000 Hz acquisition rate, corresponding
wire sensors have recently become commercially available for this to 0.01 to 0.05 s), and received by NeuXus as chunks. The artifact
purpose, most EEG-fMRI studies are still conducted without the use of reduction is then performed, as illustrated by the flowchart in Fig. 1. The
additional sensors. main parameters that need to be defined for NeuXus artifact reduction
Despite the variety of methods proposed for artifact reduction in algorithm are listed in Table 1; some of those are also indicated in Fig. 1.
EEG-fMRI, most of them cannot be used in real time (online) since they The reduction was run in an HP laptop with an AMD Ryzen 7 3700 U
rely on the analysis of the complete signals and are often computa­ processor with Radeon Vega Mobile Gfx 2.30 Hz integrated graphics and
tionally intensive. The online reduction of EEG artifacts would be 8GB of RAM. Since the LSTM estimation step requires intensive calcu­
essential for applications such as to monitor brain activity during the lations, the PA reduction was called with an option that accelerates it by
experiment (for example, epileptic activity) or to provide neurofeedback using Numba, a just-in-time (JIT) compiler that translates Python code
(for example, for training targeted brain activity). A few methods have to optimized machine code at runtime (Lam et al., 2015).
been proposed to reduce EEG artifacts in real time (Shaw, 2017; Wu
et al., 2016; Mayeli et al., 2016; Garreffa et al., 2003; Steyrl et al., 2018; 2.2.1. Gradient artifact
Purdon et al., 2008; Levitt et al., 2023). The only commercially available The GA reduction starts (from the start of the execution or when a
tool (BrainVision RecView) uses AAS to reduce GA and PA, yet it is specific marker is received, for example marking the start of the fMRI
proprietary (from Brain Products GmbH, Germany). The filtering sequence) by concatenating the incoming chunks into TR-length seg­
method proposed in (Shaw, 2017) filters the GA, but introduces tem­ ments (Fig. 1– Top; Table 1: tr). The segments are averaged to cancel the
poral distortions, and needs further validation before becoming publicly signal components not tied to the TR (physiological signal, and artifacts
available. The PCA method proposed in (Wu et al., 2016) reduces the PA except the GA) and create an average GA template. After a minimum
and is not available as well. The ICA method proposed in (Mayeli et al., number of segments have been averaged (Table 1: min_wins), the part of
2016) still requires a prior GA and PA reduction (it was validated using the template which the current chunk contributed to is selected and
RecView) and introduces a delay from 2 to 4 s. Finally, a few subtracted from the chunk, removing its GA, and a marker is thrown to
sensor-based methods have been proposed to reduce GA and PA in real mark the start of the subtraction. After a maximum number of segments
time ( (Garreffa et al., 2003; Steyrl et al., 2018; Purdon et al., 2008; have been averaged (Table 1: max_wins), as a new segment enters into
Levitt et al., 2023)); however, these require additional (mostly custom) the average, the earliest segment is removed. In this way, the template is
hardware not available for most EEG-fMRI studies. built based on the latest TRs, which makes it adaptive to changes in the
Here, we present a fully open-source tool for the real-time reduction GA (e.g. due to head motion).
of gradient and pulse artifacts in simultaneous EEG-fMRI studies, which
is fast, applicable to any EEG-fMRI setup, and publicly available as part 2.2.2. Pulse artifact
of NeuXus, a real-time EEG processing toolbox developed in Python After GA reduction, the chunks are downsampled to 250 Hz. This
(Legeay et al., 2022). The method is based on the AAS technique and step includes an antialiasing lowpass filter of type Chebyshev type I and
uses a long short-term memory (LSTM) artificial neural network for order 8, such that frequencies above the new Nyquist frequency (125
detecting R peaks in the ECG. We first briefly introduce the NeuXus Hz) are prevented from aliasing and contaminating the band of interest.
toolbox, then describe the real-time artifact reduction method, and then The ECG signal is bandpass-filtered from 0.5 to 30 Hz using a Butter­
the evaluation of its performance, done retrospectively (on previously worth filter of order 4, to remove noise and improve the subsequent R
collected data), and compared with a well-established offline method peak detection.
(EEGLAB’s AAS) and a commercially available online method (Brain­ The PA reduction then starts when it receives the marker from the
Vision RecView). start of the GA subtraction, by concatenating the GA-subtracted,
downsampled and ECG-filtered chunks into a detection window. As
2. Materials and methods the window fills, each concatenated chunk is selected and output
without change (Fig. 1: the last selection is represented by the gray
2.1. NeuXus EEG processing toolbox bracket), until a specific position in the detection window (Table 1: a
stride and a margin from the end), after which the chunks are held until
NeuXus is an open-source toolbox developed in Python that allows the window is complete.
the user to process EEG data in real time, by creating a script, instanti­
ating a set of classes, and executing it with the file neuxus (e.g., in the

2
G. Caetano et al. NeuroImage 280 (2023) 120353

Fig. 1. Flowchart of the real-time artifact reduction algorithm in simultaneous EEG-fMRI recordings. The current chunk (time and channel dimensions were switched
to ease the visualization) (in blue) is averaged with previous chunks (in gray) into a GA template, which is used to subtract the GA from the current chunk. The result
is downsampled, ECG-filtered, and concatenated into a detection window. The R peaks are detected and used as triggers to average the EEG across cardiac cycles into
a PA template. Finally, this is subtracted from the last part of the detection window (up to a margin). The GA reduction will keep sending chunks, which will be
pushed to the detection window, and at every certain amount of push (stride), they will trigger a new detection and subsequent PA template average and subtraction.
The parameters “tr”, “max_cycle” and “min_wins” (as described in Table 1) are indicated.

Table 1
Parameters used in the artifact reduction algorithm.
Parameter Description Default value

GA reduction tr Time of repetition of the MRI sequence 1.26 (s)

min_wins Minimum number of averages before subtraction 7
max_wins Maximum number of averages 30
PA reduction det_win_len1 Length of the detection window 500
stride Number of points to enter the filled detection window before new detection 50
margin Duration at the end of the detection window in which data is prevented from being segmented into cardiac cycles and subtracted 0.1 (s)
thres2 Value above which estimated R peak probabilities turn into positive labels 0.05
min_counts Minimum number of neighboring positive labels for a positive label to be considered 5
min_cycle Minimum cardiac cycle duration 0.4 (s)
max_cycle Maximum cardiac cycle duration 1.5 (s)
min_wins Minimum number of averages before subtraction 10
max_wins Maximum number of averages 20
numba Value to specify whether to use Numba for the detection True
filter_ecg Value to specify whether to filter the ECG in the detection window False
1
The length of the detection window is specified in points, instead of duration, since the LSTM detector expects a fixed number of points, independently of the
duration it represents. For our data’s sampling rate of 250 Hz, the length and the stride correspond to 2 s and 0.2 s, respectively.
2
The threshold and minimum counts were taken from the original network without change.

3
G. Caetano et al. NeuroImage 280 (2023) 120353

2.2.2.1. R peak detection. When the window is complete, the cardiac subtraction (Fig. 1: blue bracket). Therefore, only at a detection is a
cycles are found by detecting the R peaks on the corresponding ECG segment output, which ensures it has correctly entered the template and
signal. The detection algorithm is based on Laitala et al. (Laitala et al., been subtracted by it. An example of the real-time GA and PA reduction
2020) (including the LSTM network architecture) and adapted to is shown in Fig. 2. NeuXus also allows the clean data to be sent by using
real-time execution. The LSTM network is composed by two bidirec­ the Lab Streaming Layer (LSL), to be monitored in a viewer (e.g. in
tional LSTM layers and a dense layer. Each layer combines the input StreamViewer) or saved in LabRecorder (https://github.com/sccn/labst
with trainable parameters (132 737 in total) and passes it to an activa­ reaminglayer).
tion function. In the LSTM layers, the input is scanned forward and
backwards, and the activation function is a hyperbolic tangent. In the 2.3. Data acquisition
dense layer, the activation function is a sigmoid, and the output is an
array with the same size as the input to the network. The data used to train the NeuXus LSTM network and evaluate the
The detection starts by normalizing the ECG signal between − 1 and NeuXus artifact reduction algorithm was acquired on a 3T Vida MRI
1, and passing it through the LSTM network, which estimates the scanner (Siemens, Germany), using an MR-compatible EEG system
probability of an R peak at each time point. This probability is then (Brain Products, Germany) with a 32-channel BrainCap, a SyncBox (to
averaged with the probabilities estimated at that point for all windows it synchronize the fMRI and EEG acquisitions), and Triggerbox (to receive
has previously belonged to. The probabilities of the points too close to a marker from the scanner marking the start of the acquisition). Data
the end of the window (within a predefined margin from the end; acquisition was approved by the Comissão de Ética para a Investigação
Table 1: margin) are not considered since the network is assumed not to Clínica of Hospital da Luz, in accordance with the Declaration of Hel­
have enough points beyond that to make an informed decision. These sinki, and all subjects provided written informed consent. All data re­
probabilities are thresholded (Table 1: thres) to yield an excessive ported here will be available upon reasonable request. Validation data is
number of positive R peak labels. The labels then pass through two available at: https://doi.org/10.6084/m9.figshare.22561555.v2. The
filtering steps: (1) only the positive labels surrounded by a minimum cap was placed on the head according to the 10–20 system, and one ECG
number of positive labels are selected (Table 1: min_counts); and (2) if electrode was placed on the back. The data was acquired using Brain­
multiple positive labels are found within one cardiac cycle (Table 1: Vision Recorder (Brain Products, Germany) at a sampling rate of 5000
min_cycle), then only the one with the highest probability is kept. The Hz. The fMRI data was acquired with a 64-channel head RF coil using a
default values of thres, min_counts and min_cycle were kept from Laitala 2D-EPI sequence (TR/TE = 1260/30 ms, in-plane acceleration with
et al. (Laitala et al., 2020). In our dataset, a small variation in the value GRAPPA factor 2, SMS factor 3) to collect 60 axial slices from the whole
of thres (in the order of 10− 2 ) did not lead to a different partition be­ brain with 2.2 mm isotropic resolution.
tween the probabilities of negative and positive labels (in the order of
10− 5 and 10− 1 , respectively). The parameter min_cycle may be tweaked 2.3.1. Data used to train the LSTM network for R peak detection
depending on the subject and task, but min_counts should be tweaked The NeuXus LSTM network was trained using ECG data collected
with care since it coincides with a training parameter (the number of from 6 female volunteers (21 to 41 years) during a total of 32 min per
positive labels the LSTM network is trained to find per peak). Finally, the subject while they were performing a series of cognitive tasks concom­
filtered positive R peak labels are used to segment the EEG signals of the itantly with simultaneous EEG-fMRI scanning. The network was trained
detection window into cardiac cycle segments, and the window segment in 10 epochs, with 40 batches of 256 windows spanning 2 s (which sets
that was previously held is selected for subtraction (Fig. 1: blue bracket, the detection window to 2 s as well), randomly collected from the ECG of
until the margin). this dataset.

2.2.2.2. Average artifact subtraction. After every detection, the cardiac 2.3.2. Data used to evaluate the GA and PA reduction
cycle segments are averaged to create a PA template. To accommodate The NeuXus artifact reduction algorithm was evaluated retrospec­
variations in heart rate, a maximum cardiac cycle duration is set for the tively on data collected from three groups of healthy subjects, each
template (Table 1: max_cycle), such that it may include a different under a different condition: resting state (RS), eyes open / eyes closed
number of segments for each part of the cycle. When a detection window (EO/EC) task, and motor imagery (MI) task. The number of subjects,
segment is selected for subtraction, the corresponding part in the tem­ runs and the duration of each run are specified in Table 2. In the RS
plate, which must have been averaged using a minimum number of condition, subjects were instructed to fixate a white cross on a black
cardiac cycles (Table 1: min_wins) (Fig. 1: dark-orange part), is selected screen for 2 min. For each subject, a second RS EEG acquisition was
and subtracted from the segment, removing its PA. The PA-reduced performed inside the MRI scanner but without concurrent fMRI to
segment is then output, and a marker is thrown to mark the start of capture only the PA without the GA, and evaluate the PA correction. In
the PA subtraction. If the corresponding part in the template has a sec­ the EO/EC task, subjects were instructed to fixate a white cross on a
tion that has not yet been averaged using a minimum number of cycles, black screen for 2 min and close their eyes for another 2 min. In the MI
that section will not be subtracted, and the output segment will be task, subjects looked at a screen showing a first-person perspective of a
partially PA-reduced. In fact, it is preferable not to subtract the PA at person on a boat with two rows (a scenario called NeuRow (Vourvo­
times when the template is not accurate enough. This may happen when poulos et al., 2016)). The paradigm is illustrated in Fig. 3: at each trial
encountering an abnormally long cycle, either naturally or because the (18 trials per run), one of the arms starts rowing and the subjects are
detector missed an R peak (in which case the values averaged incorrectly instructed to imagine they are the ones moving it.
into a distant part of the template are prevented from subtracting).
Nonetheless, these cases are uncommon and can be minimized by 2.4. Evaluation
decreasing the minimum number of averages or increasing the
maximum number of averages (Table 1: max_wins) (in the resting state 2.4.1. R peak detection
dataset 14% of the points were not PA-reduced in the manner described NeuXus’ LTSM R peak detection algorithm was evaluated on the EO/
above). EC dataset by using the F1 score through comparison with the ground
The new chunks that enter the detection window will shift the ones truth (manual peak detection, performed with the tool available at: htt
there and remove the earliest. Every time there has been a certain shift ps://github.com/LaSEEB/interactiveQRS). Its performance was
(Table 1: stride), the LSTM detector is applied to the ECG, the cardiac compared with that of the original LTSM offline algorithm (Laitala et al.,
cycle segments are updated, and the next window segment is selected for 2020) (referred here as Laitala), and the commonly used Pan-Tompkins
algorithm (Pan and Tompkins, 1985). The Kruskal-Wallis test was used

4
G. Caetano et al. NeuroImage 280 (2023) 120353

Fig. 2. Example of the real-time artifact

reduction in simultaneous EEG-fMRI re­
cordings, showing: i) GA reduction on the EEG
C3 channel and ECG signals after 7.56 s (six
1.26 s TR’s), and ii) PA reduction on the EEG C3
signal 2.4 s afterwards (the time taken to fill a
2 s detection window, in green, and slide it two
times with a stride of 0.2 s), when two cardiac
cycles were successfully detected (the PA min­
imum windows was lowered to 2 for the illus­
tration). It’s worth noting the PA-reduced
chunks were output until a 0.1 s margin from
the last detection (the yellow PA-reduced signal
does not reach the end of the detection last
window), and the chunks after will be PA-
reduced and output at the next detection.

2.4.2. Artifact reduction

Table 2
The proposed real-time GA and PA reduction algorithm (referred
Datasets used to evaluate the GA and PA reduction algorithm.
here as NeuXus) was evaluated in comparison with the only commer­
Task Number of Number of Runs per Duration cially available real-time tool (RecView), and the well-established and
Subjects Subject [s]
commonly used offline reduction method based on the same algorithm
Resting state (RS) 8 1 120 (EEGLAB AAS correction, in MATLAB) (Delorme and Makeig, 2004). In
1 (without fMRI) 120 RecView, the MRI artifact and the Pulse-artifact filters were used to
Eyes open / eyes closed 11 1 240
(EO/EC)
reduce the artifacts. In the MRI filter, the TR was set to 1260 ms (Table 1:
Motor imagery (MI) 13 3 323 tr), and in the PA filter, the number of pulses in the average was set to 30
(Table 1: PA: max_wins). In EEGLAB, the fMRIb plugin was used (Ian­
netti et al., 2005). In the GA reduction, the TR and number of windows
in the GA average were set to 1260 ms and 20 (Table 1: GA: max_wins),
respectively, and the options to perform OBS and ANC were disabled. In
the PA reduction, the number of windows in the PA average was set to
30, and the reduction method was set to “mean”, to perform exclusively
AAS.
The GA and PA reduction stages were evaluated individually using
the RS data, by comparing the EEG spectral power in frequency bands
corresponding to the artifact (artifact bands) with that in bands
considered free from artifact (background bands), before and after each
reduction stage, following the method described in detail in (Abreu
et al., 2016). The spectral power was computed using the fast Fourier
transform (FFT) as implemented in MATLAB. For the GA, the artifact
bands were defined as 0.04 Hz-wide intervals centered around 1/TR
(1/1.26 Hz) and its harmonics up to 30 Hz, whereas the background
bands were defined by taking 50% of the interval between consecutive
Fig. 3. Schematic of the Motor Imagery paradigm with the corresponding
artifact bands (Fig. 4- top). For the PA, the artifact bands were defined
NeuRow display. around the PA peaks, which were identified by fitting a Lorentzian in the
vicinity of the average heart rate and its harmonics, as described in
detail in (Abreu et al., 2016), whereas the background bands were
to determine statistically significant differences between these methods.
defined in the same way as for the GA (Fig. 4- bottom).
Whenever the null hypothesis was rejected (at a significance level of
The overall artifact reduction (including both GA and PA) was
p<0.05), pairwise comparisons between methods were performed using
evaluated using the EO/EC and MI datasets, by analyzing its impact on
Dunn’s test.
typical EEG measures of reactivity to eyes closure (increase in alpha
power over the occipital cortex) and motor cortex activation (even-

5
G. Caetano et al. NeuroImage 280 (2023) 120353

Fig. 4. Illustration of the artifact (green) and background (red) bands, for GA (top) and PA (bottom), on the frequency spectrum of a representative example (subject
1 and channel Oz of the RS dataset). The bands are shown superimposed on the spectra before artifact reduction, and after reduction by each of the methods (NeuXus,
RecView, EEGLAB). A zoomed-in section of the top row is displayed underneath to clarify the differences of the frequency spectra between methods.

related reduction in alpha power over the motor cortex), respectively. 2.4.3. Timing
For the EO/EC dataset, the ratio of the EEG power in the alpha band Finally, to evaluate the method’s ability to perform in real time, the
(8–12 Hz) in channel Oz between EO and EC epochs (μEOα and μECα , time that each data point takes to go through each stage of the real-time
respectively) was calculated before and after the reduction, according artifact reduction in NeuXus was measured for one illustrative RS run.
to: Before and after each stage (GA reduction, downsampling, ECG filtering
μECα − μEOα and PA reduction), each chunk’s data points were timestamped with the
RatioEO/ECα % = × 100% (1) current time (measured by time.perf_counter()). In practice, only the
μEOα
chunk rows (time instances) were timestamped, since for every time
For the MI dataset, for each subject, the data of each run was cor­ instance the number of columns (channels) is constant and the values
rected individually, and then concatenated and preprocessed by band­ are assumed to be concurrent. The time instances and timestamps were
pass filtering from 1 to 40 Hz, interpolating bad channels (using EEGLAB saved for each moment in the pipeline. After NeuXus execution, the time
tools), and referencing to the channels’ average. The left arm trials were instances from before and after each stage were matched, and the cor­
analyzed and channel C4 was selected (as it corresponds to the left responding timestamps subtracted, to yield the time each data point
primary motor cortex and is thus commonly used to measure EEG ac­ took in that stage. These times were then used to calculate the median
tivity associated with movement on the right side of the body) and the and 25% and 75% percentiles. By tracking data points (through their
event-related desynchronization (ERD) was computed in the alpha band time instances) instead of complete chunks, the variable length of each
(as this is commonly used as an EEG measure of motor activation). For chunk does not affect the timing, and the points held in the PA reduction
this purpose, the event-related spectral perturbation (ERSP) was (until a detection) only receive their final timestamp when they are
computed as follows (Makeig, 1993): the signal in each trial was output, ensuring their time in that stage is correctly measured.
transformed into a power signal in time-frequency domain using
short-time Fourier transform (STFT) (as implemented in EEGLAB), the 3. Results
power signals were averaged across trials, and the result normalized by
the baseline period according to: 3.1. R peak detection
P(f , t) − PBaseline (f )
ERSP% (f , t) = × 100% (2) The results of the R peak detection are shown in Fig. 5. NeuXus’
PBaseline (f )
LSTM network classification achieved a higher F1 score than Pan-
where f is frequency, t time, P(f.t) is power in frequency f at time t, and Tompkins (p < 0.05).
PBaseline is the power in frequency f averaged across the baseline period
(5 s previous to the trial onset). The ERD was obtained by averaging the 3.2. Artifact reduction
ERSP across all frequencies f in the alpha band (8–12 Hz) and all time
points t during the task period (5 s after trial onset). As with the R peak The results of the evaluation of the GA reduction in the RS dataset are
detection, the Kruskal-Wallis test (at a significance level of 0.05) with presented in Fig. 6. In the artifact bands, there were no significant dif­
post-hoc analysis with Dunn’ tests were performed to determine sig­ ferences between methods. In the background bands, both NeuXus and
nificant differences between methods. RecView yielded a smaller power reduction than the offline correction
(3 vs. 1 and 8%, respectively, on average across channels and subjects, p

6
G. Caetano et al. NeuroImage 280 (2023) 120353

Fig. 5. Evaluation of the R peak detection: Left) Example of the R peak detections in a 5 s window, in which NeuXus misses one peak, and Laitala and Pan-Tompkins
miss two peaks. Right) F1 scores for NeuXus, Laitala and Pan-Tompkins methods across all subjects (each color represents one subject).

Fig. 6. Evaluation of GA reduction. Power change in the artifact bands (top) and background bands (bottom), for each method (NeuXus, RecView, EEGLAB): Left)
Topographies (showing the distributions across channels) of the average across subjects; Right) Boxplots (showing the distributions across subjects) of the average
across channels. Different colors represent different subjects.

< 0.05). In general, the reductions followed a similar pattern across results for the MI dataset are presented in Fig. 10. There were no sig­
channels (on average across subjects), except for RecView in the back­ nificant differences in ERD between methods.
ground bands; in this case, a greater decrease was found near CP2, The computation timing results are presented in Table 3. NeuXus
probably due to one subject deviating considerably from the median (as exhibited a median execution time of 0.2304s with Numba, and 0.6407 s
is also evident in the boxplots of channel averages). without Numba (Numba decreased the execution time by 64%). With
The results of the evaluation of the PA reduction in the RS dataset are Numba, the median execution time for GA reduction, downsampling,
presented in Fig. 7. In the artifact bands, both NeuXus and EEGLAB ECG filtering and PA reduction was of 0.0013, 0.0045, 0.0020 and
yielded a stronger power reduction than RecView (64 and 58 vs. 39%, 0.2226s, respectively.
respectively, on average across channels and subjects, p < 0.05). In the
background bands, there were no significant differences between 4. Discussion
methods. All reductions followed a similar pattern across channels and
subjects. This work proposes and validates the NeuXus open-source tool for
The results of the GA and PA reduction in an illustrative time real-time reduction of the gradient and pulse artifacts in simultaneous
segment of the RS dataset are presented in Fig. 8, showing how the EEG-fMRI acquisitions using conventional hardware setups. We showed
reduction impacts the signal in time domain across the different that NeuXus performs at least as well as the offline reduction (EEGLAB)
methods The results in the EO/EC dataset are presented in Fig. 9, and the only commercially available online tool (RecView), in terms of
showing no significant differences in the alpha power reactivity to eyes artifact power reduction with background signal preservation in resting
closure between NeuXus and the other methods, only between EEGLAB state as well as its impact on EEG measures of alpha power reactivity to
and RecView (136 vs. 61%, on average across subjects, p < 0.05). The eyes closure and desynchronization during motor imagery. NeuXus’

Fig. 7. Evaluation of PA reduction. Power change with PA reduction in the artifact bands (top) and background bands (bottom), for each method (NeuXus, RecView,
EEGLAB): Left) Topographies (showing the distributions across channels) of the average across subjects; Right) Boxplots (showing the distributions across subjects) of
the average across channels. Different colors represent different subjects.

7
G. Caetano et al. NeuroImage 280 (2023) 120353

Fig. 8. Effects of artifact reduction on an illustrative time segment of resting-state data (RS dataset): example of the GA and PA-cleaned signals using NeuXus,
RecView and EEGLAB for one illustrative time segment of one subject.

Fig. 9. Effects of artifact reduction on alpha power reduction in eyes closed relative to eyes open conditions (EC/EO dataset): Left) Power frequency spectrum during
EO (blue) and EC (red) conditions for one illustrative subject; the individual alpha band is indicated in gray. Right) Distributions across all subjects of the alpha power
reduction (averaged across channels). Different colors represent different subjects.

Fig. 10. Effects of artifact reduction on the event-related desynchronization (ERD) in the alpha band during motor imagery (MI dataset): Left) ERSP plots obtained
for channel C4 in the left MI trials, in a representative subject, using the three artifact reduction methods. Right) Distributions across all subjects of the ERD, i.e., the
ERSP averaged across the alpha frequency band (8–12 Hz) and the task time interval (0.5–4 s). Different colors represent different subjects.

real-time ability was also validated by achieving total execution times

Table 3
under 250 ms.
Time for a data point to go through each stage of the artifact reduction pipeline,
with and without Numba. Since the downsampling reduces the points, only the
ones kept were used for the measurement. 4.1. Benchmarking of the GA reduction

Without Numba With Numba

In the GA reduction, NeuXus performed close to the other methods in
GA reduction (s) 0.0044 (0.7%) 0.0013 (0.6%) the artifact bands, and better than EEGLAB in the background bands.
Downsampling (s) 0.0101 (1.6%) 0.0045 (2.0%)
NeuXus differs from the EEGLAB method in at least three aspects. Firstly,
ECG filtering (s) 0.0020 (0.3%) 0.0020 (0.9%)
PA reduction (s) 0.6243 (97.4%) 0.2226 (96.6%) it does not align TR segments. EEGLAB upsamples the EEG and segments
Total (s) 0.6407 (100%) 0.2304 (100%) each TR at limits slightly shifted to maximize the correlation between
that and the first TR. This step helps to correct the differences of the GA
in each TR that may be due to the desynchronization between the EEG
and fMRI sampling rates, or the inexact time the scanner may take to
acquire one volume. In real time, this step would cost more time,
memory and increase the algorithm complexity, and since SyncBox was

8
G. Caetano et al. NeuroImage 280 (2023) 120353

used and there was no reason to suspect the time to acquire a volume real time. The Pan-Tompkins searches for high and sharp peaks but fails
deviated from the TR, it was not pursued. Hence, it is recommended to when it finds them on the T wave (augmented by the magnetic field),
use a device to synchronize the acquisitions and input the accurate value even if the location, within the context of the cycle, differs from previous
of the TR. detections.
Secondly, NeuXus does not apply a highpass filter to the data before
building the template. EEGLAB filters the baseline oscillations (< 1 Hz) 4.3. Benchmarking of the PA reduction
to improve the artifact estimation (since they are not expected for the
GA), and then subtracts the template to the unfiltered data, preserving In the PA reduction, NeuXus performed better than RecView in the
them, which depending on the experiment, may be useful to analyze. artifact bands, and close to the other methods in the background bands.
The filter was not implemented in the real-time GA reduction since it Similarly to RecView, NeuXus used the R peaks from its real-time
introduces a delay on the data, which would mismatch the template detection, whereas EEGLAB used manually marked peaks, which could
from the unfiltered data, affecting the subtraction. In an offline scenario, have favored EEGLAB. However, this was not reflected in the results.
a filter with half the order could be applied twice, one forward and NeuXus differs from EEGLAB and RecView by not attempting to capture
another backwards, to cancel the delays. However, in real time the the PA in a particular region of the cardiac cycle. In EEGLAB, this region
complete recording is not available, and applying this strategy on in­ is around the point 0.21 s after the R peak, for an extent based on the
dividual chunks is not viable since it would distort each one at the edges. average cycle, and in RecView it is around the R peak (within a default
A possible solution would be to use a highpass filter with linear phase, interval between 0.1 s before and 0.6 s after the peak). Instead, NeuXus
which applies a constant delay (group delay) to all frequency compo­ averages the full cardiac cycles into the template, and just uses the part
nents (avoiding phase distortion), and then shift the filtered signal back that has been averaged a minimum number of times for subtraction. By
that delay. Linear phase could be achieved with an FIR filter or doing so, there is no assumption regarding the position of the PA relative
approximated with an IIR filter (more computational efficient). Since to the R peak, and the criteria to subtract is based on average strength
this would increase the algorithm complexity, the results were close to alone. It must be noted, however, that if the region closely after the R
the best performing method and the baseline below 1 Hz was not peak has no PA content, it will still likely be used for the subtraction,
considered to be of wide relevance, it was not implemented. It is worth since it will have been averaged a high number of times (since most
noting that a filter can be applied before GA reduction itself, removing cycles include it). The same will happen for other regions, although the
the baseline from the data that becomes part of the template, but also ones most distant from the R peak can be excluded by increasing the
that is subtracted by it. The baseline would not be preserved, and each minimum number of averages (Table 1: min_wins).
point would have a small delay, but the subtraction would not be mis­ Our rationale in the choice of an offline AAS method for the PA
matched, since the filtered data would be used for both the template and reduction was to compare our online method with the closest possible
the subtraction. method performed offline. In the future, it would also be interesting to
Thirdly, NeuXus does not scale the amplitudes of the template to the compared it with PCA-based methods (Niazy et al., 2005). Nevertheless,
chunk. EEGLAB multiplies the template by a constant to minimize the Niazy et al., 2005 compared his offline PCA-based method to an offline
least squares between the template and the data. Although it was tested, AAS method, and so we believe that the performance of our online AAS
it was not found to affect the results. RecView applies a TR alignment method relative to an offline PCA-based method can be inferred from
step as well, by building a user-defined number of templates, only from this.
chunks that are closely aligned, and select the closest to the current The method used to evaluate artifact reduction has some limitations,
chunk for the subtraction. particularly for PA. In fact, while the approach of separating the spec­
After GA reduction, the data are downsampled and the ECG signal is trum into artifact and background bands works well for GA, which is
bandpass-filtered. The filter introduces a delay that misaligns the ECG more spectrally confined, it is less precise for PA, which is more
relative to the EEG. For the QRS frequency range of 8–50 Hz (Ter­ spectrum-wide due to temporal variations in heart rate. Nevertheless,
eshchenko and Josephson, 2015), the average delay is estimated to be we believe that this approximate assessment is useful for the purpose of
0.0012 s. This means that the PA might be segmented at slightly shifted comparing methods, since all of them suffer from the same imprecisions.
R peaks. As with GA, applying the bandpass filter forward and back­
wards to every chunk would distort each chunk’s edges. As an alterna­ 4.4. Timing
tive, it could be applied to the concatenated ECG signal in the R peak
detection (Table 1: filter_ecg). Not only is the ECG longer at this stage NeuXus takes a median of 0.2304s to output a clean data point with
(and so has less edges than the chunks that comprise it), but each edge Numba, and 0.6407 s without. This time was measured from the
also has a low contribution to the artifact reduction (the first does not beginning of the GA reduction until the end of the PA reduction. How­
coincide with the current chunk, and the second lies within the margin). ever, before reaching the GA reduction, the EEG signal must travel
Nevertheless, the PA contains mostly lower frequencies and is not through the acquisition hardware and software. In this experimental
exactly timed by the R peaks anyway. Therefore, we considered that the setup, this includes the amplifier hardware filters, A/D converter, USB
slight delay and resulting mismatch with the EEG would not affect the communication, Windows OS, RDA server (in BrainVision Recorder),
artifact subtraction, and kept the filter before the PA reduction. ethernet connection (TCP/IP), and finally the RDA client (in NeuXus).
Brainproducts reports an average delay of 0.050 to 0.100 s in this
4.2. Benchmarking of the R peak detection pipeline. The acquisition and processing delays should be considered to
estimate how close the experiment can be performed to real time. Still,
In the R peak detection, NeuXus achieved a higher F1 score than Pan- using Brainproducts hardware, there is the possibility of accessing data
Tompkins. NeuXus detection is a close real-time adaptation of Laitala directly from the USB drivers without the use of recording software, and
but differs in three aspects. NeuXus applies a bandpass filter instead of a decrease the acquisition time (https://pressrelease.brainproducts.com/r
smoothing filter to the ECG before the prediction (following the EEGLAB eal-time-eeg/).
offline detection). Moreover, it excludes the probabilities too close to the Most of the processing time is taken in the PA reduction. The original
end of each prediction window. Without this step, the LSTM was shown LSTM network that our tool was based on was created, trained, and used
to predict premature R peaks, or not predict the R peak at all, in cases with Tensorflow (Abadi et al., 2023). The reported time for prediction
when the actual R peak, or the following T wave, respectively, had not was in the order of seconds (even using an Nvidia GTX 1070 GPU, and
yet entered the window. Finally, NeuXus excludes the positive labels too thus restricting the range of computers able to run it). By reproducing
close to others using a rewritten code passage to be clearer and work in the LSTM network in NumPy, we were able to reduce the time to 0.15 s.

9
G. Caetano et al. NeuroImage 280 (2023) 120353

To reduce it further, Numba was used to compile the network at run­ Acknowledgements
time. Numba works by compiling a specified function at the first time it
is called, to generate machine code that is run from the second time We acknowledge the Fundação para a Ciência e Tecnologia (FCT) for
onwards. The first call is slow, and so when using Numba, NeuXus makes financial support through Projects PTDC/CCI-COM/31485/2017 (MF),
a dummy prediction at the start of the execution, when the user must PTDC/EMD-EMD/29675/2017 - LISBOA-01–0145-FEDER-029675
wait (in the test computer, 20 to 60 s). However, in the actual pre­ (GC), PD/BD/150356/2019 (IE), CEECIND/01073/2018 (AV), UIDB/
dictions, the time is reduced to 0.0263 s. Most of the remaining PA 50009/2020 (GC).
reduction delay (median of 0.2 s) is deliberately imposed by holding the We also would like to thank Juho Laitala and his colleagues (Lam
data until the next detection (to correctly classify its position in the et al., 2015) for sharing their work on robust ECG R peak detection as an
cardiac cycle). This delay can be explained by observing that the data is open-source repository.
selected for subtraction from an interval corresponding to a stride plus a
margin (0.2 + 0.1 = 0.3 s) to a margin (0.1 s) from the end of the References
detection window. Furthermore, since most of the time is spent pre­
dicting R peaks and holding data (almost 96.6%), which do not depend Pan, J., Tompkins, W.J., 1985. A real-time QRS detection algorithm. IEEE Trans. Biomed.
Eng. BME-32 (3) https://doi.org/10.1109/TBME.1985.325532.
on the number of channels, NeuXus time is expected to increase only Makeig, S., 1993. Auditory event-related dynamics of the EEG spectrum and effects of
slightly when using a cap with more electrodes. We confirmed this by exposure to tones. Electroencephalogr. Clin. Neurophysiol. 86 (4), 283–293. https://
timing it using a 64-channel cap, which yielded a median of 0.2431s. doi.org/10.1016/0013-4694(93)90110-H. Apr.
Bonmassar, G., Anami, K., Ives, J., Belliveau, J.W., 1999. Visual evoked potential (VEP)
To reduce NeuXus’ processing time, the holding time can be reduced measured by simultaneous 64-channel EEG and 3T fMRI. Neuroreport 10 (9).
by decreasing the margin (although it may yield misclassified R peaks), https://doi.org/10.1097/00001756-199906230-00018.
and/or the stride, which will trigger predictions more frequently Hoffmann, A., Jäger, L., Werhahn, K.J., Jaschke, M., Noachtar, S., Reiser, M., 2000.
Electroencephalography during functional echo-planar imaging: detection of
(although each prediction has a duration, as described above). Another epileptic spikes using post-processing methods. Magn. Reson. Med. 44 (5) https://
way would be to train the LSTM with shorter samples, shortening the doi.org/10.1002/1522-2594(200011)44:5<791::AID-MRM17>3.0.CO;2-2.
detection window, and hence the detection time. Currently, 2 s were Allen, P.J., Josephs, O., Turner, R., 2000. A method for removing imaging artifact from
continuous EEG recorded during functional MRI. Neuroimage 12 (2). https://doi.
used (Laitala used 4 s), but indeed 1 s may be enough to train the
org/10.1006/nimg.2000.0599.
network, since it captures an average cardiac cycle, and reduce the Bonmassar, G., et al., 2002. Motion and ballistocardiogram artifact removal for
detection time by half. interleaved recording of EEG and EPs during MRI. Neuroimage 16 (4), 1127–1141.
https://doi.org/10.1006/NIMG.2002.1125. Aug.
Garreffa, G., et al., 2003. Real-time MR artifacts filtering during continuous EEG/fMRI
5. Conclusion acquisition. Magn. Reson. Imaging 21 (10), 1175–1189. https://doi.org/10.1016/j.
mri.2003.08.019.
This work presents a real-time GA and PA reduction method for EEG Delorme, A., Makeig, S., 2004. EEGLAB: an open source toolbox for analysis of single-
trial EEG dynamics including independent component analysis. J. Neurosci. Methods
data recorded during fMRI, which is fully open source and publicly 134 (1). https://doi.org/10.1016/j.jneumeth.2003.10.009.
available, fast, and hardware independent, offering the first tool Negishi, M., Abildgaard, M., Nixon, T., Constable, R.Todd, 2004. Removal of time-
combining all these properties. The method performs at least as well as varying gradient artifacts from EEG data acquired during continuous fMRI. Clinic.
Neurophysiol. 115 (9) https://doi.org/10.1016/j.clinph.2004.04.005.
the offline correction based on the same algorithm and the only avail­ Iannetti, G.D., et al., 2005. Simultaneous recording of laser-evoked brain potentials and
able hardware-independent real-time reduction (proprietary) software. continuous, high-field functional magnetic resonance imaging in humans.
It presents a necessary short delay that can be configurable. The clean Neuroimage 28 (3), 708–719. https://doi.org/10.1016/J.
NEUROIMAGE.2005.06.060. Nov.
data can then be used to monitor or perform neurofeedback training. Niazy, R.K., Beckmann, C.F., Iannetti, G.D., Brady, J.M., Smith, S.M., 2005. Removal of
The method is part of the open-source Python platform NeuXus, which FMRI environment artifacts from EEG data using optimal basis sets. Neuroimage 28
can be easily used and extended. (3). https://doi.org/10.1016/j.neuroimage.2005.06.067.
Jing, M., Sanei, S., 2006. Scanner artifact removal in simultaneous EEG-fMRI for
epileptic seizure prediction. In: Proceedings of the International Conference on
CRediT authorship contribution statement Pattern Recognition. https://doi.org/10.1109/ICPR.2006.1026.
Wan, X., Iwata, K., Riera, J., Kitamura, M., Kawashima, R., 2006. Artifact reduction for
Gustavo Caetano: Conceptualization, Methodology, Software, simultaneous EEG/fMRI recording: adaptive FIR reduction of imaging artifacts.
Clinic. Neurophysiol. 117 (3) https://doi.org/10.1016/j.clinph.2005.07.025.
Validation, Formal analysis, Investigation, Writing – original draft, Ritter, P., Villringer, A., 2006. Simultaneous EEG-fMRI. Neurosci. Biobehav. Rev. 30 (6)
Writing – review & editing. Inês Esteves: Methodology, Investigation, https://doi.org/10.1016/j.neubiorev.2006.06.008.
Data curation. Athanasios Vourvopoulos: Software, Resources, Masterton, R.A.J., Abbott, D.F., Fleming, S.W., Jackson, G.D., 2007. Measurement and
reduction of motion and ballistocardiogram artefacts from simultaneous EEG and
Investigation. Mathis Fleury: Investigation. Patrícia Figueiredo: fMRI recordings. Neuroimage 37 (1), 202–211. https://doi.org/10.1016/J.
Conceptualization, Methodology, Resources, Writing – review & editing, NEUROIMAGE.2007.02.060. Aug.
Supervision, Project administration, Funding acquisition. Xu, P., Chen, H., Liu, Z., Yao, D., 2007. A new method based on sparse component
decomposition to remove MRI artifacts in the continuous EEG recordings. Sheng Wu
Yi Xue Gong Cheng Xue Za Zhi 24 (2). https://doi.org/10.1109/
Declaration of Competing Interest iembs.2005.1616849.
Purdon, P.L., Millan, H., Fuller, P.L., Bonmassar, G., 2008. An open-source hardware and
software system for acquisition and real-time processing of electrophysiology during
The authors have no conflicts of interest to declare. high field MRI. J. Neurosci. Methods 175 (2). https://doi.org/10.1016/j.
jneumeth.2008.07.017.
Data availability Mullinger, K., Debener, S., Coxon, R., Bowtell, R., 2008. Effects of simultaneous EEG
recording on MRI data quality at 1.5, 3 and 7 tesla. Int. J. Psychophysiol. 67 (3)
https://doi.org/10.1016/j.ijpsycho.2007.06.008.
NeuXus is an open-source software package, and is available at: Debener, S., Mullinger, K.J., Niazy, R.K., Bowtell, R.W., 2008. Properties of the
https://github.com/LaSEEB/NeuXus. In particular, the version used at ballistocardiogram artefact as revealed by EEG recordings at 1.5, 3 and 7 T static
time of publishing is preserved in the release tag v.0.0.4 (https://github. magnetic field strength. Int. J. Psychophysiol. 67 (3) https://doi.org/10.1016/j.
ijpsycho.2007.05.015.
com/LaSEEB/NeuXus/tree/v0.0.4). An example of how to use the arti­ Ryali, S., Glover, G.H., Chang, C., Menon, V., 2009. Development, validation, and
fact reduction tool can be found in the subfolder: /examples/mri-arti­ comparison of ICA-based gradient artifact reduction algorithms for simultaneous
fact-correction/. A validation dataset is provided at: https://doi. EEG-spiral in/out and echo-planar fMRI recordings. Neuroimage 48 (2). https://doi.
org/10.1016/j.neuroimage.2009.06.072.
org/10.6084/m9.figshare.22561555.v2. Sartori, E., et al., 2009. Gradient artifact removal in co-registration EEG/fMRI. In: IFMBE
Proceedings. https://doi.org/10.1007/978-3-642-03882-2_304.

10
G. Caetano et al. NeuroImage 280 (2023) 120353

Freyer, F., Becker, R., Anami, K., Curio, G., Villringer, A., Ritter, P., 2009. Ultrahigh- Jorge, J., Grouiller, F., Gruetter, R., van der Zwaag, W., Figueiredo, P., 2015. Towards
frequency EEG during fMRI: pushing the limits of imaging-artifact correction. high-quality simultaneous EEG-fMRI at 7 T: detection and reduction of EEG artifacts
Neuroimage 48 (1). https://doi.org/10.1016/j.neuroimage.2009.06.022. due to head motion. Neuroimage 120, 143–153. https://doi.org/10.1016/J.
Sun, L., Rieger, J., Hinrichs, H., 2009. Maximum noise fraction (MNF) transformation to NEUROIMAGE.2015.07.020. Oct.
remove ballistocardiographic artifacts in EEG signals recorded during fMRI Rothlübbers, S., Relvas, V., Leal, A., Murta, T., Lemieux, L., Figueiredo, P., 2015.
scanning. Neuroimage 46 (1). https://doi.org/10.1016/j.neuroimage.2009.01.059. Characterisation and reduction of the EEG Artefact caused by the helium cooling
Moosmann, M., Schönfelder, V.H., Specht, K., Scheeringa, R., Nordby, H., Hugdahl, K., pump in the MR environment: validation in epilepsy patient data. Brain Topogr. 28
2009. Realignment parameter-informed artefact correction for simultaneous EEG- (2) https://doi.org/10.1007/s10548-014-0408-0.
fMRI recordings. Neuroimage 45 (4). https://doi.org/10.1016/j. Abreu, R., et al., 2016. Ballistocardiogram artifact correction taking into account
neuroimage.2009.01.024. physiological signal preservation in simultaneous EEG-fMRI. Neuroimage 135.
Yan, W.X., Mullinger, K.J., Geirsdottir, G.B., Bowtell, R., 2010. Physical modeling of https://doi.org/10.1016/j.neuroimage.2016.03.034.
pulse artefact sources in simultaneous EEG/fMRI. Hum. Brain Mapp. 31 (4) https:// Vourvopoulos, A., Ferreira, A., Badia, S.B.I., 2016. NeuRow: an immersive VR
doi.org/10.1002/hbm.20891. environment for motor-imagery training with the use of Brain-Computer Interfaces
Liu, Z., de Zwart, J.A., van Gelderen, P., Kuo, L.W., Duyn, J.H., 2012. Statistical feature and vibrotactile feedback. In: PhyCS 2016 - Proceedings of the 3rd International
extraction for artifact removal from concurrent fMRI-EEG recordings. Neuroimage Conference on Physiological Computing Systems. https://doi.org/10.5220/
59 (3). https://doi.org/10.1016/j.neuroimage.2011.10.042. 0005939400430053.
Huang, C.H., Ju, M.S., Lin, C.C.K., 2012. A robust algorithm for removing artifacts in EEG Mayeli, A., Zotev, V., Refai, H., Bodurka, J., 2016. Real-time EEG artifact correction
recorded during FMRI/EEG study. Comput. Biol. Med. 42 (4) https://doi.org/ during fMRI using ICA. J. Neurosci. Methods 274, 27–37. https://doi.org/10.1016/j.
10.1016/j.compbiomed.2011.12.014. jneumeth.2016.09.012. Dec.
Nierhaus, T., Gundlach, C., Goltz, D., Thiel, S.D., Pleger, B., Villringer, A., 2013. Internal Wu, X., Wu, T., Zhan, Z., Yao, L., Wen, X., 2016. A real-time method to reduce
ventilation system of MR scanners induces specific EEG artifact during simultaneous ballistocardiogram artifacts from EEG during fMRI based on optimal basis sets
EEG-fMRI. Neuroimage 74. https://doi.org/10.1016/j.neuroimage.2013.02.016. (OBS). Comput. Methods Programs Biomed. 127, 114–125. https://doi.org/
Mullinger, K.J., Havenhand, J., Bowtell, R., 2013. Identifying the sources of the pulse 10.1016/j.cmpb.2016.01.018. Apr.
artefact in EEG recordings made inside an MR scanner. Neuroimage 71. https://doi. Shaw, S.B., 2017. Real-time filtering of gradient artifacts from simultaneous EEG-fMRI
org/10.1016/j.neuroimage.2012.12.070. data. 2017 International Workshop on Pattern Recognition in Neuroimaging, PRNI
Xia, H., Ruan, D., Cohen, M.S., 2014. Removing ballistocardiogram (BCG) artifact from 2017. Institute of Electrical and Electronics Engineers Inc. https://doi.org/10.1109/
full-scalp EEG acquired inside the MR scanner with Orthogonal Matching Pursuit PRNI.2017.7981510. Jul.
(OMP). Front. Neurosci. 8 (8 JUL), 218. https://doi.org/10.3389/ Li, J., et al., 2017. A unified canonical correlation analysis-based framework for
FNINS.2014.00218/ABSTRACT. Jul. removing gradient artifact in concurrent EEG/fMRI recording and motion artifact in
Chowdhury, M.E.H., Mullinger, K.J., Glover, P., Bowtell, R., 2014. Reference layer walking recording from EEG signal. Med. Biol. Eng. Comput. 55 (9) https://doi.org/
artefact subtraction (RLAS): a novel method of minimizing EEG artefacts during 10.1007/s11517-017-1620-3.
simultaneous fMRI. Neuroimage 84, 307–319. https://doi.org/10.1016/J. Steyrl, D., Krausz, G., Koschutnig, K., Edlinger, G., Müller-Putz, G.R., 2018. Online
NEUROIMAGE.2013.08.039. Jan. reduction of artifacts in EEG of simultaneous EEG-fMRI using reference layer
Luo, Q., Huang, X., Glover, G.H., 2014. Ballistocardiogram artifact removal with a adaptive filtering (RLAF). Brain Topogr. 31 (1) https://doi.org/10.1007/s10548-
reference layer and standard EEG cap. J. Neurosci. Methods 233, 137–149. https:// 017-0606-7.
doi.org/10.1016/J.JNEUMETH.2014.06.021. Aug. Abreu, R., Leal, A., Figueiredo, P., 2018. EEG-informed fMRI: a review of data analysis
Frigo, G., Narduzzi, C., 2014. EEG gradient artifact removal by compressive sensing and methods. Front. Hum. Neurosci. 12 https://doi.org/10.3389/fnhum.2018.00029.
Taylor-Fourier transform. In: Proceedings of the IEEE MeMeA 2014 - IEEE Zhang, S., Hennig, J., Levan, P., 2019. Direct modelling of gradient artifacts for EEG-
International Symposium on Medical Measurements and Applications, Proceedings. fMRI denoising and motion tracking. J. Neural Eng. 16 (5) https://doi.org/10.1088/
https://doi.org/10.1109/MeMeA.2014.6860079. 1741-2552/ab2b21.
Acharjee, P.P., Phlypo, R., Wu, L., Calhoun, V.D., Adali, T., 2014. Gradient artifact Laitala, J., et al., 2020. Robust ECG R-peak detection using LSTM. In: Proceedings of the
removal in concurrently acquired EEG data using independent vector analysis. In: ACM Symposium on Applied Computing. https://doi.org/10.1145/
Proceedings of the ICASSP, IEEE International Conference on Acoustics, Speech and 3341105.3373945.
Signal Processing - Proceedings. https://doi.org/10.1109/ICASSP.2014.6854727. Bullock, M., Jackson, G.D., Abbott, D.F., 2021. Artifact reduction in simultaneous EEG-
Ferreira, J.L., Aarts, R.M., Cluitmans, P.J.M., 2014. Optimized moving-average filtering fMRI: a systematic review of methods and contemporary usage. Front Neurol 12.
for gradient artefact correction during simultaneous EEG-fMRI. In: Proceedings of https://doi.org/10.3389/fneur.2021.622719.
the ISSNIP Biosignals and Biorobotics Conference, BRC. https://doi.org/10.1109/ Legeay, S., Caetano, G., Figueiredo, P., Vourvopoulos, A., 2022. NeuXus: a biosignal
BRC.2014.6880955. processing and classification pipeline for real-time brain-computer interaction. In:
Neuner, I., Arrubla, J., Felder, J., Shah, N.J., 2014. Simultaneous EEG-fMRI acquisition at MELECON 2022 - IEEE Mediterranean Electrotechnical Conference, Proceedings,
low, high and ultra-high magnetic fields up to 9.4T: perspectives and challenges. pp. 424–429. https://doi.org/10.1109/MELECON53508.2022.9842925.
Neuroimage 102 (P1). https://doi.org/10.1016/j.neuroimage.2013.06.048. R. Abreu, J. Jorge, and P. Figueiredo, “EEG Quality: the Pulse Artifact,” EEG - fMRI, pp.
Jorge, J., van der Zwaag, W., Figueiredo, P., 2014. EEG-fMRI integration for the study of 167–188, 2022, doi: 10.1007/978-3-031-07121-8_8.
human brain function. Neuroimage 102 (P1). https://doi.org/10.1016/j. M. Abadi et al., “TensorFlow: large-Scale Machine Learning on Heterogeneous
neuroimage.2013.05.114. Distributed Systems”, Accessed: Mar. 22, 2023. [Online]. Available: www.tensorflow
Tereshchenko, L.G., Josephson, M.E., 2015. Frequency content and characteristics of .org.
ventricular conduction. J. Electrocardiol. 48 (6), 933. https://doi.org/10.1016/J. Levitt, J., Yang, Z., Williams, S.D., Lütschg Espinosa, S.E., Garcia-Casal, A., Lewis, L.D.,
JELECTROCARD.2015.08.034. 2023. EEG-LLAMAS: a low-latency neurofeedback platform for artifact reduction in
Lam, S.K., Pitrou, A., Seibert, S., 2015. Numba: a LLVM-based Python JIT Compiler. In: EEG-fMRI. Neuroimage 273. https://doi.org/10.1016/J.
Proceedings of the LLVM-HPC2015: 2nd Workshop on the LLVM Compiler NEUROIMAGE.2023.120092. Jun.
Infrastructure in HPC - Held in conjunction with SC 2015: The International C. Mulert and L. Lemieux, EEG - fMRI: physiological basis, technique, and applications.
Conference for High Performance Computing, Networking, Storage and Analysis. 2023. doi: 10.1007/978-3-540-87919-0.
https://doi.org/10.1145/2833157.2833162.

11