Supplement – Algorithm Details

There are three main parts to the MATLAB algorithm for quantifying RSA: R-peak detection, R-R interval correction, and RSA extraction. We note that the algorithm was originally developed for real-time operation on a microprocessor. Therefore, numbers were used with easy fractional equivalents (e.g. 26/32 rather than 80%).

R-Peak Detection

R-peak detection involves distinguishing the R-wave peak from other ECG features (e.g., the T-wave) and noise. A modified version of the Pan-Tompkins algorithm (Pan et al., 1985; Hamilton et al. 1986)was used. The pre-processing steps of our algorithm are similar to the Pan-Tompkins method with adjustments made to accommodate our ECG sampling rate of 256 Hz. A 33-point FIR filter with a passband extending from 5-18Hz was used to filter the data. After this initial processing, a 5-point derivative filter was used to extract the QRS slope information. The differentiated signal was squared point by point, and moving window integration (MWI) performed over a 32-point window. A training period was used to establish an estimate of the MWI peak height, which was initialized to the maximum value of the MWI signal during this period. The R-peak detection algorithm implemented in this work also used an estimate of the time between peaks (i.e., R-R interval), which was initialized to 752 ms (80 bpm). After the training period, R-peak detection was accomplished in much the same manner as described in (Pan et al., 1985; Hamilton et al. 1986).When a peak was detected in the MWI signal, a search back was performed in the filtered signal over a 25-point window starting 65 samples back from the current index; the R-peak was considered to be the maximum point in this window.

The R-peak detection algorithm implemented in this work differs from the Pan-Tompkins and Hamilton-Tompkins algorithms (Pan et al., 1985; Hamilton et al. 1986)in determining validity of detected peaks. In our approach, the R-R interval time was used as the primary factor for determining the validity of the new peak. An estimated R-R interval was determined from the previous six values, and three thresholds were established- a low limit (29/32 of the estimated value), a high limit (37/32 of the estimated value), and a missed limit (twice the low limit). A 200 ms refractory period (i.e., the time period in which it would not be physiologically possible for another heart beat to occur) was also included as a hard limit. If the new peak fell within the refractory period, it was not considered.

The MWI peak value was used as a secondary factor for determining R-peak validity. Again, an estimate of the MWI peak value was determined from the previous six values. If the new peak was detected between the low limit and high limit (i.e., the expected time), it was accepted provided that the MWI peak value was at least 1/8 of the estimated peak value. However, if the new peak value was detected outside of this range (i.e., a less likely time), a higher threshold was used. If the peak was detected after the upper limit, the MWI peak value had to be at least 1/4 of the estimated peak value. If the peak was detected before the low limit was reached, the MWI peak value had to be at least 3/4 of the estimated peak value if it was very close to the previous peak (<300 ms), 1/2 if it was in the interval corresponding to a likely T-wave (300-360 ms), and 1/4 if it was longer than the expected T-wave interval. An additional T-wave check was also performed as described in [1] in the case that the peak was detected between 300-369 ms after the previous peak.If the peak had a slope profile similar to a T-wave, then it was not accepted regardless of height. Once a new peak was accepted, the MWI peak height and R-R interval estimates were updated with the new value, and the R-R interval correction routine called.

If a peak was not detected before the missed limit was reached, a search back was performed. If a peak was found at the estimated peak location, it was accepted. If there was not a peak at the estimated value, searches were performed ahead of and behind this value, and the closest peak accepted as long as it met a minimum height criterion of 1/32 of the MWI peak height estimate. As before, once the peak was identified and accepted, the MWI peak height and R-R interval estimates were updated with the new value, and the R-R interval correction routine was called. If no peak was found, a peak was inserted at the estimated value and the R-R interval correction routine was called. In this case, the MWI peak buffer was reset to the value used from training, and the R-R interval buffer was reset to the corrected interval estimate. Additionally, the algorithm was prevented from inserting another peak for a full beat cycle.

