EEG-Beats: Automated analysis of heart rate variability (HVR) from EEG-EKG

Heart rate variability (HRV), the variation of the period between consecutive heartbeats, is an established tool for assessing physiological indicators such as stress and fatigue. In non-clinical settings, HRV is often computed from signals acquired using wearable devices that are susceptible to strong artifacts. In EEG (electroencephalography) experiments, these devices must be synchronized with the EEG and typically provide intermittent interbeat interval information based on proprietary artifact-removal algorithms. This paper describes an automated algorithm that uses the output of an EEG sensor mounted on a subject’s chest to accurately detect interbeat intervals and to calculate time-varying metrics. The algorithm is designed for raw signals and is robust to artifacts, resulting in fine-grained capture of HRV that is synchronized with the EEG. An open-source MATLAB toolbox (EEG-Beats) is available to calculate interbeat intervals and many standard HRV time and frequency indicators. EEG-Beats is designed to run in a completely automated fashion on an entire study without manual intervention. The paper applies EEG-Beats to EKG signals measured with an EEG sensor in a large longitudinal study (17 subjects, 6 tasks, 854 datasets). The toolbox is available at https://github.com/VisLab/EEG-Beats.


Introduction 23
Heart rate variability is an indicator of many physiological and behavioral factors (Acharya et al., 2006) 24 such as stress (Kim et al., 2018), fatigue (Vicente et al., 2016), and performance (Spangler et al., 2020) 25 in normal subjects (Nunan et al., 2010) (Shaffer & Ginsberg, 2017. HRV measures can be extracted 26 from signals generated from a variety of sensor types including electrocardiogram (EKG), 27 electroencephalography (EEG), photoplethysmogram (PPG), and blood pressure monitors (ABP). 28 HRV measures have been standardized (Task Force of the European Society of Cardiology and the 29 North American Society of Pacing and Electrophysiology, 1996). Many open-source and proprietary 30 tools have been developed for detecting peaks in cardiac signals and computing these measures. 31 The benchmark study by Vest et al. (Vest et al., 2018) compares five HRV toolboxes: PhysioNet HRV 32 Toolkit (Goldberger et al., 2000), Kubios (Tarvainen et al., 2014), the Kaplan toolbox (Kaplan & 33 Staffin, 1998) parameter settings and caution against using default parameters when analyzing raw EKG. Areas of 37 particular concern include preprocessing, signal quality assessment for noisy segment removal, and 38 detection of arrhythmias. 39 Traditionally, the focus of EKG in EEG experiments has been the removal of cardiac interference 40 (Tamburro et al., 2019). Our interest in HRV arose from the possibility of using HRV indicators 41 secondary measures of subject physiological state during EEG experiments. By placing a single EEG 42 sensor on the chest, it is possible to extract HRV as a matter of routine during EEG data analysis. 43 However, our initial efforts to analyze a large EEG/EKG study using several available toolboxes 44 including Kubios and the PhysioNet Cardiovascular Signal Toolbox (PNC) was unsuccessful because 45 peak detection often failed without extensive manual intervention. 46 In addition to large variations in both baseline signal levels and peak amplitudes, EEG signals also tend 47 to have large signal bursts due to muscle activity and loose detectors. To address the difficulty of 48 consistent large-scale extraction of HRV measures from EEG, the approach proposed in this paper uses 49 a top-down, divide-and-conquer strategy to peak-finding in contrast to typical sliding window 50 approaches such as Pan-Tompkins (Pan & Tompkins, 1985). While straightforward and robust to 51 typical EEG artifacts, but the extraction method is applicable only to normal cardiac signals. The 52 method detects positions of the peaks corresponding to heart beats and but not QRS complexes or other 53 feature information such as arterial fibrillation or widespread arrhythmias. 54 This paper describes the HRV algorithm and an open-source implementation in a MATLAB toolbox. 55 The toolbox, EEG-Beats, extracts inter-beat intervals with associated HRV measures and includes 56 utilities for statistical analysis and visualization. EEG-Beats can be run in automated fashion on an 57 entire EEG study. The results are validated by comparison with PNC, and examples of how the toolbox 58 might be used are provided. 59

Methods and Materials 60
An RR interval is the distance between successive peaks without regard to the shape of the associated 61 QRS complexes. Normal RR intervals, referred to in the literature as NN intervals, are not distinguished 62 from other RR intervals in this paper. This section describes the automated EEG-Beats processing 63 algorithms for extracting interbeat (RR) intervals, calculating HRV indicators from RR intervals, and 64 analyzing the relationships of the extracted indicators. All explanations,computations,and examples 65 in this paper use the EEG-Beats default parameter values. When values are given, the names of the 66 corresponding settable parameters in parentheses appear in parentheses. EEG-Beats uses the statistics 67 of an entire signal to determine peaks and should not be applied in an online fashion. Signals are 68 recommended to be at least 5 minutes in duration. 69 In this paper, the central tendency is generally indicated by the median and distribution characteristics 70 by the interquartile range (IQR). Thresholds are based on the robust standard deviation, defined as 71 1.4826 times the median absolute deviation from the median. The robust outlier criteria specifies that 72 points that are 1.5 (maxWhisker) × IQR outside the first and third quartiles (mid quartiles) of a 73 distribution are considered to be outliers. 74

Signal preparation 75
EEG-Beats downsamples the signal to 128 Hz (srate) and applies a [3,20] Hz (filterHz) finite impulse 76 response (FIR) filter to the EKG signal. The 3 Hz high-pass filter effectively removes trend, while the 77 20 Hz low-pass filter removes much of the high-frequency noise. The algorithm then subtracts the 78 median signal and truncates the amplitude so that the signal is within 15 (truncateThreshold) robust 79 standard deviations of the median signal. Truncation eliminates extreme signal spikes, which in the 80 case of a loose detector can be thousands of times larger than the normal signal for EEG sensors. 81

Heartbeat detection 82
EKG acquired using EEG sensors has two distinctive characteristics that make setting of thresholds 83 problematic for traditional heartbeat algorithms. First, the signal can undergo spikes in amplitude 84 thousands of times larger than the base signal when a sensor makes poor contact. Second, the amplitude 85 and signal directions can vary over time, as illustrated by the signal clip shown in Fig. 1

93
EEG-Beats determines the peak locations of the heart beats by first orienting the signal so that the 94 peaks are always in the positive direction and troughs (if present) are to the right of their associated 95 peaks. Reorientation simplifies the subsequent steps in the algorithm. The algorithm then uses a process 96 of successive refinement to identify the largest peaks and then refine the partition to include smaller 97 peaks. Finally, there is a clean-up phase that removes extraneous peaks. The following subsections 98 describe these steps in more detail. 99

Determining a consensus direction 100
Because large artifacts can obscure the true direction of the heartbeats, EEG-Beats uses a consensus 101 algorithm to determine how to reorient the signal if necessary. By default, the algorithm divides up the 102 signal into two-second (flipIntervalSeconds) intervals and finds the direction of the maximum absolute 103 value in each interval. Only intervals whose most extreme value is more than 1.5 (threshold) robust 104 standard deviations from the overall median are considered when determining whether to negate the 105 signal, based on the dominant direction of the extreme points. A similar algorithm determines whether, 106 assuming that a peak is associated with a trough, the trough is before or after the peak. If necessary, 107 EEG-Beats flips/reverses the signal so that most of the peaks are up with following troughs. 108

Selecting eligible peaks 109
To determine whether a particular point could really be a heartbeat, the algorithm assumes that eligible 110 points must have a maximum amplitude of at least 1.5 (threshold) robust standard deviations from the 111 median and must be at least 500 (rrMinMs) ms away from points that have already been determined to 112 be heartbeats. The rrMinMs and rrMaxMs parameters determine the allowed range of the RR interval 113 sizes. By default these are set to 500 ms and 1500 ms, corresponding to heart rates of 120 bpm and 40 114 bpm, respectively. While appropriate for normal EEG recordings, these settings should be adjusted for 115 recordings during which the subject is performing strenuous or stressful activities. 116 In addition to position and amplitude criteria, EEG-Beats uses a sharpness to determine whether a point 117 is an eligible peak. In the case of single peak, the signal at an eligible peak must fall below the median 118 signal within 100 (0.5*qrsDuration) ms on either side of the peak value. In the case of a peak followed 119 by a trough, an eligible peak's trough (as determined by MATLAB's findpeaks function) must be 120 within 100 (0.5*qrsDuration) ms of the peak. These checks determine whether a large deviation in 121 signal amplitude actually corresponds to a sharp peak. EEG-Beats applies the algorithm twice (once 122 assuming all peaks are single and once assuming all peaks have following troughs) and then combines 123 the results as described below. 124

