BrainBeats as an Open-Source EEGLAB Plugin to Jointly Analyze EEG and Cardiovascular Signals

Informed consent was obtained from each participant, and the Ural Federal University ethics committee approved the experimental protocol.

1. BrainBeats requirements

1. Install MATLAB and EEGLAB on the computer. EEGLAB can be downloaded at https://github.com/sccn/eeglab and unzipped (or cloned for Git users) anywhere on the computer. See the GitHub page for more details on installation. 
2. Add the path to the EEGLAB folder in MATLAB's home panel by clicking Set Path button. Select the eeglab folder with the unzipped file and click Save > Close.
3. Launch EEGLAB by typing eeglab in MATLAB's command window.
4. Install the BrainBeats plugin by clicking File > Manage EEGLAB extensions. Type brainbeats in the search bar, select the BrainBeats plugin in the list and click Install/Update.
5. Load the sample dataset into EEGLAB. Click File > Load existing data. Navigate to the EEGLAB folder, go to the plugins folder, go to the BrainBeats folder and open to the sample_data folder. Select the file dataset.set.
   NOTE: This dataset corresponds to sub-032 (resting state with eyes opened) from an open-source multi-subject dataset^62,63. This was selected because it contains simultaneous EEG (64 channels), ECG (one channel), and PPG (one channel), at a 1000 Hz sampling rate. ECG and PPG signals were collected using the auxiliary inputs of the actiCHamp system. EEG data were recorded with active electrodes placed according to the 10-20 system, with FCz as the online reference and Fpz as the ground electrode, and with the impedance maintained below 25 kOm. ECG was recorded using one active electrode placed on the right wrist, the reference electrode on the left wrist, and the ground on the left inner forearm at 3 cm distally from the elbow. PPG was recorded from the left index finger. EEG, ECG, and PPG data are time-synchronized since they were recorded simultaneously with the same amplifier. See references^62,63 for more details.
6. Check the Save Outputs checkbox to save everything in the corresponding .set file at the same location as the original file that was loaded into EEGLAB.
   NOTE: For this tutorial, we merged the EEG, ECG, and PPG data into one EEGLAB dataset, loaded the 3D boundary element method (BEM) electrode coordinates, and downsampled the signal to 250 Hz to reduce the file weight (for user download) and accelerate computing time. Because this dataset had no abnormal EEG channel, we artificially modified channel TP9 to illustrate BrainBeats' bad channel detection and removal algorithm. Similarly, we artificially simulate a large EEG artifact at the beginning of the file and a high-frequency muscle-like artifact in the temporal channels from 3-6 s to illustrate BrainBeats' artifact removal feature.

2. Heartbeat-evoked potentials (HEP) and oscillations (HEO)

1. Open BrainBeats' first general user interface (GUI) to select the main parameters. In EEGLAB, click Tools > BrainBeats > 1^st level (subject level). Select Heartbeat-evoked potentials HEP as the analysis to run, ECG as the heart data type, and click on the button to display the list of channels to select the ECG channel labeled ECG (or type it directly in the text box). Keep the options Visualize outputs (recommended) and Save outputs selected and click Ok.
2. A second GUI window pops up based on the previous choices (i.e., HEP analysis with ECG data). Select preprocessing parameters for both EEG and ECG signals. Turn off preprocessing by unchecking the boxes Preprocess ECG and Preprocess EEG to see if the data has already been preprocessed prior to launching BrainBeats. Change Power line noise à 50 Hz in the Preprocess EEG section since these EEG data were recorded in Russia. Click OK to launch.
3. A warning message will appear, asking for confirmation to remove the PPG channel that has been detected. This is because the toolbox is not designed to analyze both ECG and PPG simultaneously (or other auxiliary channels) and keeping it in the dataset will lead to severe errors (e.g., artifact removal, poor ICA decomposition, etc.). Click Yes. BrainBeats starts performing some checks, setting some default parameters, and separates the ECG from the EEG data to preprocess the ECG signal and compute the RR intervals.
4. ECG and RR time series are preprocessed using validated algorithms from the Physionet Cardiovascular Signal toolbox^46,64 adapted to fit BrainBeats' data formatting, increase clarity, parameter tuning, and computing time (see references for validation of the algorithms). The plugin outputs the RR intervals, timestamps, filtered ECG signal, R-peaks indices, and heart rate (HR). Tune these parameters via the GUI or command line.
   NOTE: The ECG signal is bandpass-filtered using a customized, validated filter (1-30 Hz) and scans the signal to identify the QRS complex and R-peaks using the Pan-Tompkins (P&T) method⁶⁵, implementing some signal processing operations including differentiation, squaring, integration, and smoothing for best performance. The P&T energy threshold is estimated based on the sample rate and the smoothed ECG values to avoid disruption from large bumps. If the RR interval variability exceeds 1.5 times the median, it searches for missed peaks. The mean R-peak sign is calculated over 30 s segments, and peak points are refined through a refractory period check, managing flatline conditions and ensuring consistent detection.