R-R Interval Correction

The R-R interval correction routine was implemented to remove erroneous intervals that had met the acceptance criteria described in the previous section, as well as adjust actual long or short intervals that did not preserve the overall pattern of the R-R interval stream. This was accomplished by comparing the current R-R interval to a running estimate based on the previous six correct R-R intervals. The selection of bounding limits may be application specific. For this study, intervals were flagged for correction that fell outside of (26/32 – 44/32) of the running estimate. It was assumed that the need for correction resulted from either a missed R-peak (the interval was too long) or something other than a R-peak detected (the interval was too short). Regardless of the source of the error, a correction was attempted as detailed in Table S1. In order for the correction to be accepted, it was required to fall within (27/32 - 42/32) of the running estimate. Otherwise, the original interval was preserved. Additionally, a running counter was incremented for each correction and decremented for each interval accepted without correction to a minimum of zero. If this counter reached 10, the algorithm was forced to accept the next beat without attempting a correction.

In addition to the corrections in Table S1, hard limits were also enforced for the minimum and maximum acceptable R-R interval. These were selected assuming a maximum heart rate of 200 bpm and a minimum heart rate of 40 bpm. If a corrected interval was below the minimum limit, it was increased to the lower limit before being accepted, and the added time subtracted from the next R-R interval. Similarly, if the corrected interval was above the maximum limit, it was lowered to the upper limit before being accepted, and the surplus time added to the next R-R interval.

For this study, where real time operation was not required, further corrections were applied to the R-R interval stream if an interval z-score (as computed for the entire file) exceeded +/- 5. The values were averaged with an adjacent interval, resulting in a less extreme change.

RSA Computation

After corrections were applied, cubic spline interpolation was used to resample the discrete R-R interval stream at 5 Hz so that the frequency-domain analysis could be applied. The Porges-Bohrer method was used for measuring the RSA response (Lewis et al., 2012). The corrected, resampled R-R interval stream was first filtered to the RSA band frequency appropriate for young children (0.24-1.04 Hz). The data was high-pass filtered using a 21-point, third order polynomial filter before applying a bandpass filter (26-point Kaiser-window FIR filter with a passband from 0.24-1.04Hz). The natural log of the variance in 30 second, non-overlapping windows of filtered data was then computed.

Table S1: Correction algorithm actions

Current
Interval / Previous Interval / Action
Correct
(26/32 – 44/32) / Correct / Accept
Short / Add previous and current intervals if sum < 37/32
Long / Accept
Extra Long / Insert an interval with value = (current + previous) / 3
Short
(< 26/32) / Correct / Accept previous and flag current interval as short
Short / Add previous and current intervals if sum < 37/32
Long / Average previous and current intervals
Extra Long / Average previous and current intervals
Long
(> 44/32) / Correct / Accept previous and flag current interval as long
Short / Average previous and current intervals
Long / Accept previous and flag current interval as long
Extra Long / Insert an interval with value = (current + previous) / 3
Extra Long
(> 54/32) / Correct / Insert an interval with value = (current + previous) / 3
Short / Average previous and current intervals
Long / Insert an interval with value = (current + previous) / 3
Extra Long / Insert an interval with value = (current + previous) / 3

[1]Pan J, Tompkins WJ. A Real-Yime QRS Detection Algorithm. IEEE Trans Bio-Medical Eng Eng 1985;BME-32:230–6. doi:10.1109/TBME.1985.325532.

[2]Hamilton P, Tompkins W. Quantitative investigation of QRS detection rules using the MIT/BIH arrhythmia database. IEEE Trans Biomed Eng 1986;33:1157–65. doi:10.1109/TBME.1986.325695.

[3]Lewis GF, Furman SA, McCool MF, Porges SW. Statistical strategies to quantify respiratory sinus arrhythmia: are commonly used metrics equivalent? Biol Psychol 2012;89:349–64. doi:10.1016/j.biopsycho.2011.11.009.