Successively subdividing to find the peaks 125
The main algorithm uses a divide-and-conquer strategy similar to that used in binary search.  Beats divides the signal into 31 (consensusIntervals) equal-size intervals and finds the maximum point 127 (the fence post) in each interval. After adding the first point and the last point in the signal as outer 128 fence posts, EEG-Beats eliminates internal fence posts that are not eligible peaks before beginning 129 divide and conquer. 130 The divide-and-conquer phase proceeds as follows. For each pair of consecutive fence posts, the 131 algorithm finds the maximum between the first two fence posts in the list. If this maximum point is an 132 eligible peak, the algorithm inserts the point as a fence post between the original two fence posts. 133 Regardless of whether the point is an eligible peak, the algorithm zeros out the signal within a 200 134 (qrsDuration) ms window around this point and continues the process. When there are no eligible 135 peaks between the first two fence posts, the algorithm removes the first fence post and repeats the 136 process with the next pair, until no fence posts pairs remain to be processed. 137

Cleaning up 138
Since the signal was truncated, the actual peaks may be slightly off. EEG-Beats makes a pass through 139 the peak list to adjust the actual peak positions to their true maximum positions and then combines the 140 peaks from the two methods (single-peak versus peak-trough). Most peaks from the two methods will 141 be coincident. For unmatched peaks that are within 100 (0.5*qrsDuration) ms of each other, EEG-142 Beats adjusts the position of the peak with the smaller amplitude to be that of the peak with the larger. 143 A second clean-up step is to remove extraneous peaks by determining whether removing the peak 144 produces a more consistent representation. For that to be the case, a peak's neighbors must be within 145 1500 (rrMaxMs) ms and the interbeat interval with peak removed is closer to the median interbeat 146 interval than the largest of the interbeat intervals on either side of that peak. EEG-Beats also marks 147 outlier peaks whose absolute amplitude value is either less than 0.5 (minPeakAmpRatio) or greater than 148 2 (maxPeakAmpRatio) times the absolute value of the median peak amplitude. 149