5. Next, BrainBeats identifies abnormal RR intervals or spikes within RR intervals using a forward-backward search and physiological thresholds. The signal quality index (SQI) is calculated⁴⁶; check this if the system displays warnings where more than 20% of the RR times series contains RR artifacts (outside of physiological limits or with an SQI below 0.9). A plot displays the filtered ECG signal, identified R-peaks, NN intervals, and interpolated artifacts (see Figure 2).
   NOTE: The RR artifacts are interpolated by default using the shape-preserving piecewise cubic method to obtain the normal-to-normal (NN) intervals but remove them (not recommended) or use another interpolation method (linear, cubic, nearest neighbor, previous/next neighbor, spline, cubic convolution, or modified Akima cubic) if needed. When several ECG channels are present, the RR intervals are estimated for each, and the channel with the least number of RR artifacts is selected for the following steps.
6. Scroll through zoomed-in 30 s windows of the R-peaks for closer inspection by pressing the Right/left Arrows. If the data contains several ECG/PPG channels, use the channel with the best signal quality index for the RR intervals. BrainBeats does not support both ECG and PPG signals simultaneously at the moment. For tutorial purposes, we chose a sample dataset that contains both datatypes.
7. Once done with the ECG signal, BrainBeats bandpass filters the EEG data at 1-40 Hz using a nonlinear causal minimum-phase FIR filter by default to reduce smearing activity between pre- and post-heartbeat periods, preserve causality, and avoid undesired group delays⁴⁴. This is particularly important for users examining the pre-heartbeat period. If the lowpass filter is set to a value above the power line frequency (e.g., 80 Hz lowpass with power line frequency at 50 Hz), use a sharp notch filter to remove the line noise artifact. The EEG data are then re-referenced to infinity using the REST algorithm (best suited for HEP analysis⁴⁹) unless less than 30 channels are detected (in which case they cannot be reliably re-referenced, and a warning is generated to let users know).
8. BrainBeats then detects, removes, and interpolates abnormal EEG channels (Figure 3). Check that the default parameters are flat lines greater than 5 s (clean_flatlines algorithm), a maximum high-frequency noise standard deviation of 10, a window length of 5 s (to better capture slow-frequency artifacts⁵²), a minimum correlation between neighboring channels of .65, and a maximum tolerated portion of 33% (clean_channels algorithm).
   NOTE: The number of RANSAC samples is set to 500 by default to increase convergence and replicability of bad channel rejection (although it increases the computation time).
9. Next, insert R-peaks as event markers into the EEG data to mark each heartbeat and segment the data around these markers without baseline removal (per guidelines^43,66; Figure 4). Since NN intervals have different lengths and EEG must be segmented at a constant length, estimate the minimum epoch size cutoff following R-peak events using the 5^th percentile of the interbeat-interval (IBI) data (i.e., the value below which 5% of the shortest IBIs fall, displayed as a dashed red line on a histogram; see Figure 5).
   NOTE: This 5^th percentile value is a good compromise to preserve as many epochs as possible while ensuring they are not too short since the period of interest for HEP/HEO analysis is 200-600 ms post heartbeat^49,67.
