Binary symbolic dynamics analysis to detect stress-associated changes of nonstationary heart rate variability

Psychological stress may have harmful physiological effects and result in deteriorating health. Acute psychological stress acts also on cardiac autonomic regulation and may lead to nonstationarities in the interbeat interval series. We address the requirement of stationary RR interval series to calculate frequency domain parameters of heart rate variability (HRV) and use binary symbolic dynamics derived from RR interval differences to overcome this obstacle. 24 healthy subjects (12 female, 20–35 years) completed the following procedure: waiting period, Trier Social Stress Test to induce acute psychological stress, recovery period. An electrocardiogram was recorded throughout the procedure and HRV parameters were calculated for nine 5-min periods. Nonstationarities in RR interval series were present in all periods. During acute stress the average RR interval and SDNN decreased compared to rest before and after the stress test. Neither low frequency oscillations (LF), high frequency oscillations (HF) nor LF/HF could unambiguously reflect changes during acute stress in comparison to rest. Pattern categories derived from binary symbolic dynamics clearly identified acute stress and accompanying alterations of cardiac autonomic regulation. Methods based on RR interval differences like binary symbolic dynamics should be preferred to overcome issues related to nonstationarities.

www.nature.com/scientificreports/ of the hypothalamic-pituitary-adrenal (HPA) axis as a delayed response providing long-term effects 13,14 . The sympathetic activation causes an increased release of transmitters and hormones in the central and peripheral nervous system 4,13 . This leads to changes in the organism that are necessary to facilitate a 'fight, fright or flight' response, such as elevating the metabolic rate, the blood pressure and respiration and increasing the blood flow to the heart and skeletal muscle 13 . Activation of the HPA-axis as the second stage, provides energy for a longer period of time affecting the individual's behavioral, neuronal and hormonal response to stress 13 . Especially the 'alarm' responses of acute stress, for example, the sympathetic activation and parasympathetic withdrawal of autonomic nervous functioning during acute stress, are suggested to be assessable by parameters of heart rate variability (HRV) 15,16 . Sympathetic activity of cardiac autonomic regulation is often calculated using low frequency (LF) power of spectral analysis of HRV whereas parasympathetic activity may be assessed using high frequency (HF) power or the root of the mean squares of differences between adjacent RR intervals (RMSSD) in the time domain 17 .
Methodological issues are rarely addressed in the context of the assessment of stress-related changes by means of frequency domain parameters of HRV. The application of power spectral analysis to a physiological time series requires stationary conditions of the time series. I.e., the underlying physiological system producing the RR interval series should be as constant as possible to meet this condition. Especially during stress related responses of autonomic regulation this prerequisite is rarely met. As a solution, the analysis of short term recordings using durations considerably shorter than the standard of 300 s was recently suggested 15 . However, the shorter the time series the more the different HRV measures deviate from the calculations carried out over the 300 s duration 18,19 . E.g. the amount of LF oscillations is likely to be underestimated and spectral leakage gets more prominent for shorter timer series 20 . Hence, especially the cardiac sympathetic response to stress cannot be reliably assessed using the LF component. The bias in the calculation of LF may then lead to less pronounced differences in the time course of the acute stress response.
The quantification of HRV by means of parameters derived from symbolic dynamics analysis provides solutions to overcome two of the main obstacles: (1) the coarse grained description of the RR interval series may be chosen in such a way that the nonstationarity condition does not apply anymore. (2) Appropriately chosen parameters reflecting the variability within the symbolic series require fewer data to yield proper results. It has been shown that a binary description based on the differences between adjacent RR intervals still contains sufficient information to capture the alterations of cardiac autonomic activity during a graded head-up tilt test procedure 21,22 . The occurrence of specific pattern categories could be linked to parasympathetic and sympathetic functioning of cardiac autonomic regulation.
In this study, we first explore the stationarity of RR interval series as a prerequisite for the calculation of frequency domain parameters. Data captured from 24 healthy participants during acute stress induced by the Trier Social Stress Test (TSST) is used for this purpose. Perceived stress is quantified using a visual analogue scale (VAS). The RR interval series is used to determine the existence of nonstationarities. HRV is quantified using time and frequency domain parameters as a standard set of HRV parameters. Furthermore, three pattern categories derived from binary symbolic dynamics are used to quantify HRV. The time course of each parameter is used to gain insight in alterations of cardiac autonomic regulation during acute stress.
Visual analogue scale (VAS). Across the procedure, changes in the VAS-scores as a marker for the subjectively perceived stress, were significant (p Skillings-Mack < 0.001). The median VAS-score increased from the waiting period T1 (5) to the TSST test-period, where the maximum was observed at the end of the TSST at T5 (59) directly after the arithmetic task. Subsequently, the VAS-scores decreased continuously from T5 to T9 (end of recovery period) back to the level baseline (Table 1). Correspondingly, the VAS-score at T4 and T5 were statistically different from all other analysis periods before and after the stress exposure. Gender differences were not observed. Figure 1 shows an example of a RR interval series of one subject throughout the procedure. The stress response on the RR interval series during analysis periods T4 and T5 reflects subjectively perceived stress as quantified by VAS. The decrease of RR intervals indicates elevated physiological stress. Qualitatively, in this particular example the recovery from the stress exposure takes a few minutes as can be seen in the transition from T5 to T6 which shows a lengthening of the RR intervals.