Combining peaks 150
The final stage uses successive combination to get a single representation of heartbeat. EEG-Beats sets 151 the group (single-peak or peak-trough) with the greatest number of peaks as the base representation 152 and considers the remaining peaks one at a time starting from left to right. EEG-Beats adds the peaks 153 to the base representation as long as they aren't closer than 500 (rrMinMs) ms to peaks in the existing 154 representation. Peaks before the first peak and after the last peak in the base representation are handled 155 slightly differently, with new peaks added from inside outward. EEG-Beats also reports, but does not 156 remove peaks whose amplitudes are likely to be too high or too low to be beats using the standard 157 robust outlier criteria for boxplots. Peaks whose amplitude satisfies the robust outlier criteria are 158 reported as high or low amplitude peaks. 159

Saving the results 160
EEG-beats has an EEGLAB plug-in that allows users to run the analysis on single recordings. 161 However, EEG-Beats is meant to be run using scripts on an entire study in a directory tree. Example 162 scripts show how to provide the source root directory for a study and the root directory for saving the 163 results. The script automatically processes all of the EEG .set files in the source directory tree and 164 produces plots for each recording similar to those shown in Fig. 1. EEG-Beats also saves a structure 165 (ekgPeaks) containing the peaks found for each recording. This structure can be used as input to 166 produce and save a structure containing the RR measures (rrInfo) calculated from the interbeat intervals 167 for the entire study. EEG-Beats also provides scripts to do automated analysis of variance and 168 visualizations of RR measure distributions if the user provides a suitably formatted metadata file as 169 described below. 170

