Complexity of cardiac signals for predicting changes in alpha-waves after stress in patients undergoing cardiac catheterization

The hierarchical interaction between electrical signals of the brain and heart is not fully understood. We hypothesized that the complexity of cardiac electrical activity can be used to predict changes in encephalic electricity after stress. Most methods for analyzing the interaction between the heart rate variability (HRV) and electroencephalography (EEG) require a computation-intensive mathematical model. To overcome these limitations and increase the predictive accuracy of human relaxing states, we developed a method to test our hypothesis. In addition to routine linear analysis, multiscale entropy and detrended fluctuation analysis of the HRV were used to quantify nonstationary and nonlinear dynamic changes in the heart rate time series. Short-time Fourier transform was applied to quantify the power of EEG. The clinical, HRV, and EEG parameters of postcatheterization EEG alpha waves were analyzed using change-score analysis and generalized additive models. In conclusion, the complexity of cardiac electrical signals can be used to predict EEG changes after stress.


Results
Study patient demographics. A total of 117 patients were enrolled in the study. All patients tolerated the procedure well, and no clinical complications were noted during the index procedure or admission. Among them, thirty-three patients were excluded because of poor data recording and incomplete EEG recording, and the remaining 84 patients were selected for the final analysis. The clinical information of all patients (n = 84) are presented in Table 1. The mean age was 64.1 ± 11.8 years, and 69 were men. Sixty-five patients underwent both cardiac catheterization and coronary artery intervention (including stent implantation) in the index procedure (revascularization group). Nineteen patients underwent cardiac catheterization alone (nonrevascularization group) owing to the lack of significant stenosis in the index procedure.
Serum neurotransmitter measurement and ECG/EEG parameters in the revascularization and nonrevascularization groups. The differences in the serum neurotransmitter levels, HRV metrics and EEG parameters are presented in Table 2. In terms of serum neurotransmitters, the dopamine level was high in the revascularization group (231.60 ± 224.09 vs. 148.96 ± 201.66; p = 0.047). The differences in the orphanin FQ and serotonin levels between the groups were not statistically significant. In the EEG reading, the pre-procedure alpha activity, post-procedure alpha activity, and difference in the pre/post procedure alpha activity were comparable in both groups. The HRV metrics of pNN20 (pre-) and log-pNN20 (pre-) in the revascularization group were significantly different (0.4974 ± 0.2553 vs. 0.3567 ± 0.2267.66; p = 0.026 and − 0.87 ± 0.68 vs. − 1.27 ± 0.77; p = 0.026). The other HRV parameters were comparable in both groups.
Scientific RepoRts | 5:13315 | DOi: 10.1038/srep13315 Experimental framework. Figure 1 shows the proposed framework, which consists of multiple analysis steps, including an overall variation analysis method (time-and frequency-based analysis), statistical method, and clinical data. Each ECG and EEG segment, with a time series of 4 min, was used for signal processing and statistical analysis. The novel ideal of the procedure is emphasized, and details of the procedure are presented in the Methods section. This procedure is completely based on the linear and nonlinear analysis of the pre-and post-test heartbeat. To improve the procedure and enhance the accuracy of predicting human's relaxing states, we included the post-catheterization EEG alpha waves and clinical data. In particular, we used STFT when there was a change in the encephalic electricity.
The statistics on the heartbeat and brainwaves were obtained using change-score analysis. Before model fitting, the continuous ECG and EEG parameters were analyzed using Spearman's rank correlation to remove a confounding variable, and GAMs were used to increase the prediction accuracy.
Data classification and definitions. ECG parameters were computed using two metrics, i.e., linear and nonlinear methods, as shown in Supplementary Table S1. Some HRV metrics showed strong interactions between HRV and EEG, with a correlation coefficient of 0.8 (see Supplementary Table S4 and S5 for details). The correlation between MSE and conventional HRV analysis was nonsignificant for all patients. The EEG alpha activity was determined using STFT. The experiment involved the use of 52 continuous variables to investigate the role of relaxing states (post-test alpha activity). To enhance personalization and improve the investigation accuracy, the change score analysis derived 11 discrete variables of clinical data and 3 continuous serum neurotransmitter variables. Before multiple regression model fitting, the continuous parameters of ECG and EEG were analyzed using Spearman's rank correlation to remove confounding variables from the test for continuous variables, which resulted in 23 continuous variables. Supplementary Table S2 shows the pretest, and the difference between the pre-and post-test was considered after determining the Spearman's rank correlation. To determine the prognostic value for the recognition of relaxed states, we proposed the following three methods: change-score analysis, change-score analysis using GAMs, and change-score analysis using GAMs without considering the EEG data obtained before stress (cardiac catheterization).
Change-score analysis. Table 3 shows the 84 patients who were analyzed using the change-score analysis. Six multiple regression parameters were evaluated. The pre-alpha activity, pre-meanNN, and Area 6-20 (pre-) 14.7 ± 6.8 15.7 ± 6.5 14.4 ± 6.9 P = 0.4938  Table 2. Effect of the revascularization treatment on the autonomic activities, brain waves and serum neurotransmitter.
difference-meanNN were the three main predictive factors used to recognize relaxing states on the basis of difference-alpha activity. Moreover, the pre-test alpha activity could influence the occurrence of difference-alpha activity because of its direct effect on the difference-meanNN and its indirect effect on pre-meanNN. Furthermore, the age, pre-log-Slopes 1-5, and pre-Slopes 1-5 were minor predictive factors for the difference-alpha activity. The results of the six regression parameters shed light on the mechanism of the major predictive factors of the difference-alpha activity, and the multiple R-squared value was 0.6509. This proposed method achieved an overall accuracy of 80.7% in recognizing changes in the EEG alpha waves.
Change-score analysis using generalized additive models. The GAMs were used to detect the nonlinear effects of continuous covariates and identify appropriate cutoff points for discretizing a continuous covariate. GAMs enable the direct observation of the partial effect of a continuous covariance; therefore, the cutoff points can be easily selected. Supplementary Table S3 shows the continuous ECG parameters that were analyzed using the GAMs. The GAM of all patients is shown in supplementary Figs S1 to S19. After stepwise variable selection, five GAM variables were included in the change-score analysis, and the multiple R-squared value increased to 0.7485, as shown in Table 4. The pre-alpha activity and difference-meanNN were the predominant predictive factors of difference-alpha activity. Furthermore, the difference-meanNN, difference-the ratio of LF over HF and GAM's HF were minor x values for difference-alpha activity. The proposed method achieved an overall accuracy of 86.5% in recognizing changes in EEG alpha waves.   Change-score analysis using generalized additive models without electroencephalography data obtained before stress (cardiac catheterization). This model was developed for wearable sensor nodes that, based on change-score analysis, identified linear and nonlinear indices of HRV, hypertension, and neurotransmitter levels without considering the electroencephalography data obtained before stress. As shown in Table 5, 12 multiple regression parameters were evaluated. Slopes 1-5 were the main predictive factors of the difference-alpha activity. In contrast to the change-score analysis and the change-score analysis using GAMs, nonlinear analyses are useful for predicting the difference-alpha activity. Furthermore, the difference-meanNN and difference-the ratio of LF over HF were minor predictive factors for the difference-alpha activity. In addition, the serum neurotransmitter levels, hypertension, and serotonin were included in the model. After stepwise variable selection, the multiple R-squared value was 0.4988. The proposed method achieved an overall accuracy of 70.6% in recognizing changes in EEG alpha waves. Although the prediction accuracy was relatively low, continuous EEG parameters need not be considered. Hence, HRV analyses are sufficient for assessing the autonomic function.