Nonstationarities in the RR interval series.
Stationary segments of the RR interval series RR i during stress periods T4, T5 and recovery period T7 are depicted in the middle row of Fig. 1. As expected, the stress response during T4 and T5 leads to several short stationary segments indicating that each entire RR intervals series during T4 and T5 is nonstationary. Unexpectedly, also the recovery period T7 shows three stationary segments, i.e. the RR interval series is nonstationary also during T7. These findings are supported by the results of the Restricted weak stationarity (RWS) test because the randomly chosen subsequences show different means (p < 0.001) and different variances (p < 0.05) during T7. The entire group also showed a varying amount of stationary segments in the course of the procedure (p Friedman < 0.01). A median of four stationary segments was found in T1, T2, T5, T7, T8 and T9. T3 and T4 had five stationary segments whereas T6 had only three segments. Accordingly, the RWS analysis also indicated nonstationarities in the different analysis periods: in 205 out of 216 analysis periods the average RR interval varied significantly in the randomly chosen subsequences indicating nonstationarities. The variance varied significantly in the subsequences With respect to gender differences only the median RR interval during stress period T4 showed a difference: male subjects had a higher median RR interval compared to female subjects (645 ms vs. 582 ms, p < 0.05). However, this difference did not lead to gender differences in any HRV parameter in the time and frequency domain nor did it lead to gender differences in the symbolic dynamics parameters.
The frequency domain parameters reflected the perceived stress to a lesser extent. LF was lowest during the stress test compared to the second part of the waiting period and compared to the end of the recovery period (T2: 7.45 ln ms 2 , T5: 6.93 ln ms 2 , 7.34 ln ms 2 , p Friedman < 0.05). However, the first part of the waiting period (T1)  Fig. 2. Short RR intervals during periods T4 and T5 indicate stress. Middle row: 5-min RR interval series during the stress periods T4, T5 and recovery period T7. The red lines indicate the median RR interval during stationary segments found by the heuristic segmentation algorithm. Bottom row: differences of successive RR intervals of the analysis periods T4, T5 and T7. The straight red lines indicate that each series is stationary. The symbolic dynamics parameters were also able to reflect changes in the course of the experimental procedure. The pattern category P0V% derived from acceleration and deceleration of RR intervals was significantly lower during the waiting period and recovery period compared to the stress test (T1: 24.4, T4: 44.4, T9: 24.5, p Friedman < 0.001, see Table 2). At the same time pattern categories P1V% and P2V% were higher during the waiting period and recovery period compared to the stress test (P1V%: T1: 60.0, T4: 44.3, T9: 58.1, p Friedman < 0.001; P2V% T1: 15.7, T4: 10.7, T9: 15.9, p Friedman < 0.001). The symbolic parameters using a threshold also reflected the course of the experimental procedure. Pattern category P0V τ % was low during the waiting period and the recovery period and increased during the stress test (T1: 38.0, T4: 76.3, T8: 34.3, p Friedman < 0.001). Pattern categories P1V τ % and P2V τ % were lowest during the stress test compared to the waiting period and the recovery period (P1V τ %: T3: 38.8, T4: 17.4, T6: 43.6, p Friedman < 0.001; P2V τ %: T1: 19.7, T4: 5.2, T8: 21.3, p Friedman < 0.001).