Available RR measures for HRV 171
In clinical settings, HRV is often analyzed for short-term variations in time windows ranging from 30 172 seconds to 5 minutes.

Analysis and visualization 201
In addition the visualization of the peak locations and distribution of peaks versus RR intervals (e.g. 202 Fig. 1), EEG-Beats provides a visualization of the RR intervals as a function of time overlaid on (e.g., 203 Fig. 2) or as a subplot with the EKG signal and peaks (Fig. 4). RR intervals detected as outliers using 204 various strategies are marked on these RR interval graphs using different symbols. 205 EEG-Beats also has several summary visualizations and statistical analyses designed to be run on a 206 study as a whole. The user must provide metadata in a MATLAB structure with one row for each 207 dataset. The structure must have a fileName field containing the name of the file that the EKG was 208 extracted from. That fileName is matched with the fileName field in the ekgPeaks structure to perform 209 analysis. 210 EEG-Beats can compare peak representations beat-by-beat for two different toolboxes (e.g., Table 3)  211 or for output from EEG-Beats for different peak-finding parameter settings. The EEG-Beats RR 212 measure compare functions allow comparison of results from different packages or from using different 213 EEG-Beats artifact removal settings (e.g., Table 4). 214 EEG-Beats allows users to form feature vectors representing HRV in each segment (5 minutes by 215 default) by selecting a list of the RR measures from those available from Table 1. EEG-beats  216 normalizes the features by z-scoring and uses MATLAB's tsne function (t-distributed stochastic 217 neighbor embedding) to project these high-dimensional feature vectors into a 2D or 3D space (Van Der  218 Maaten & Hinton, 2008). The idea behind t-SNE is that high-dimensional vectors projecting to the 219 same low-dimensional clusters are likely to also be similar in the original high-dimensional space. 220 EEG-Beats provides the ability to plot these projections using colors and/or shapes to distinguish 221 different metadata values. 222 EEG-Beats can also plot the boxplots of any RR measure of Table 1 segregated by any of the metadata  223 variables (e.g., Fig. 6). Finally EEG-Beats can apply an analysis of variance using MATLAB's anovan 224 function for any RR measure using factors corresponding to any combination of fields in the user-225 provided metadata structure (e.g., Tables 5 and 6). 226

Data used for testing 227
This study uses the NCTU RWN_VDE data collection from a longitudinal experiment conducted at 228 National subjects wore an Actigraph activity monitor for the duration of the study and completed several 231 subjective assessments of stress, fatigue, and sleep quality on a daily basis. Central to the study was 232 the use of the Daily Sampling System (DSS), which processed this information and automatically 233 invited subjects to the laboratory to undergo the experiment based on putative subject state. 234 The 17 subjects were recorded on up to 9 days, ideally representing low, normal, and high fatigue

Manual verification of heart beat locations 246
To validate EEG-Beats peak-finding, we manually reviewed the EEG-Beats peak positions for all 854 247 EEG data recordings of the NCTU RWN_VDE study using the pan and zoom features for MATLAB 248 figures with the peak and point plots of Fig. 1. A summary is shown in Table 2. 249 Total peaks indicates the total number of peaks detected in the data corpus. EEG-Beats accurately 251 detected peaks using the default parameter settings in a large portion (626) of the recordings, many of 252 which had substantial artifacts including dramatic variations in amplitude, dropouts, and high-253 amplitude bursts. Peaks identified as having too high or too low of an amplitude use the standard robust 254 outlier criteria described in the methods. Another 181 datasets had a few errors, while 34 datasets had 255 a "larger" number of errors. Fig. 2 uses an overlay plot to display Dataset 800, an example of a dataset 256 with a "larger" number of errors. 257