10. Segment EEG data from -300 ms before the R-peaks to the 5^th percentile value post-R-peak, with the R-peak at time 0. Reject epochs shorter than 550 ms or containing more than one R-peak (which would bias the ERP/ERSP), per guidelines^49,67. Epochs containing large EEG artifacts are detected using root-mean-square (RMS) and signal-to-noise ratio (SNR) metrics and MATLAB's isoutlier function (Figure 6). Artifactual epochs are removed.
11. Perform blind source separation using the default Infomax algorithm, accounting for data rank for best performance⁵³. If desired, choose the preconditioned independent component analysis (PICARD) algorithm for fast computation^54,54 by choosing the option Picard (fast) or the modified Infomax algorithm by choosing the option Modified Infomax (long but replicable) for the field ICA method (although this involves much longer computation times). Use the ICLabel plugin⁵⁶ to automatically classify ICs to extract non-brain artifacts (ocular components are removed with 90% confidence, whereas muscular, line noise and channel noise are removed with 99% confidence; Figure 7).
12. Keep the option Visualize outputs selected in the first GUI window to visualize the final EEG time series, the grand average HEP (Figure 8 top), HEP for each heartbeat (Figure 8 bottom), and the grand average heartbeat-evoked oscillations (HEO; Figure 9). Examine HEOs in terms of event-related spectral perturbation (ERSP, i.e., changes in EEG power across heartbeats; Figure 9 top) and inter-trial phase coherence (ITC, i.e., consistency of the phase angle across heartbeats; Figure 9 bottom).
    NOTE: ERSP is computed using a default 3-cycle wavelet (with a Hanning-tapered window applied, pad ratio of 2) and with the number of cycles in the wavelets used for higher frequencies expanding slowly up to 20% of the number of cycles in the equivalent FFT window at its highest frequency (1 minus 0.8). This controls the shapes of the individual windows measured by the function and their shapes in the resulting time/frequency panes. An arbitrary baseline is removed for illustration purposes, and ERSP is computed for frequencies 7-25 Hz to capture the typical HEO effect described in the literature, namely 300-450 ms post heartbeat in the alpha band (8-13 Hz) over frontocentral electrodes^17,67. Lower frequencies cannot be estimated due to the short epoch size defined by the interbeat intervals. Nonparametric (permutation) statistics are applied to visualize the HEO for a p-value of 0.05, corrected for false discovery rate (FDR, i.e., type 1 error or family-wise error). These plots are generated for tutorial purposes or single-trials analysis.
13. Preprocessing plots are generated to visualize the different steps. To turn it off, uncheck the box Visualize preprocessings. The final EEG data (cleaned and segmented around the R-peaks) do not include the ECG data since it would bias the ERP/ERSP analysis. To preserve the heart channel in the final output, check the box Keep heart channel.
    NOTE: Pause here if needed before processing the next file (next condition or participant).
14. BrainBeats supports EEGLAB's history function. At the end of all operations, type eegh in MATLAB's command window to print the command line to repeat all the above steps via a single command line, with the parameters that were selected manually in the GUI, allowing easy automation. Find preprocessing outputs (e.g., signal quality index of the cardiovascular time series, NN intervals, removed EEG channels, segments, and components, etc.) in the EEGLAB structure: EEG.brainbeats.preprocessings. All parameters are also exported in EEG.brainbeats.parameters.
15. For advanced users, perform all the above steps with default parameters with the following command lines (see the tutorial script in the BrainBeats repository for more options):
    eeglab; close; % Launch EEGLAB without the GUI
    main_path = fileparts(which('eegplugin_BrainBeats.m')); cd(main_path);
    EEG = pop_loadset('filename','dataset.set',
    'filepath',fullfile(main_path, 'sample_data')); %Load the sample dataset
    EEG = brainbeats_process(EEG,'analysis',
    'hep','heart_signal','ECG', 'heart_channels',
    {'ECG'},'clean_eeg',true); % Launch BrainBeats 1^st level to process the file for HEP analysis with default parameters
16. The steps above performed HEP/HEO from the ECG signal. Use the following steps for the PPG signal.
17. In the following steps, perform the same operations but using a PPG signal. Load the same dataset again (see step 1.5) since the previous operations overwrote it, and open BrainBeats' first GUI again to select the main parameters. Click Tools > BrainBeats > 1^st level (subject level). Select Heartbeat-evoked potentials (HEP) as the analysis to run, PPG as the heart data type, and click on the button to display the list of channels to select the PPG channel. Keep the options Visualize outputs (recommended) and Save outputs selected and click OK.
18. The second GUI window pops up in a manner similar to step 2.2 . The only difference is that we can now fine-tune the parameter to process the PPG signal instead of the ECG signal. Click OK to run with the default parameters.
19. A warning message appears, asking for confirmation to remove the extra ECG channel that has been detected. Again, this is expected. Click Yes. By default, the toolbox will preprocess the PPG signal, detect the pulse waves to obtain the RR intervals, identify RR artifacts, if any, and interpolate them (Figure 10). Steps 2.7. to 2.12. are performed, and the same plots and outputs are generated but based on the R waves detected from the PPG signal (see Figure 11 and Figure 12).
    NOTE: R waves are detected using the slope of the signal within a specified window. Potential pulses are then flagged when the slope exceeds a dynamic threshold, which is adjusted based on the detection history and signal characteristics. The algorithm then searches within an eye-closing period to pinpoint the maximum slope, and subsequently, the onset of the pulse wave is determined through thresholding. The R-wave peaks are identified as the valleys near the onset, and their locations are recorded. The algorithm iterates through the entire signal, continuously adjusting detection thresholds and identifying R-wave peaks, which are then used to calculate the RR intervals.