Discussion
The impact of different kinds of stress (e.g. mental stress, psychosocial stress) on physiological functioning such as the heart rate and heart rate variability has been investigated in numerous studies [23][24][25][26][27][28][29] . The Trier Social Stress Test as a standardized psychosocial stress procedure has also been investigated with respect to changes of HRV during the procedure 15,16,30 . However, issues arising from nonstationarities in the analyzed time series caused Table 1. Visual analog scale and heart rate variability parameters in the course of the experimental procedure. The first row shows the median, the second row shows the 25%-and 75%-percentile. The lowermost row of each parameter lists significant differences to other times (p < 0.05). # p Skillings-Mack < 0.001, *p Friedman = 0.05, **p Friedman < 0.01, ***p Friedman < 0.001 .   T1  T2  T3  T4  T5  T6  T7  T8  T9   VAS #   5  6  39  36  59  16  10  6  3   1-10  2-12  20-52  22-57  30-78  10-31  2-17  0-10  0-8   -T1  T1, T2  T1, T2, T3 T1, T2, T4  T1, T2, T3, T4, T5  T1, T2, T3, T4,  T5, T6   T2, T3, T4, T5,  T6, T7   T2, T3, T4, T5,  T6, T7 RR intervall (ms)*** 855  813  750  616  672  839  882  853  885   730-920  724-887  699-793  573-682  641-771  785-925  815-935  819-956  789-960   --T1, T2  T1, T2  T1, T2, T4  T1, T2, T3, T4, T5  T1, T2, T3, T4,  T5, T6   T1, T2, T3, T4,  T5, T6   T1, T2, T3, T4,  T5, ---T1, T2, T3 T1, T2, T3  T4, T5  T4, T5  T4, T5  T1, T4, T5 RMSSD (ms)***  T2, T3  T3  T1, T2, T3, T4 T1, T2, T4, T5  T1, T2, T4, T5  T1, T2, T4, T5  T1, T2, T4, T5 LF/HF** www.nature.com/scientificreports/ by the stress testing procedure are rarely addressed 31 . In this study, we showed that the stress testing procedure imposed nonstationarities on the RR interval series. We observed nonstationarities in the RR interval series during the stress test. The nonstationarities during stress are obviously caused by the changing demands ('nonstationary' conditions) during the stress procedure giving rise to changes of cardiac autonomic regulation and, hence, irregular trends in the time series. Surprisingly, nonstationarities were also observed during quiet rest in the waiting and recovery period although this condition would be called a 'stationary' condition. These nonstationarities in the time series may also have been caused by irregular trends arising from e.g. different depths of relaxation. I.e. although the resting condition seems to be 'stationary' it may still change during its course. We note that we did not control for breathing nor did we give any instructions with respect to relaxation during the resting periods. However, acute mental stress may lead to alterations of breathing patterns 23 and also cardiorespiratory interaction 32 . Hence, alterations of breathing patterns during acute mental stress and speech may have contributed to alterations of cardiac autonomic regulation and may also lead to nonstationarities. Furthermore, very low frequency fluctuations linked to e.g. vagal baroreflex sensitivity may also lead to nonstationarities of the RR interval series 33 because the shortest segments contained 40 RR intervals approximately equivalent to the threshold between very low frequency and low frequency oscillations in the frequency domain. The Fourier transformation requires stationarity of the underlying time series. Hence, nonstationarities in the RR interval series have an impact on the calculation of the HRV parameters in the frequency domain. A comparison between parameters calculated from Fourier analysis and parameters calculated from the wavelet transformation, which can be used for nonstationary time series, showed only small differences 34 . Nevertheless, the differences may attenuate the variance across the stress testing procedure for parameters of the Fourier transformation. The LF parameter showed relatively little variance across the waiting, stress and recovery periods. During stress period T5 LF was lowest. However, the decrease of LF is contrary to what would be expected: the stress periods T4 and T5 elicit sympathetic activation as reflected by the decrease of the median RR interval and SDNN. At the same time parasympathetic activity decreased as indicated by RMSSD 17 . As LF is affected by sympathetic as well es parasympathetic influence, LF in this particular case does not reflect the sympathetic activation but seems to reflect only the decrease of parasympathetic activity. LF% as an equivalent to LF expressed in normalized units should reflect solely sympathetic activity 17 . Indeed, LF % showed an increase during T5 indicating the sympathetic activation correctly. Of note is that during a tilt testing procedure eliciting cardiac sympathetic activation the LF parameter also performed worse compared to LF% 35 . Hence, LF% may be better suited as a parameter reflecting sympathetic activation although it did not reflect the increase in stress during stress period T4.
HF was consistently lower during T5 compared to all waiting and recovery periods indicating diminished vagal activity during stress. This result is in accordance to RMSSD in the time domain which also reflects parasympathetic activity. Furthermore, HF% was also low during stress period T5 but it was not different from the waiting periods T1 and T2. Hence, in this case HF is superior compared to HF% because the results of HF are Table 2. Symbolic dynamics parameters in the course of the experimental procedure. The first row shows the median, the second row shows the 25%-and 75%-percentile. The lowermost row of each parameter lists significant differences to other times (p < 0.05). ***p Friedman < 0.001.
The parameters derived from symbolic dynamics are not biased by nonstationary RR interval series because the associated series of differences between successive RR intervals is stationary although the underlying RR interval series may be nonstationary. P0V% and P0V τ % consistently increased during stress periods T4 and T5 compared to waiting and recovery periods. It has been shown that P0V% and P0V τ % may be interpreted in terms of sympathetic activity, i.e. the higher P0V% and P0V τ % the higher the sympathetic activity 21,22,37 . Hence, these parameters indicate the time course of low sympathetic activity during waiting and recovery periods as well as sympathetic activation during acute stress. On the other hand, P1V% and P1V τ % consistently decreased during stress and were lower compared to waiting and recovery periods. As these parameters indicate parasympathetic activity, i.e. the higher P1V% and P1V τ % the higher the parasympathetic activity, also the course of parasympathetic activity could be clearly captured by these parameters. P2V% and P2V τ % also decreased during stress periods T4 and T5 compared to the waiting and recovery periods. However, these pattern categories could not be unambiguously linked to sympathetic or parasympathetic activity 21,22,37 . Hence, these parameters remain unclear with respect to the interpretation in terms of cardiac autonomic regulation.
We note that gender differences were only observed during the stress period T4: although the perceived stress was similar, male subjects had a higher median RR interval compared to female subjects. However, this difference did not affect e.g. the amount of nonstationarities nor did it affect any parameter in the time and frequency domain or parameters from the symbolic dynamics analysis. Recent studies showed that female subjects had higher perceived stress and lower HRV parameters compared to male subjects 28 . These differences could be attributed to gender differences in e.g. coping styles and emotion regulation strategies 38 . In the present study gender differences were not observable most likely due to the relatively homogenous study population (University students) and the participation dates of female subjects that had to be in the second half of their menstrual cycle to decrease gender differences 39 .
In conclusion, nonstationarities in RR interval series occur during transient states like e.g. acute stress but our analysis showed that also during resting and quiet states nonstationarities have to be expected. From this perspective, the quantification of HRV with frequency domain parameters is limited in almost any case because the requirement of stationarity is rarely met. The nonstationarities obviously biased the LF and HF and, as a result, these parameters were less informative than e.g. the RR interval series or SDNN because the latter parameters clearly indicated acute stress whereas LF and HF did not. Methods based on the differences of RR intervals like e.g. parameters derived from binary symbolic dynamics are more informative compared to frequency domain parameters. They do not depend on nonstationarities and they can also be interpreted in terms of cardiac autonomic regulation and, hence, allow physiologically meaningful interpretations.