264
Dataset 800 has a very noisy signal with a number of very low and very high peak amplitudes as 265 indicated by the large number of points on the vertical axes in the right graph of Fig. 2. Although EEG-266 Beats misidentifies some peak positions for this dataset, all of the misidentifications occur for peaks 267 on one side or the other RR intervals marked either as bad neighbor RRs or as outlier peak RRs on the 268 left graph. Using the three RR outlier removal criteria removes most of the RR outliers in datasets with 269 peak errors. 270 Most of EEG-Beats misidentifications occur when a peak is overlaid with or is located close to an 271 artifact as illustrated by the example in Fig. 3. The vertical arrow at 620 seconds in the left graph of 272 Fig. 3 marks a large-amplitude artifact that EEG-Beats labeled as a peak. Once that peak was marked 273 as an eligible peak, it became a fencepost during peak detection, and the correct peaks to the right and 274 left were too close to be detected in successive refinement. 275 The right graph of Fig. 3 plots peak amplitude versus the length of the RR intervals immediately to the 276 left (blue asterisk) and to the right (block square) of the peak. This graph, which provides a connection 277 between peaks and RR intervals, can often provide a quick summary of the issues that a given dataset 278 might have. If the plot forms an oval shape with many RR interval lengths falling close to the central 279 quartiles (solid gray horizontal lines), the data is well-behaved. The dashed gray lines indicate the 280 robust outlier thresholds, the traditional thresholds for marking outliers on box plots. This particular 281 dataset has a lot of high negative amplitude peaks falling outside the dashed gray lines, indicating the 282 possible presence of high-amplitude signal artifacts. By default, EEG-Beats clips the amplitude on 283 these plots at 3.0×IQR outside the mid quartiles, so points clustered along the vertical edges of the 284 graph may represent very large or very small amplitude peaks. 285  In peak vs. RR plots, RR markers for RR intervals adjacent to a given peak fall on the same vertical 294 line. Uniform heart rate signals have peak vs. RR plots with closely-spaced, vertically-aligned asterisk-295 square pairs. Similarly, each RR interval is represented by a horizontally-aligned asterisk-square pair 296 corresponding to the peaks at either end of that interval. Widely separated horizontal pairs indicate a 297 dramatic difference in amplitudes between the bounding peaks of the RR interval. 298 A tip-off that there might be a problem with the peak at 620 seconds is indicated by the horizontal 299 arrow in the right graph of Fig. 3. The lone blue asterisk represents a large amplitude peak with a very 300 long interbeat interval to its right. Further, that large RR has a large negative amplitude (< −1000) peak 301 to the left and a much smaller negative amplitude (> −650) peak to the right. As a result of incorrectly 302 selecting the artifact as a peak, the RR interval to the right is much longer than that of its neighbors. 303 Thus this RR interval was correctly caught as invalid under the bad-neighbor criterion. 304 EEG-Beats is able to correctly identify peaks in the face of many artifacts such as the one that appears 305 near at 611 seconds in the left graph of Fig. 3. This artifact is of lower amplitude than the adjoining 306 heartbeat peaks, so it is not considered. However, EEG-Beats is often able to handle artifacts that are 307 much higher than the adjacent beats, provided they don't meet the sharpness criteria or fall too closely 308 to a real peak. 309 EEG-Beats has both overlay and subplot time visualizations of the RR intervals. Fig. 4 shows an 310 example of the subplot visualization for Dataset 11 (Subject S01). The top view shows the RR intervals 311 versus time as in the overlay plot, while the middle graph shows the EKG signal with peaks overlaid. 312 The bottom row has an expanded view of the EKG for a small portion of the signal, as well as the peak 313 amplitude versus RR plot on the right. The vertical arrow in the bottom graph marks an extra, low-314 amplitude EEG peak adjacent to a very small RR interval. The RR interval is marked both with a blue 315 circle (RR amplitude outlier) and a green cross (RR bad-neighbor). 316 On the surface this dataset falls into the "larger number of incorrect peaks" category as the dataset in 317 Fig. 2. In fact, this dataset just had a few incorrectly identified peaks ─ all of which were associated 318 with outlier RR amplitudes. The "bad-neighbor" criteria had many "false positives". While many of 319 the subjects in the study have slowly varying interbeat intervals, some subjects such as Subject S01 320 had heart rates that were quite variable over very short periods of time. This short-time scale 321 variability is quite clearly visible in the inset shown in the bottom row of Fig. 4. 322 with value near 218 seconds is associated with the presence of the extraneous peak. The graph in the lower right corner is 330 the peak versus RR interval graph for the dataset. The extraneous peak appears in the lower right graph as a blue asterisk 331 in the lower left corner (small amplitude peak to the right of a very small RR interval).