20. For advanced users, perform all the above steps with default parameters with the following command lines (see the tutorial script in the BrainBeats repository for more options):
    eeglab; close; % Launch EEGLAB without the GUI
    main_path = fileparts(which('eegplugin_
    BrainBeats.m')); cd(main_path);
    EEG = pop_loadset('filename','dataset.set',
    'filepath',fullfile(main_path, 'sample_data')); %Load the sample dataset
    EEG = brainbeats_process(EEG,'analysis',
    'hep','heart_signal','PPG', 'heart_channels',
    {'PPG'},'clean_eeg',true); % Launch BrainBeats 1^st level to process the file for HEP analysis with default parameters

3. Extracting EEG and HRV features

1. Load the same dataset again (see step 1.5; Click File > Load existing dataset > Select dataset.set) since it was overwritten by the previous operations and open the main GUI again to select the main parameters (step 2.1; Click Tools > BrainBeats > 1^st level). Select Extract EEG & HRV features for the analysis type, ECG for the heart signal type, and select ECG in the list of electrode labels. Click OK.
2. The second GUI window pops up like in step 2.2 but with different parameters for EEG preprocessing and extracting HRV and EEG features. Set these parameters as described below. Click OK to launch with default parameters.
   1. In the HRV section, click on the button freq. options to select the method to compute HRV power (default set to normalized Lomb-Scargle periodogram), the window overlap (default set to 25%), and for performing a second-level normalization (unset by default; see note below for more detail).
   2. In the EEG features section, click on the button freq. options to fine-tune some parameters, such as the overall frequency range on which to compute the power spectral density (PSD; default = 1-40 Hz), the units (decibels, µV²/Hz, or normalized by total power), the window type (default = hamming), the window overlap (default = 50%), the window length (default = 2 s), and the types of frequency bounds for each band.
      NOTE: HRV power is computed by default using the normalized Lomb-Scargle periodogram, which does not require resampling (therefore better preserving original information) and best deals with non-uniformly sampled data, missing data, and noise (typical with NN intervals)⁶⁸. The normalized version scales the power by the variance of the signal, providing results that are less sensitive to varying noise levels, more focused on the relative strength of periodic components, and more comparable across different recordings or subjects. Other methods available include the non-normalized Lomb-Scargle periodogram, the Welch method, and the Fast Fourier Transform (FFT). Resampling is automatically performed for the Welch and FFT methods to create the necessary regularly sampled time series. A second-level normalization can be applied by dividing the power of each frequency band by the total power, providing a more intuitive measure of the relative contribution of each frequency component to the overall power. It is disabled by default since it is only meaningful when all four bands are available, requiring at least 24 h of signal. These algorithms are adapted from the Physionet Cardiovascular signal processing toolbox⁴⁶.
3. The same warning message appears, asking for confirmation to remove the extra PPG channel that has been detected. Again, this is expected. Click Yes. BrainBeats will start preprocessing the ECG data and extracting the NN intervals identically, as in step 2.4. Then, it extracts heart-rate variability (HRV) features from the NN intervals in the time (SDNN, RMSSD, pNN50), frequency (ULF, VLF, LF, HF, LF:HF ratio, total power), and nonlinear (Poincare, phase-rectified signal averaging, fuzzy entropy, and fractal dimension) domains.
   NOTE: BrainBeats automatically checks the file length to ensure that minimum requirements are met (e.g., ULF-HRV power requires 24 h of data), sends warning messages if not, and does not export these features to prevent unreliable estimations. BrainBeats follows guidelines and recommendations for estimating HRV metrics^27,45.
4. BrainBeats preprocesses the EEG data as in step 2.7 . (see Figure 3). Large artifacts are detected automatically in the continuous data using the artifact subspace reconstruction (ASR) algorithm^50,52 (default SD criterion set to 30 and using 80% of available RAM to increase speed). These large artifacts are removed from the EEG data (see Figure 13). To adjust these parameters in the GUI, select the fields Threshold to reject bad segments with ASR and Available RAM to use for ASR.
   NOTE: The EEG and cardiovascular time series do not need to be time synchronized for the features mode since the features are estimated on each signal separately. Thus, EEG artifacts can be removed directly from the EEG data (in red, Figure 13), unlike for the HEP mode, where epochs containing artifacts were rejected for both time series because time synchronization with ms accuracy is essential for that method.
5. Perform ICA using the same algorithms and parameters as for HEP (see step 2.11.), except that this time the heart components are removed if detected with 99% confidence (they were preserved for HEP/HEO since we do not want to remove relevant cardiac-related signals).
6. The box Frequency domain is checked in step 3.2., so BrainBeats extracts the following frequency domain features: the average power spectral density (PSD) for the delta (1-3 Hz), theta (3-7 Hz), alpha (8-13 Hz), beta (13-30 Hz), and gamma (30+ Hz) frequency bands, the individual alpha frequency (IAF), and alpha asymmetry on all available (symmetric) electrode pairs. Check these features.
   NOTE: PSD converted to decibels (dB) facilitates the comparison of results across recordings and subjects. The frequency bounds can be set to the conventional bounds (e.g., predefined 8-13 Hz for the alpha band) or to the individualized bounds, which are detected from the distribution of the power spectral density to account for inter-individual differences⁶⁹ (e.g., 7.3-12.6 Hz for the alpha band). The algorithm was designed for the alpha band and does not perform as well for other bands, especially when peaks are not present in the power spectral distribution. The individual alpha frequency (IAF) is estimated using the alpha center of gravity to better deal with split peaks or ambiguous peaks⁶⁹. Alpha asymmetry is computed following guidelines (2-s hamming window with 50% overlap, the logarithm of alpha power from the left channel minus logarithm of alpha power from the right channel)⁴⁷. Hence, positive values indicate greater left-than-right alpha power and vice versa. Alpha asymmetry can be normalized by dividing alpha power from each electrode by the alpha power summed across all electrodes⁴⁷. The symmetric pairs are obtained using theta distances, requiring the EEG data to contain electrode coordinates.
7. The box Nonlinear domain is checked in step 3.2., BrainBeats extracts the fuzzy entropy and the fractal dimension for each EEG channel.
   NOTE: Nonlinear-domain features are thought to capture nonlinear, complex dynamics of the brain that are missed by spectral measures and show particular promise for investigating interactions between various body systems^70,71,72. Fuzzy entropy is more reliable and robust than its alternatives (sample and approximate entropy) but requires longer computation times (especially with long EEG time series with high sampling rates). To address this issue, when EEG signals are longer than 2 minutes long with a sampling rate greater than 100 Hz, they are automatically downsampled (or decimated when the factor is not an integer) to 90 Hz (i.e., corresponding to a Nyquist frequency of 45 Hz, to match the default lowpass filter and avoid line noise artifacts as much as possible). Furthermore, parallel computing is activated by default when estimating EEG features, which reduces computing time, especially when many EEG channels are available.
8. Select Visualize outputs in the first GUI (see step 3.1.), to generate a plot displaying the power spectral density (PSD) for HRV and EEG data (Figure 14) along with scalp topographies displaying some EEG features (Figure 15).
   NOTE: You can also find some preprocessing outputs in EEG.brainbeats.preprocessing and all the parameters used in EEG.brainbeats.parameters. We encourage users to report these outputs in scientific publications to increase the replicability of findings.
9. Keep the Save outputs box checked in the first GUI window to save all features exported into the EEGLAB .set file in EEG.brainbeats.features and save them in a .mat file in the same folder where the dataset was loaded.
10. BrainBeats supports EEGLAB's history function. At the end of all operations, type eegh in MATLAB's command window to print the command linethat will allow you to repeat all the above steps via a single command line, with the parameters that were selected manually in the GUI, allowing easy automation and replication of the operations.
11. For advanced users, perform all the steps above with the following command:
    eeglab; close; % Launch EEGLAB without the GUI
    main_path = fileparts(which('eegplugin_BrainBeats.m')); cd(main_path);
    EEG = pop_loadset('filename','dataset.set',
    'filepath',fullfile(main_path, 'sample_data')); % Load the sample dataset
    EEG = brainbeats_process(EEG,' analysis',
    'features','heart_signal','ECG', 'heart_channels',{'ECG'},'clean_eeg',true);
12. The previous steps extracted HRV features from ECG signal. Use the following steps to extract HRV features from PPG signal (EEG features are the same).
13. Load the same dataset again (step 1.5.) since it was overwritten by the operations and open the main GUI again (step 2.5.). Select Extract EEG & HRV features for the analysis and select PPG for the heart signal type and PPG for the channel name. Click OK.
14. The 2^nd GUI window now shows the parameters for preprocessing PPG and for extracting the HRV and EEG features. Click OK to run with the default parameters. Parameters are described in step 2.17.
15. A warning message will appear, asking for confirmation to remove the detected ECG channel. This is to be expected since the toolbox is not designed to analyze both ECG and PPG simultaneously (or other auxiliary channels) and keeping it in the dataset will lead to serious errors (e.g., artifact removal, poor ICA decomposition, etc.). Click Yes.
16. BrainBeats preprocesses the PPG signal and estimates the NN intervals as in step 2.5. Check the extracted HRV features from the NN intervals, the same as in step 3.2 . except that the NN intervals have now been obtained from the PPG signal. The EEG signals are preprocessed as in step 3.2. BrainBeats plots the PSD (Figure 16). The only difference here is the PSD estimated from the NN intervals obtained from the PPG as opposed to ECG.
17. For advanced users, perform all the steps above with the following command:
    eeglab; close; % Launch EEGLAB without the GUI
    main_path = fileparts(which('eegplugin_BrainBeats.m')); cd(main_path);
    EEG = pop_loadset('filename','dataset.set','filepath',
    fullfile(main_path, 'sample_data')); %Load the sample dataset
    EEG = brainbeats_process(EEG,'analysis',
    'features','heart_signal', 'PPG','heart_channels',{'PPG'},'clean_eeg',true);

4. Extract heart artifacts from EEG signals.

1. Load the sample dataset (see step 1.5.).
2. Open the main GUI window by clicking Tools > BrainBeats > 1^st level (subject level) and select Extract heart artifacts from EEG signals for the analysis type, ECG for heart signal type and select ECG in the list of electrode labels. Click OK.
3. The 2^nd GUI window shows the preprocessing parameters. Set the Power line noise à 50 Hz (Europe), edit the confidence level if needed, check the box Boost mode (beta), and click OK to run with the default parameters since the EEG signals from the sample dataset are not preprocessed.
   NOTE: the confidence level to detect heart components is set to 80% by default, which may be too low or too high for some datasets. Increasing this value will increase the chances of detecting heart components but decrease the reliability of that detection. The Boost mode (beta) is optional and aims to improve classification performance by smearing the heart signal into the EEG signals.
4. A warning message will appear, asking for confirmation to remove the extra PPG channel that has been detected. This is to be expected since the toolbox is not designed to analyze both ECG and PPG simultaneously (or other auxiliary channels) and keeping it in the dataset will lead to serious errors (e.g., artifact removal, poor ICA decomposition, etc.). Click Yes.
5. The ECG signal is bandpass-filtered to remove slow frequency drifts below 1 Hz and high-frequency noise above 20 Hz (with a noncausal zero-phase FIR filter). Preprocess the EEG signals as in step 3.4.
6. Perform independent component analysis (ICA) using the preconditioned ICA for real data algorithm (PICARD). To change this option, choose the standard Infomax algorithm or modified Infomax algorithm for replication from the GUI in step 4.3 . Then, classify the independent components automatically with ICLabel. If a component is classified as a heart component with 80% confidence, it is, by default, detected and extracted automatically from the EEG data.
7. Keep the box Visualize outputs in the first main GUI (step 4.2.) to visualize the scalp topography of the removed component (Figure 17 Left) and the final EEG time series (in blue, Figure 17 Right) after extracting the heart component (in red, Figure 17 Right).
   NOTE: the ECG channel is preserved for visualization to confirm the extraction of ECG-related components, but it is removed after this step since it no longer contains any relevant information.
8. For advanced users, perform these steps using the following command lines:
   eeglab; close; % Launch EEGLAB without the GUI
   main_path = fileparts(which('eegplugin_BrainBeats.m')); cd(main_path);
   EEG = pop_loadset('filename','dataset.set','filepath',
   fullfile(main_path, 'sample_data'));
   EEG = brainbeats_process(EEG,'analysis','rm_heart',
   'heart_signal','ECG', …
   'heart_channels',{'ECG'},'clean_eeg',true,'vis_cleaning',false,…
   'conf_thresh',.8,'boost',true);