Methods
This study is part of a multi-faceted research project, investigating different aspects of stress-responsive processes to gain a comprehensive approach to acute psychological stress perception and response. We investigated psychological parameters of subjective stress perception, physiological stress parameters such as HRV-parameters, biochemical stress parameters (salivary cortisol levels and salivary alpha-amylase activity) as well as epigenetic parameters such as salivary microRNAs (miRNAs). For further details on other aspects of this project, we would like to refer to our recent publications 8,9 . Subjects. 24 healthy subjects (12 female, age: 20-35 years) were recruited among a population of University students. The female subjects had to be in the second part of their menstrual cycle at the date of participation because sex differences were observed in stress responsibility which apply mainly to the first half of the menstrual cycle 39 . Furthermore, female participants did not take any hormonal contraception. All participants were in good physical and psychological health, had no history of psychiatric diseases, were non-smoking, taking no drugs, alcohol or medication, and did not perform any type of meditation or relaxation-exercises regularly (more than once a week). The female participants did not take any hormonal contraception and were in the second half of their menstrual cycle at the date of participation. To minimize interference with circadian variations of cortisol levels, all tests were carried out between 3 and 5 p.m 9 .
All experiments were conducted in accordance with the Declaration of Helsinki. The study was approved by the ethical committee of Witten/Herdecke University, Germany (96/2015) and registered in the German Register for Clinical Studies DRKS which is linked to the WHO-Register (Registration-ID: DRKS00010134). The participants were informed in a written and oral format about the study aims and procedures, particularly their participation in a psychosocial stress test, and the timetable of the experiment was explained. To prevent reduced stress reactions due to prior mental adaption to the expected task, participants were not informed about the specific details of the TSST. Written informed consent of each participant was obtained before the onset of the experiment. They were debriefed at the end of the experiment receiving full information about the procedure 9 .  10 was applied as a reliable and standardized acute psychosocial stress test, mainly following the TSST-protocol by Birkett, 2011. Time-points of data collection were modified to suite the individual requirements of the study (Fig. 2). The stress test was split in three periods: a waiting period of 30 min, the actual stress test period of 20 min, and a recovery period of 60 min. During the whole procedure, the participants were not allowed to use any electronic media. During the waiting and recovery period, the participants could relax being located on their own in a quiet atmosphere with comfortable seating. The actual TSST test period took place in the social laboratory, a simply equipped room containing a chair for the participant during the preparation phase as well as an office desk and two chairs for the 'experts' . There, the participants had to prepare (10 min) and then deliver (5 min) an oral presentation applying for an individually ideal job-offer. Afterwards, they had to perform a mental arithmetic task (5 min), sequentially subtracting 13 from 1,022. The speech delivery and mental arithmetic task took place in front of a panel of two persons as 'experts' , who followed a strict protocol. They wore white lab coats, exhibited unemotional neutrality, avoided any oral or mimic feedback, and just adverted if there was still time remaining or if a mistake was made during the mental arithmetic task, instructing the participant to start again from 1,022. Furthermore, dummies of a video camera and a microphone were installed and the participants were told to be recorded during their speech delivery and mental arithmetic task.
Heart rate variability. An electrocardiogram (ECG) was recorded continuously throughout the procedure using a portable Holter recorder (TOM Medical MK3, Graz, Austria). Time markers were set concurrently to the VAS at nine time points (T1, … T9) to enable proper identification of the previous 5 min as respective analysis periods (Fig. 1).
The Holter device's sampling rate of the ECG was 4,096 Hz. Hence, the internally detected times of R-peaks had a precision < 1 ms. The times of the R-peaks and the ECG at a sampling rate of 256 Hz were saved on a memory card. The detected times of R-peaks were visually checked and corrected in case of false detections due to e.g. artifacts (< 1% of all detected R-peaks). Subsequently, the RR interval series was calculated as the temporal distance between successive R-peaks. Times of ventricular and supraventricular beats were replaced by appropriately interpolated times 41 . The resulting RR interval series RR i (i = 1, . . . , N) served as the basis for the HRV analysis.
The RR interval series of each 5 min analysis period was quantified as follows. The average RR interval, its standard deviation (SDNN) and the root of the squared mean difference between successive RR intervals (RMSSD) were calculated as basic parameters in the time domain. The median length of the RR interval series varied between 338 (T9) and 483 (T4) RR intervals. The calculation of the frequency domain parameters was carried out using a re-sampled time series at 4 Hz. The re-sampled time series was detrended and a Hanning window was applied. Low (LF: 0.04-0.15 Hz) and high frequency (HF: 0.15-0.4 Hz) oscillations and the fraction LF/HF were quantified using a fast Fourier transformation (2048 data points with zero padding) 17 . The total spectral power was adjusted to the variance of the RR interval series and, hence, HF and LF were expressed in ms 2 . The proportions of LF and HF in relation to the total power, LF% and HF%, were also calculated because these quantities tend to minimize the impact of changes in total power 17 . Hence, they may better reflect changes of cardiac autonomic regulation 21,35 . Symbolic analysis. The RR interval series RR i (i = 1, . . . , N) is transformed into a binary symbolic series by two different approaches. The first approach simply reflects the succession of acceleration and deceleration of heart rate. I.e. the difference series �RR i = RR i − RR i−1 (i = 2, . . . , N) is calculated and the symbolic sequence is created according to the sign of each difference 37 : www.nature.com/scientificreports/ The 0 s represent decelerations of the heart rate whereas the 1 s represent accelerations (see Fig. 3, left column).
In this approach, no parameter has to be chosen. In the second approach the difference series ΔRR i is transformed into a binary series according to a predefined threshold. The binary coding represents whether the absolute value of the difference ΔRR i is below or above the threshold τ: This transformation reflects whether the succession of RR intervals has only small changes (0 s) or also contains larger changes (1 s, see Fig. 3, right column). In this study, the threshold τ is set to 35 ms, i.e. approximately 5% of the grand average RR interval. This threshold resulted in binary series with sufficient changes between 0 and 1 s reflecting changes in the dynamics of the RR interval series.
We calculated the relative frequency of each pattern category for both symbolic sequences (P0V%, P1V%, P2V% and P0V τ %, P1V τ %, P2V τ %). It has been shown that the categories P0V%, P0V τ % and P1V%, P1V τ % properly reflect sympathetic and parasympathetic modulations of cardiac autonomic regulation, respectively 22,37 . Analysis of nonstationarities in the RR interval series. Methods such as the Fourier transformation to calculate LF and HF require stationarity of the analyzed time series. A time series is stationary if its statistical characteristics (e.g. the mean, the standard deviation and all higher moments) are invariant with respect to time translation. However, especially in time series like the RR interval series during stress this prerequisite is seldom met. To quantify the amount of nonstationarities in each analysis period T1, …, T9 two different approaches are used. The first approach uses a heuristic segmentation of the time series. Stationary segments are created utilizing the pooled variance of two adjacent segments in such a way that the mean value between adjacent segments is maximized 42,43 . We set the minimum length of stationary segments to 40 RR intervals because shorter segments could lead to a segmentation caused by fluctuations in the low frequency band. The amount of segments in each analysis period is used as an indicator of nonstationarities. In a stationary analysis period the RR interval series must not be segmented. www.nature.com/scientificreports/ The second approach is called the restricted weak stationarity (RWS) test and checks the mean and the variance of randomly chosen subsequences of the time series under consideration 44 . In case of stationarity all subsequences should have the same mean and variance. This approach quantifies separately the probability of different means and different variances in the subsequences. In stationary RR interval series these probabilities should be p > 0.05. As suggested by the authors we used 8 randomly subsequences containing 50 RR intervals for this approach.
Statistical analysis. All statistical procedures are descriptive. Non-parametric statistical procedures were used because of the low number of subjects. The distributions of LF, HF and LF/HF showed skewed distributions and, hence, they were transformed taking the natural logarithm. The distribution of each parameter was quantified by the median and the interquartile range (25% and 75%-percentile). Non-parametric statistical procedures were consistently used. The Friedman test for repeated measures was used to assess changes of the parameters at the times T1-T9. In case of missing values (VAS), the Skillings-Mack test was used as a replacement for the Friedman test to take advantage of all available data 45 . If the Friedman test (or the Skillings-Mack test) showed significant changes of a parameter, pair-wise differences between different times were checked post-hoc including adjustment for multiple comparisons 46 . Gender differences were tested using the Mann-Whitney U-test. A p < 0.05 was considered statistically significant. www.nature.com/scientificreports/