332
Unlike the example of Fig. 2, where the green crosses were clearly outliers from a slowly varying RR 333 interval amplitude trend, the green crosses in this case are not clearly distinguishable as outliers. The 334 bad neighbor criteria for eliminating outlier RRs works well for subjects with heart rates that are 335 relatively steady or are slowly varying over time scales larger than the "neighborhood". However, it is 336 not clear how useful or applicable such a criterion is when analyzing individuals with more unusual 337 heart rate patterns. 338 Nine of the datasets had extended periods of dropouts. An example appears in Fig. 5 for Dataset 587 339 (Subject S17). The RR artifact detection works well to remove outliers, although this dataset should be 340 trimmed before using in analysis. The peak vs. RR interval amplitude shown in the right graph clearly 341 exposes the abnormality, as the points are aligned vertically at large amplitudes. The vertical line to 342 the right would expand into a normally shaped ovoid if the signal before 100 seconds were removed 343 from the data record. 344 345 Fig. 5. RR interval overlay visualization for Dataset 587 (subject S17) of the NCTU RWN_VDE study. Caption as in Fig.   346 2. This dataset was one of the nine datasets identified as having extended periods of dropouts.

347
A final category of problematic datasets listed in Table 2 are the ten datasets in which EEG-Beats 348 picked the wrong direction for the peak orientation (either flipping when it shouldn't have or not 349 flipping when it should have). Flipping errors are usually not apparent from the EKG signal, but 350 show up in the non-ovoid shapes of peak-interval plots as illustrated in Fig. 6. 351 interval distribution changed very little since the either choice gave a consistent peak spacing, but the 359 correct orientation resulted in more accurate peak detections. (With the correct orientation, the dataset 360 moved from having a larger number to having no errors during manual verification.) The second 361 column shows the results from Dataset 339, also from Subject S05 but from a different session. This 362 dataset had peaks and troughs of roughly equal amplitude. It turned out that the amplitude of the 363 "peaks" was more variable than that of the troughs, resulting in the long tail in amplitude. When flipped 364 that dataset had no peak detection errors, while in the original orientation, it had a few errors on manual 365 review. 366 The third column of Fig. 6 (Dataset 294 Subject S06) shows an example of the third class of issues 367 EEG-Beats encountered in making a flip decision. The peaks in these datasets looked like three-lobe 368 triplets with almost equally deep troughs. The dominant amplitude within the triplets in the direction 369 EEG-Beats picked as up sometimes changed, resulting in three different spacing groups (hence three 370 horizontal lines in the plot). It happened that the dominant amplitude in the troughs was consistently 371 the central trough, so in that orientation gave a better distribution. The fourth column of Fig. 6 (Dataset 372 295 Subject S06) has the same alternation between unidirectional and peak-trough during the 373 recording. Here the two clusters are less prominent because the differences in amplitude aren't as great. 374 Several observations can be made about the graphs in Fig. 6. Since these issues did not appear in 375 recordings from other sessions with these subjects, we can conclude that the issues occurred because 376 of improper seating or electrode replacement. Notice that regardless of the orientation, EEG-Beats was 377 able to provide a fairly consistent RR distribution. The two columns associated with each subject 378 correspond to recordings made on different days ─ one with an active task and the other with a resting 379 task. Heart rates appeared to be lower in the resting tasks for both subjects. 380

Comparison of heart beat detection with the PNC toolbox 381
As a point of reference, we compared peak locations determined by EEG-Beats and the PhysioNet 382 Cardiovascular Signal Toolbox (PNC) for the first 5-minute block in each dataset. When used with its 383 default parameters, PNC failed on almost all datasets due to rejection because of atrial fibrillation 384 detection. The results in this section bypassed PNC's atrial fibrillation detection. 385 Table 3 compares peak detection results from EEG-Beats and PNC using the default parameters for 386 both toolboxes. EEG-Beats and PNC detect peaks at exactly the same places in most cases. The PNC 387 predictions that are 1 to 3 frames off from EEG-Beats are not located at the exact frame positions of 388 the signal maxima, which can happen when amplitudes are thresholded. EEG-Beats compensates for 389 this problem by adjusting the peaks to the exact peak maxima in the unthresholded signals at the end 390 of peak finding. EEG-Beats detects beats in all datasets but one that it correctly identifies as containing 391 just noise. PNC detect beats in this non-signal dataset, but was unable to detect beats in 81 other 392 datasets. Dataset 11 of Fig. 4 was such a dataset. Most traditional beat detection algorithms are not 393 prepared to deal with the artifacts and the varying peak amplitudes that are common to EEG monitoring 394 of EKG. The peak matches were computed after the ten datasets that were found to have drop-outs 395 were removed. 396 for computation of the frequency measures, however, the assumed frequency resolution differs as does 420 approach to trend removal. PNC removes the mean of the signal prior to spectral computation, while 421 by default EEG-Beats removes a polynomial (cubic) trend, resulting in smaller overall total power and 422 smaller low frequency power for EEG-Beats. However, the non-dimensional power ratios for the two 423 methods are similar. The EEG-Beats trend removal options are settable. The time measure predictions 424 differ slightly, most likely due to the additional RR artifact removal steps applied by PNC. 425

EEG-Beats visualizations 426
EEG-Beats provides scripts for creating boxplots of any RR measure segregated by the values of a 427 metadata variable such as subject or task. The user must provide the study metadata in a structure with 428 one line per dataset. The fieldnames of the structure can used to specify how to select and group metadata variables for the boxplots. Fig. 6 shows an example that was generated by EEG-Beats for the 430 LFnu measure with subject IDs as the group variable. All datasets were used with no adjustments. 431 The user must specify a particular baseline block of data for the scaling. In the case of NCTU 441 RWN_VDE, a 5-minute block of resting data acquired at the beginning of each experimental session 442 was an obvious baseline. The scaling reduces variability of the indicators across subjects. However, 443 for the RWN_VDE data collection, the sessions were scheduled at different levels of the nominal 444 subject fatigue level, and scaling by a session-specific value reduces the effect of daily fatigue level on 445 the subject's measures. A better approach for this data collection might be to scale by the average of 446 the baselines over all of a subject's sessions, but this is not supported as an automated EEG-Beats 447 feature. 448

EEG-Beats factor analysis 449
The EEG-Beats toolbox includes scripts for study-wide analysis of variance (ANOVA). As with the 450 boxplots, EEG-Beats relies a metadata structure provided by the user and can analyze an arbitrary 451 number of pairs of factors identified by the metadata structure fieldnames. The NCTU RWN_VDE 452 data set has subject, task, fatigue level (DSS measure), gender, and replicate number as possible group 453 variables. Table 5 shows an example of the output of ANOVA analysis using task and fatigue levels 454 as factors. Most of the RR measures in Table 5 show a significant statistical dependence on task and 455 on nominal fatigue level at the 0.01 significance level. Significance is usually improved with scaling. 456 The RR measures were computed from data in which RR artifactual values were removed using the 457 default settings. 458 Caution should be exercised when interpreting these results. Subject-task analysis of variance also 461 showed highly-significant dependence on those factors. However, when subject-fatigue two-way 462 analysis of variance was performed, subject dependence was found to be very highly significant, but 463 fatigue level was not. The interaction between the two was significant, however. This result has a 464 physical explanation. In this experiment, subjects were invited into the lab when the results of the daily 465 sampling system indicated that they were in a particular fatigue state. How these measurements reflect 466 on performance is highly individualized, so the fatigue-level designations are not independent of 467 subject. 468 EEG-Beats also allows three-way analysis of variance. Table 6 shows the results of subject-task-fatigue 469 analysis. Subject and task are highly significant factors for all RR measures, but fatigue has no effect. 470 This points to the importance of analyzing heart rate variability for subjects individually. The results also show that task is an important factor. In a detailed longitudinal study, Spangler et al. 474 (Spangler et al., 2020) has recently shown that performance and HRV indicators were strongly 475 individualized, but that there were consistent session-level, task-level, and block-level dependencies 476 common across subjects. 477

Discussion/Conclusion 478
No automated algorithm completely survives its first encounter with real data. In large-scale 479 computation, it is important to recognize that even adaptive approaches will fail in cases when the data 480 is bad enough. Further, even when an algorithm works well, a myriad of parameter settings may affect 481 final results. EEG-Beats focuses not only on automating well, but also on providing tools for quickly 482 assessing failures and evaluating dependencies of parameter choices on end results for an entire 483 collection. 484 EEG-Beats has utilities for beat-by-beat comparisons of two versions of the RR intervals for the same 485 dataset. Scripts provide overall summaries of agreements (such as those produced for Table 3 for the  486 EEG-Beats vs. PNC beat comparisons). These scripts can be also used on two versions of the RR 487 intervals produced by EEG-Beats for different values of any parameters to see how much the end 488 results are affected by changes in algorithm parameters. For example, we might try a lower setting of 489 the default 500 ms value of rrMinMs, (which corresponds to a heart rate of 120 bpm) to see whether 490 the default setting correctly captures all of subject records. The 500 ms settings may not be appropriate 491 for subjects recorded during vigorous activity. The point is that any suspect parameter setting can be 492 quickly tested and the results summarized, either for individual datasets or for the entire study. 493 Another example is calculating differences in measure values produced by different algorithms for the 494 same set of RR intervals as illustrated by the RR measure comparison of results for EEG-Beats and 495 PNC as presented in Table 4. Also presented in Table 4 is a comparison of EEG-Beats results with and 496 without its RR interval outlier removal. 497 The EEG-Beat visualizations are designed to enable researchers to assess quickly whether something 498 has gone wrong for a particular dataset. We use the very-large icon preview in the file browser to 499 quickly spot outliers. The RR interval versus EEG peak amplitude plots are particularly useful for 500 assessing whether a dataset might have issues. Non-ovoid distributions, distributions with multiple 501 clusters, or distributions with long trailing or leading tails merit closer inspection (Fig. 6). Almost all 502 problematic datasets can by spotted by looking at these previews. 503 The development of EEG-Beats was motivated by the potential for using heart rate variability as a low-504 cost secondary measure of subject state in EEG experiments. We started looking at EEG for EKG when 505 a careful analysis of interbeat interval information provided as output by the wrist monitors and bio-506 harness detectors also used in these experiments showed inconsistencies and numerous dropouts. We 507 developed EEG-Beats after encountering substantial difficulties in applying standard sliding-window 508 approaches to peak-finding due to EEG artifacts. EEG-Beat's top-down divide-and-conquer approach 509 to peak-finding is able to handle a variety of difficult artifactual signals in an automated fashion, but it 510 is not applicable for online applications. Further, because EEG-Beats focuses on peak detection rather 511 than detection of QRS complexes, it is not suitable for clinical applications and is designed to be used 512 on recordings from normal subjects. 513 We have shown good agreement with the well-benchmarked PhysioNet Cardiovascular Signal Toolbox 514 (PNC) for cases in which PNC can detect peaks. EEG-Beats is organized into a peak finding stage 515 (eeg_beats) that produces a structure containing detailed peak information and a computational stage 516 (eeg_ekgstats) that takes this structure and outputs a structure containing the RR measures. Scripts are 517 provided to run these functions on an entire study in an automated fashion and to perform the analyses 518 demonstrated in this paper. EEG-Beats also has an EEGLAB plugin and an associated GUI is under 519 development. EEG-Beats is freely available at https://github.com/VisLab/EEG-Beats. Structures 520 containing NCTU RWN_VDE EKG signals, metadata, and heartbeats are being released as 521 supplemental material as part of this paper to allow other researchers to compare their algorithms. 522

Conflict of Interest 523
The authors declare that the research was conducted in the absence of any commercial or financial 524 relationships that could be construed as a potential conflict of interest. 525