Discussion
In movement control, the brain commands the motor system to produce the desired movements. The brain appears to control the motor system, but the afferent input from the motor system may strengthen or weaken the synaptic connections, thereby reorganizing the brain and leading to functional and structural changes in the brain. We posited that the brain interacts similarly with the heart. In our previous study 17 , we demonstrated that there are correlations between the signal complexity of cardiac and cerebral activities. In the present study, we used cardiac catheterization as a stress event to evaluate the predictive power of post-stress brain activity by considering the ECG complexity parameters. As mentioned previously, EEG alpha waves represent a relaxing index 28 . Therefore, we developed three methods, which primarily involve ECG and EEG parameters, for assessing the human relaxing state after stress. A numerical study of estimating the human relaxing state was developed in the literature 4 . The present method requires purely parametric values, and the human relaxing state can be calculated dynamically. This method mainly requires HRV metrics and the instantaneous power spectrum of EEG during the estimated relaxing states. Because this method was further incorporated with clinical data, such as the personal profile, diagnosis or serum neurotransmitter, it should be able to accurately characterize the relaxation response.
Two nonlinear methods, MSE and DFA, were used to quantify the heartbeat complexity and variability of heartbeat time series to provide information about the autonomic function. Because the connection between brain and cardiac activity has been demonstrated to be highly complex, HRV is a conventional metric for the pathophysiology and psychopathology 29 . Although conventional HRV metrics are reliable for predicting human relaxing states, the prediction accuracy is limited in the presence of frequent arrhythmias 30 . In this study, MSE and DFA were used to derive nonstationary and nonlinear dynamic changes from the heart rate time series 17,31,32 . The short-fitted slope (Slopes 1-5) of the curve obtained using MSE was considered when the change-score analysis was used in the change-score analysis. Moreover, when we used generalized additive models to increase the prediction accuracy, the short-fitted area (Areas 1-5) and long-fitted area (Areas 6-20) were used. Compared with conventional linear HRV indices, the nonlinear indices of the employed method had higher prediction ability for cardiac catheterization.
Previous studies have indicated that alpha activity reflects the selective cortical inhibition that is involved in global neural integration 33,34 . This feature of alpha activity has received considerable attention in cognitive neuroscience. The cognitive literature suggests that alpha activity is potentially involved in various other physiological processes, such as relaxing, or in specific regions of the brain 35 . Therefore, we hypothesized that alpha power is correlated with the relaxing state because the effect of neural inhibition could change emotions according to the physiological and pathological states. Therefore, short-time Fourier transform was used to quantify the spectral power of alpha activity from a 4-min time series of alpha activity. Our results showed that when the post EEG alpha waves were considered in the change-score analysis, the predictive accuracy increased from 70.6% to 80.7%.
Although we employed change-score analysis to estimate the human relaxing state after stress, numerous nonlinear HRV metrics (listed in Supplementary Table 3) were not considered. Consequently, the generalized additive models were fit to detect the nonlinear effects of continuous covariates 36 . In this study, generalized additive models were used to visualize the relationship of nonlinear HRV metrics with the difference between the pre-and posttest alpha wave responses. On the basis of graphical visualization, nonlinear HRV metrics can be used to identify the appropriate cutoff points for discretizing a continuous covariate (see Supplementary Fig. S1-S19). Our results for the ratio between the LF and HF and for the HF were considered in the prediction methods after fitting the generalized additive models. A comparison between the proposed prediction model and the proposed prediction model, without the generalized additive models, showed that the time-domain HRV metrics were the only parameters that could be considered to predict human relaxing states.
In this study, an innovative aspect of the method is the use of HRV and EEG to fit the regression model and the observed data for effect estimation and outcome prediction. Three methods were used for assessing the difference between the pre-and post-test alpha wave response. A change-score analysis model identified six mathematical parameters in the HRV and EEG that may facilitate predicting post-stress alpha waves. Conventional linear HRV analysis involves two metrics, the meanNN (pre-) and meanNN (difference), which provide positive and considerably low estimates. In addition, MSE analysis involves two metrics, Slopes 1-5 (pre-) and Log-areas 1-5 (pre-). Among the clinical parameters, only age was associated with high predictability. Furthermore, the change-score analysis model using GAMs identified 11 regression metrics, and five GAM variables were included in the predictive model. All GAM estimated values of the HRV analysis were positive at specific cutoff points. In addition, five predictive parameters that had been identified in previous studies were retained. Considering the practical application on wearable sensor nodes 11 , the EEG data obtained before stress (cardiac catheterization) were not used in the change-score analysis. Five regression metrics are used in conventional linear HRV analysis. In the present study, sdNN(pre-), pNN20 (pre-), and LF/HF (difference) were new metrics that were added to the predictive model. The variables associated with clinical parameters included hypertension. An estimated value of serotonin was, for the first time, considered personalized in change-score analysis. Drugs that alter the serotonin levels are used to treat depression and generalized anxiety disorder. Therefore, it is reasonable for serotonin to be included in this model. This model shows that when we Scientific RepoRts | 5:13315 | DOi: 10.1038/srep13315 did not include the prestress alpha activity, the difference between the pre-and posttest alpha wave responses require more parameters to estimate the relaxing states. Of note, three predictable models did not include the treatment category of clinical data, which means that our proposed method can be applied to healthy subjects.
This study had numerous limitations. First, there was no control group in the entire relaxing condition, such in a home setting. All study subjects were patients with a relative anxiety status. Therefore, another study is needed to investigate the relationship between the brain and heart in subjects with a relaxing condition. Second, all ECG and EEG data were recorded in normal "free running" conditions, which implied the possible existence of interfering factors (e.g., physical activities and different relaxing states). Third, the database was based on data collected for a period of 1 h. Hence, dynamic noise or nonstationary artifacts might have affected the signal properties. Although the segments of ECC and EEG were detrended using a deterministic nonlinear method (empirical mode decomposition) discussed in previous studies 37 , we did not assess the noise features or noise level to obtain more detailed information. Finally, some parameters related to the inhibitory characteristics reflected by EEG, such as cortical inhibitory mechanisms, were not considered in this study 38 . Therefore, it is necessary to develop a specific mathematical algorithm for quantifying global neural integration.
In conclusion, the method of predictive analysis demonstrates the potential of the interrelations between heartbeat and neural communication networks by using an accurate quantification equation involving change-score analysis, which is lacking in most previous studies. Accordingly, to increase the probabilistic accuracy, we developed a method to test our hypothesis. In this method, although conventional HRV analysis enables the quantification of a heartbeat, nonlinear methods provide additional information about the heartbeat complexity and fractal correlation. The proposed method achieved an overall accuracy of 80.7 and 86.5% in recognizing alpha-wave changes after stress basis on the change-score analysis and GAMs, respectively. Therefore, the complexity of cardiac electrical signals could be used to predict changes in the EEG alpha waves after stress.

Methods
The ethics committee on human research of National Taiwan University Hospital (NTUH) approved the study (IRB-201206070RIC). All participants or their surrogates gave written informed consent. The investigation conformed to the principles of the Declaration of Helsinki.
Data collection. Patients with angina and positive stress tests who had undergone cardiac catheterization uneventfully were enrolled in this study from August 2012 to May 2014. One day before and immediately after cardiac catheterization, these patients underwent 1-hour EEG (Neuron-Spectrum-3, Neurosoft Company; Digital Neurophysiological Systems, Russia) monitoring with a sampling rate of 512 Hz and resolution of 16 bits. Continuous ECG recordings were extracted from the resting awake EEG for each patient. The RR interval (RRI) of ECG was determined using an automated detection algorithm, and the annotated file was carefully inspected and corrected for the extraction of the RRIs. Surface EEG was conducted using 19 electrodes of the international standard 10-20 system (FP2, F3, F4, FZ, C3, C4,  CZ, P3, P4, PZ, O1, O2, F7, F8, T3, T4, T5, and T6). All patient data, including artifacts, such as eye movements, blinks, and muscle activities, were saved in text files for offline analysis on a laptop, and 4-min segments of the time-series were chosen from each file. Medical history, including demography and medication, was carefully recorded. Blood was sampled during cardiac catheterization.

Time-and frequency-based analysis.
Conventional HRV metrics were calculated in the time and frequency domains according to the guidelines developed by the Task Force of the European Society of Cardiology 12 .
The patients were diagnosed using metrics in the time and frequency domains. The time domain metrics were as follows: the mean value of RRI time series (meanNN), standard deviation of the RRI time series (sdNN), root mean square of the differences between successive RRIs (rMSSD), percentage of absolute differences greater than 20 ms in normal RRIs (pNN20), and percentage of absolute differences greater than 50 ms in normal RRIs (pNN50). The frequency domain metrics were as follows: high-frequency component (HF, 0.15-0.4 Hz), low-frequency component (LF, 0.04-0.15 Hz), and the ratio between the low-frequency and high-frequency components (LF/HF), which were computed according to the average power spectrum.
Multiscale entropy analysis. MSE analysis can be used to estimate the complex pattern of a time series and evaluate the complexity of physiological signals on multiple time scales 21 .
The MSE analysis method consists of two main steps: coarse-graining of the signals into different time scales and quantifying the degree of irregularity in each coarse-grained time series using sample entropy (SpEn) 39 . Coarse-graining yields a one-dimensional discrete time series {{x 1 ,…, x i ,…, x N }. In addition, the consecutive coarse-grained time series, {y τ }, which corresponds to the scale factor τ, can be constructed. To construct a coarse-grained series, the original time series should first be divided into N/τ nonoverlapping windows of length τ. Subsequently, the data points within each window should be averaged. Generally, each element of the coarse-grained series is calculated as follows: Scientific RepoRts | 5:13315 | DOi: 10.1038/srep13315 With a scale factor of one, which corresponds to τ = 1, the coarse-grained series ( ) y j 1 is exactly the same as the original one.
In addition, the MSE curve can be divided into long-and short-term curves. The different complexities of the scales can be useful for clinical categorization and can be divided into four types of parameters: (a) when only slope is considered, the first 5 scales are defined as the short-fitted slope (Slopes 1-5), and the final 15 scales are defined as the long-fitted slope (Slopes 6-20); and (b) when only the area is considered, the first 5 scales are defined as a short-fitted area (Areas 1-5), and the final 15 scales are defined as a long-fitted area (Areas [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. Although MSE analysis has been applied to physiological signals, the entropy values are sensitive to very-low-frequency noise and nonstationary artifacts, especially trends. Accordingly, the de-trending process should be performed prior to MSE analysis 37 . In the present study, the ensemble empirical mode decomposition method was adopted to remove the trend of the RRI signals. Detrended fluctuation analysis. DFA can be used to quantify nonstationary and nonlinear dynamic changes in the heart rate time series. Recently, DFA has been widely used to analyze cardiac problems by inputting patients' RRI signals 40  The short-term (< 11 beats, alpha_1) and long-term (> 11 beats, alpha_2) fractal correlation exponents were calculated to obtain a clearer understanding of the fractal correlation property of the biological system 23 . In addition, the heartbeat dynamics were characterized by a scaling exponent ∝ , and the slope of the linear relationship was estimated according to the log-log plot of fluctuations versus box sizes.
Electroencephalography analysis using short-time Fourier Transform. In STFT, a short time window is applied to biomedical signals 41,42 and a series of Fourier transforms are performed within this window as the window slides across all data; STFT provides a time-frequency representation of the biomedical signal (Fig. 2). STFT can provide an instantaneous estimate of the time-varying energy because the Fourier transform can be adapted for analyzing a localized signal. Mathematically, STFT is defined as follows: } 4 x j ft 2 Figure 2. STFT along with a spectrogram and fixed window size can be used for localizing signals. w(t) and x(t) denote the window and EEG signal, respectively, and the EEG signals are expressed in the timefrequency domain.
where w(t) is a symmetric window function and x(t) is the signal to be transformed. STFT x (τ, f) uses a complex exponential basis function as the basis. The squared magnitude of STFT is referred to as a spectrogram: where τ is discrete and f is continuous. All analyses were performed using MATLAB (The MathWorks, Natick, MA, USA). The data were low-pass filtered with a cutoff frequency of 55 Hz to remove power line noise, and the power changes in the EEG following movement execution were estimated for all patients using a fast Fourier transform with windows of 1024 samples, a Hanning window with a width of 0.5 s, and a 70% overlap until all the EEG signals were analyzed.
The objective of using STFT is to assess the relaxing state of a subject exclusively through EEG analysis. Hence, the alpha activity is an essential metric for assessing relaxed wakefulness when visual processes are engaged by opening the eyes 43 . In the present study, we calculated only the alpha-activity-ratio related EEG signals using time-frequency analysis. The alpha-activity-ratio (AAR) obtained through EEG was investigated by comparing the alpha band (8)(9)(10)(11)(12) to the full-band EEG (1-55 Hz). Mathematically, the AAR is defined as follows: where E(t k , f i ) is the spectrum of the i-th bin in STFT, t k denotes the time indices, and N is the total time length. On the basis of the obtained AARs, a significant level (AAR th ) can be identified to quantify the intensity of the alpha waves.
In this study, we proposed a framework for calculating significant changes in the average time-frequency power density. This method consists of the following steps: (a) calculating the time-frequency power density, (b) dividing the time-frequency plane in which the alpha wave power density is calculated, (c) selecting a threshold for all channels, and (d) determining the energy change at the significant level for all channels. The process flow block diagram is presented in Fig. 3.
Serum neurotransmitter measurement. The levels of serum neurotransmitters, including orphanin FQ, dopamine, and serotonin, were measured. The orphanin FQ level was measured using an ELISA Kit (MyBioSource, San Diego, CA, USA), the serotonin level was measured using a serotonin ELISA Kit (Abcam, Cambridge, United Kingdom), and the dopamine level was measured using high-performance liquid chromatography. Statistical analysis. Statistical analysis was performed using the R 3.1.1 software (R Foundation for Statistical Computing, Vienna, Austria). A two-sided p value less than 0.05 was considered statistically significant. Data were expressed as the mean ± standard deviation. Differences in the serum neurotransmitter levels between the groups of patients were assessed using the Mann-Whitney U Test. Continuous data were assessed using the Mann-Whitney U Test, and categorical data were assessed using Fisher's exact test. The statistical analysis consisted of Spearman's rank correlation, a stepwise variable selection method, and change-score analysis.
The objective of conducting a regression analysis of the change scores (i.e., the difference between the pre-and posttest response) was to determine one or a few parsimonious regression models that fit the observed data well for effect estimation and outcome prediction. To ensure that the analysis quality was high, basic model-fitting techniques for variable selection, change-score analysis of pre-and posttest data, and regression diagnostics and remedies (e.g., ensuring that multicollinearity was not present and detecting influential cases) were applied. Before model fitting, the continuous ECG and EEG parameters were analyzed using Spearman's rank correlation to remove confounding variables. The significance level for the confounding effects was set to be greater than 0.8.
The stepwise variable selection procedure (with iterations between the forward and backward steps) was used to identify the optimal candidate final regression model. All relevant univariate significant and nonsignificant covariates (listed in Tables 1 and 2) and some of their interactions were included in the variable list to be selected. The significance level for entry and stay was conservatively set to 0.15. Subsequently, the optimal candidate final regression model was manually identified by dropping the covariates with p values greater than 0.05, one at a time, until all regression coefficients were significantly different from zero.
The change-score analysis was used for changing subject i by Δ Y i = Y i1 − Y i0 , where Y i0 and Y i1 are the measured continuous pre-and posttest response variables, respectively 44 . Furthermore, the regression analysis of the mean change between two groups can be defined as follows: where Group = 0 denotes the control group and Group = 1 denotes the treatment group; mathematically, ∝ − = ∝ ⁎ 1 2 2 . X i1 , X i2 , ..., X ik are the other covariates that can affect the mean of the continuous response variable and  i is the random error. As a general rule, when Δ Y i = Y i1 − Y i0 is specified as the response variable, the baseline response Y i0 must be introduced on the right-hand side of the regression equation.
The R 2 statistic (0 ≤ R 2 ≤ 1) for a linear regression model represents the correlation between the observed and predicted response values and indicates how much of the response variability is explained by the covariates in the linear regression model.
Finally, generalized additive models (GAMs) 45 were used to increase the prediction accuracy. The GAMs were fit to detect the nonlinear effects of continuous covariates and identify appropriate cutoff points for discretizing a continuous covariate. Thus, we obtained unbiased estimates of the covariates' effects and more accurate response predictions.