Introducing BisQ, A Bicoherence-Based Nonlinear Index to Explore the Heart Rhythm

Nonlinear frequency coupling is assessed with bispectral measures, such as bicoherence. In this study, BisQ, a new bicoherence-derived index, is proposed for assessing nonlinear processes in cardiac regulation. To find BisQ, 110 ten-minute ECG traces obtained from 55 participants were initially studied. Via bispectral analysis, a bicoherence matrix (BC) was obtained from each trace (0.06 to 1.8 Hz with a resolution of 0.01 Hz). Each frequency pair in BC was tested for correlation with the HRV recurrent quantification analysis (RQA) index Lmean, obtained from tachograms from the same ECG trace. BisQ is the result of adding BC values corresponding to the three frequency pairs exhibiting the highest correlation with Lmean. BisQ values were estimated for different groups of subjects: healthy persons, persons with arrhythmia, persons with epilepsy, and preterm neonates. ECG traces from persons with arrhythmia showed no significant differences in BisQ values respect to healthy persons, while persons with epilepsy and neonates showed higher BisQ values (p < 0.05; Mann-Whitney U-test). BisQ reflects nonlinear interactions at the level of sinus-and atrial-ventricular nodes, and most likely cardiorespiratory coupling as well. We expect that BisQ will allow for further exploration of cardiac nonlinear dynamics, complementing available HRV indices.

Among nonlinear HRV indices, those derived from recurrent quantification analysis (RQA) have been particularly appealing. It has been suggested that RQA-derived HRV features are even capable of catching "the brain's ability to subtly and appropriately respond, both cognitively and emotionally, to minor changes in environmental demands" [3].
Cardiovascular signals seem to emerge from an intricate combination of periodic, quasi-periodic, linear and nonlinear, Gaussian and non-Gaussian, stochastic and deterministic interactions [4][5][6]. Accordingly, the implementation of methods based on considering some of these interactions while ignoring the contribution of others can lead to deceptive interpretations. One example might be the observed high association between HRV indices apparently related to completely different dynamical underlying mechanisms. Thus, the standard deviation along the line-of-identity of the Poincaré plot (denoted by SD2), which was regarded as a novel nonlinear feature corresponding to long-term heart traces. Regarding physiological interpretation, the meaning of the three frequency bands typically present in the HRV spectrum is far from straightforward, whereas some frequencies in an ECG trace can have direct meaning, as the heart rate, atrio-ventricular node rate, and respiratory modulation.
Furthermore, tachograms and traditional HRV features (including RQA) are derived primarily from ECG signals, and this makes it reasonable to compare bispectra acquired from an ECG with HRV measures that, certainly, originate from the same ECG recording.
In this research, we explore whether individual frequency pairs of a bicoherence matrix are associated to some of the RQA-derived HRV features. This initial exploration suggested for us a roadmap to find an index, named bispectral quotient (BisQ), conformed from specific elements of a bicoherence matrix that are maximally correlated with a well-known nonlinear EEG measure derived from recurrent quantification analysis. Expecting a high association between both types of variables does not seem to be meaningless. If this were the case, it would be possible to suggest that the explored nonlinear HRV measure is reflecting nonlinear interaction between certain frequency pairs of the ECG. Otherwise, in the case of a low correlation, it may be concluded that bispectral analysis of the ECG allows for the exploration of aspects of the heart rhythm's nonlinear dynamics that have not been reflected in existing HRV nonlinear analysis, paving the way for new promising clinical explorations.
The initial step of this study was to unravel the relationship between components of an ECG bicoherence matrix and relevant RQA nonlinear HRV features. From the outcome of this first stage, a new nonlinear index for assessing nonlinear cardiovascular dynamics was proposed (BisQ). Further, the relationship between BisQ and other HRV indices was explored, and later we attempted a statistical analysis with the proposed index as it changes with different heart/brain conditions.

Materials and Methods
This study was divided into three stages: Finding the correlation between individual bicoherence Matrix components and RQA-derived nonlinear components in a sample of ECG traces; Introduction of the bispectral quotient (BisQ) feature and its estimation; and the exploration of BisQ among different groups of ECG recordings.

Relation Between Bicoherence Matrix Components and RQA-Derived Features
The first stage of our investigation focused on the study of the statistical association between individual elements of the ECG bicoherence matrix and specific RQA-derived HRV features.

The Data
To find the relationship between individual components of a bicoherence matrix and RQA-derived HRV features, we explored 110 data ECG traces recorded from 55 participants (two ECG traces from each participant). ECG recording duration was at least 20 min and from each individual recording, two 10-min non-overlapping ECG traces were taken.
Participants were 55 persons (20-69 years of age and both genders) who did not refer any cardiovascular condition, carried a normal lifestyle, and agreed to join the study.
From each ECG trace, HRV features (from RR interval tachograms originating from the corresponding ECG trace) and bispectral/bicoherence matrices were obtained.
Care was taken for the sample to be as heterogeneous as possible, hoping that this could allow the conclusions of correlation analysis to be based on an ample range of HRV values.
Participation in the assessment process was voluntary; all subjects signed an informed consent form prior to the recording. The study was approved by the Ethics Committee of the Cuban Center for Neurosciences.
For ECG recording, a MEDICID 5 (Neuronic SA, Cuba) neurophysiological recording system was used to acquire ECG data at a sampling rate of 200 Hz. Surface electrodes were positioned diagonally across the heart according to a standard Lead II configuration. Each of the 110 ten-minute ECG traces were stored as ASCII file (.txt) and submitted to further analysis. Data were resampled to 250 Hz for Math. Comput. Appl. 2020, 25, 45 4 of 15 meeting recommended standards as well for uniformity of all data sets included in this research [25,26]. For data resampling, the resample Matlab command was used.

HRV Analysis
ECG traces were converted to tachograms (RR series) via automatic QRS detection using the publicly available package Kubios 2.1, developed at the Physics Department of the University of Eastern Finland [2]. Under visual control, ectopic and missing r-peaks were corrected by using the correction algorithm included in the Kubios package.
For the estimation of time and frequency domain HRV indices, the specific MHRV documentation master file, created by George B. Moody (1998) was downloaded from the Physionet (site at http: //www.physionet.org/physiotools/). Sixteen time-domain and frequency-domain HRV features were estimated, and Table 1 summarizes the HRV features included.  For RQA analysis of tachograms, we used the Matlab source code developed by Ouyang, downloaded from the site https://www-mathworks.com/matlabcentral/fileexchange/46765-recurrence -quantification-analysis-rqa.
The following indices were estimated: recurrence rate (RcRt), determinism (DET), Shannon entropy (Entr) and, L mean . Details about the procedure for recursive quantification analysis can be found in [2][3][4][5][6]. Briefly, in RQA, the following vectors were created: where m is the embedding dimension and τ the embedding lag.
τ] matrix of zeros and ones. The element in the j-th row and k-th column of the RP matrix, i.e., RP(j, k), is 1 if the point u j on the trajectory is close enough to point u k , that is: where d(u j , u k ) is the Euclidean distance and r is a fixed distance threshold. The structure of the RP matrix usually shows short line segments of ones parallel to the main diagonal. The embedding dimension and lag were selected to be m = 10 and τ = 1, respectively. The threshold distance r was selected to be √ SD, where SD is the standard deviation of the intervals' time series. The recurrence rate (RcRt) is defined as the ratio of ones and zeros in the RP matrix.
Accordingly, the average diagonal line length (L mean ), on the other hand, is obtained as where N L is the number of length L lines. L mean has been shown to be strongly associated to the largest Lyapunov exponent [27].

Bispectral Analysis
In the most general form, the cross-bispectrum is defined as: where X, Y, and Z correspond to three signals.
In particular, the auto-bispectrum corresponds to the case when only one signal is considered,

Bicoherence
To get a measure of the inter-frequency phase locking independent of the strength of each signal, a normalized bispectrum is identified as bicoherence. The absolute value of the bicoherence is, or rather should be, bounded by zero and one which gives a quantitative measure of the coupling. An absolute value of one infers a perfect coupling and an absolute value of zero indicates that no coupling can be observed with this measure.
Here, we use a normalization constructed by 3-norms: Then, bicoherence is defined as This normalization was proposed in [28], and it differs from the classical normalization, which is based on 2-norms (i.e., power), rather than 3-norms. With the normalization used here, the absolute value of bicoherence is bounded by 1 and the normalization is still written as a product of 3 scales, one for each of the three signals. We prefer this normalization for conceptual reasons and note that in [28], it was found that this choice of normalization has negligible impact on the statistical properties of bicoherence. For bicoherence estimation, we used the Matlab function data2bs_univar developed by Guido Nolte's Lab, which is freely available within the METH toolbox available at https://www.uke.de/english/departments-institutes/institutes/neurophysiology-and-pathophysio logy/research/research-groups/index.html.
To calculate the bispectrum, the length of segments defines the frequency resolution. Thus, if the segment has a physical length of 100 s, the frequency resolution is ∆f = 0.01 Hz.
The average is made over all segments and epochs, which would be reasonable for continuous data where the entire data is considered as one epoch. Thus, in our calculations, the sampling rate is 250 Hz. The parameters are set to segment length is 25,000 data points (100 s), segment shift 12,500 for 50% overlapping, and epoch length is 25,000 (100 s). Because the segment length corresponds to 100 s, the frequency resolution is 0.01 Hz. Since we picked the maximum number of frequency bins to be 201, we obtain as a result analyzed frequencies spanning from 0 to 2 Hz, with a spectral resolution of 0.01 Hz. Further analysis included a frequency range bounded between 0.06 and 1.80 Hz.

Correlation Analysis
For each frequency pair at the bicoherence matrix BC(f1,f2) correlation with HRV L mean was estimated over all the 110 10-min non-overlapping ECG traces. Frequency pairs were arranged according to the obtained absolute value of correlation between B(i,j) and L mean-The choice of L mean was decided after finding that this RQA feature exhibit the highest correlation with individual Bicoherence Matrix elements. Frequency pairs were analyzed in the range 0.06-1.8 Hz. The lower limit is related to the filtering properties of the ECG amplifiers, and the higher frequency limit was selected for typical values of heart rate to be included into the analysis.

Bispectral Quotien (BisQ) Estimation
BisQ is obtained as the sum of the Bicoherence values estimated at the three frequency pairs yielding the highest correlation with L mean.
where indices a, b, and c denote the three frequency pairs exhibiting the highest correlation between bicoherence and L mean . Figure 1 illustrates the flow chart followed to obtain BisQ.

Correlation Analysis
For each frequency pair at the bicoherence matrix BC(f1,f2) correlation with HRV Lmean was estimated over all the 110 10-min non-overlapping ECG traces. Frequency pairs were arranged according to the obtained absolute value of correlation between B(i,j) and Lmean-The choice of Lmean was decided after finding that this RQA feature exhibit the highest correlation with individual Bicoherence Matrix elements. Frequency pairs were analyzed in the range 0.06-1.8 Hz. The lower limit is related to the filtering properties of the ECG amplifiers, and the higher frequency limit was selected for typical values of heart rate to be included into the analysis.

Bispectral Quotien (BisQ) Estimation
BisQ is obtained as the sum of the Bicoherence values estimated at the three frequency pairs yielding the highest correlation with Lmean. BisQ = BC(f1a, f2a) + BC(f1b, f2b) + BC(f1c, f2c), where indices a, b, and c denote the three frequency pairs exhibiting the highest correlation between bicoherence and Lmean. Figure 1 illustrates the flow chart followed to obtain BisQ.

Evaluation of BisQ among Different Groups of ECG Recordings
BisQ was evaluated in four groups of ECG recordings. Statistical comparisons were performed for finding possible differences between groups. Comparisons were performed using the Mann Whitney nonparametric test because data were expected to have a non-Gaussian distribution. All data were described using mean (± SD) values. A value of p < = 0.05 should be considered statistically significant. The data sets studied correspond to different moments and settings. This study was not

Evaluation of BisQ among Different Groups of ECG Recordings
BisQ was evaluated in four groups of ECG recordings. Statistical comparisons were performed for finding possible differences between groups. Comparisons were performed using the Mann Whitney nonparametric test because data were expected to have a non-Gaussian distribution. All data were described using mean (± SD) values. A value of p < = 0.05 should be considered statistically significant. The data sets studied correspond to different moments and settings. This study was not designed to find normative values or to demonstrate specific changes related to certain conditions. It corresponds to the first evaluation with BISQ, and obtained results are aimed solely at encouraging further specific studies.

Data
We submitted to analysis four groups of ECG recordings. For all recordings, total duration was 30 min and sampling frequency was changed to 250 Hz (if needed) for uniformity and allowing comparisons, and individual values to be analyzed were obtained as the average of three BisQ estimations from non-overlapping 10-min ECG traces.
The following recordings were downloaded from the Physionet.org database at www.physionet.org.
• Healthy (20 persons): Included 20 clinically healthy participants from the Fantasia data base. Recordings were taken in resting condition. Age ranged from 21 to 81 years Details about the database can be found in [29]. • Arrhythmia (21 persons): Data were downloaded from the MIT-BIH Arrhythmia database at www.physionet.org. Age ranged from 24 to 89 years. This database is one of the most exhaustively studied in literature [10,[12][13][14][15][16][17]. It corresponds to fully annotated half-hour two leads ECGs.

•
Pre-Term Infants (7): Included seven preterm infants in the first year of life, also available on the Physionet site. This database contains simultaneous ECG and respiratory recordings from 10 preterm infants from the Neonatal Intensive Care Unit (NICU) of the University of Massachusetts Memorial Healthcare. The preterm infants were studied at conceptional ages between 29 and 35 weeks.
A fourth group included 30-min ECG recordings from nine persons suffering from focal pharmaco-resistant epilepsy. Participants' age ranged from 16 to 39 years; a MEDICID 5 (Neuronic SA, Cuba) neurophysiological recording system was used to acquire ECG data at a sampling rate of 200 Hz and resampled to 250 Hz. Surface electrodes were positioned diagonally across the heart according to a standard Lead II configuration.
Participation in the assessment process was voluntary; all subjects signed an informed consent form prior to the recording. The study was approved by the Ethics Committee of the Cuban Center for Neurosciences (Resolution No. 7/2018). Table 2 summarizes the age/gender distribution of participants.

HRV Recurrent Quantification Analysis (RQA)
In Figure 2, the recurrence plot obtained for a 10-min RR trace is represented. RQA indices corresponding to this RP are described. From each of the 110 traces described in Section 2.1.1, RQA HRV indices were obtained.

Bicoherence Analysis
In Figure 3, a plot of a typical (Log) Bicoherence matrix is shown. As expected, the highest values are associated to frequencies close to the heart rate (1.16 Hz), which is regarded as the main frequency of the heart rhythm (70 beats per minute).

Bicoherence Analysis
In Figure 3, a plot of a typical (Log) Bicoherence matrix is shown. As expected, the highest values are associated to frequencies close to the heart rate (1.16 Hz), which is regarded as the main frequency of the heart rhythm (70 beats per minute).
Frequencies varied from zero to 1. We also explored real and imaginary part of bicoherence, as well as a wider range of frequency values. At this stage, we concentrated our research on real valued bicoherence with frequencies up to 1.8 Hz.

Bicoherence Analysis
In Figure 3, a plot of a typical (Log) Bicoherence matrix is shown. As expected, the highest values are associated to frequencies close to the heart rate (1.16 Hz), which is regarded as the main frequency of the heart rhythm (70 beats per minute).

Correlation Analysis
Each frequency pair of matrix BC was tested for association strength with each RQA HRV feature. The highest correlation was found with L mean . From the viewpoint of nonlinear systems theory, the strong association between L mean and the largest Lyapunov exponent [27] might provide a sound interpretation to this feature. All these three correlation values were highly significant (p < 0.00001).

Introducing BisQ
BisQ was proposed as the sum of the three following bicoherence values: BisQ exhibited a highly significant correlation with L mean is shown. The obtained correlation coefficient r = −0.5143 was obtained.

BisQ Correlation with Other HRV Variables
Other HRV variables also exhibited significant correlations with BisQ. Table 3 summarizes the significant correlation values obtained between BisQ and those HRV variables. From Table 3, it can be noted that BisQ exhibits significant correlation with the three RQA HRV features studied, L mean being the strongest.

Principal Component Analysis (PCA)
Since HRV variables are highly correlated, we attempted to perform PCA on HRV variables and see whether uncorrelated principal components might show the stronger association with BisQ. These results are summarized in Table 4. Principal component analysis revealed that more than 99% of the variance corresponds to PC1. Only four principal components exhibited significant correlations with BisQ, but none of them were higher than those for any of the RQA indices Table 5.

Multivariate Regression
We attempted to model BisQ as a linear combination of several HRV variables. When variables were incorporated and piecewise forward regression was explored, only the variable L mean was picked as the sole predictor. We attempted to apply multivariate regression with principal components, as predicting variables did not yield a better correlation than that obtained with L mean .

Correlation with Age
In the group of healthy persons (H), ages ranged from 21 to 81 years. In the arrhythmia group (A), subject ages ranged from 24 to 89 y, and among persons with epilepsy, ages ranged from 18 to 39 years. (see Table 2). Age diversity could hinder comparisons between groups if this factor appears as a latent variable. For this reason, we tried to clarify whether a correlation between person´s age and BisQ is present.
In none of the groups a significant correlation for BisQ with age was obtained. Nonetheless, when we found significant differences between healthy persons and persons with epilepsy, we again compared the group of persons with epilepsy with the subset of a young healthy patient, and still found differences, notwithstanding the lower number of individuals (see Table 2).

Comparing Groups
In Figure 4, a Box and Whisker Plot summarizes the values of BisQ in the four groups studied. As it can be observed from Figure 4, BisQ seems to be higher among persons with epilepsy as well as preterm infants, whereas persons with arrhythmia have similar mean and a higher dispersion as compared to healthy persons.  Comparisons between groups with the Mann-Whitney´s nonparametric test Table 5 revealed that the group of persons with epilepsy exhibited higher BisQ values than the group healthy persons as well as persons with arrhythmia. Pre-term infants group also had higher BisQ values than the healthy persons group.  Comparisons between groups with the Mann-Whitney´s nonparametric test Table 5 revealed that the group of persons with epilepsy exhibited higher BisQ values than the group healthy persons as well as persons with arrhythmia. Pre-term infants group also had higher BisQ values than the healthy persons group.
Comparing young healthy individuals with epileptic patients. Since the group of persons with epilepsy are in the 16 to 39 years age group, and the age range for healthy persons group is 21 to 81 years, the observed difference in BisQ values could be a result of age mismatch. Thus, we compared the group of epileptic patients with the subgroup of healthy young patients. The following BisQ values were obtained: 1.70 ± 0.45 for healthy young patients. 2.00 ± 0.13 for epileptic patients.
The Mann Whitney test still showed significant difference between groups (p < 0.05)

Discussion
In this study, an index derived from bispectral analysis has been introduced. As obtained, the Bispectral Quotient BisQ accounts for around 27% of the total variance of L mean , an RQA variable for nonlinear HRV analysis that has been regarded as being strongly associated to the maximum Lyapunov exponent [27]. In this regard, we interpret this result as evidence that this measure, derived from bispectral analysis is at least partially related to other nonlinear HRV measures that have been reported as parsimonious for detecting shrewd changes in cardiac activity as well as regulation from other physiological and environmental influences [3]. At the same time, the fact that at least 70% of BisQ variance is not related to any other HRV feature suggests that BisQ retains relevant information about heart rhythm nonlinear dynamics that could complement commonly used HRV features. Thus, we expect that, combined with traditional linear and nonlinear HRV measures, BisQ can become a useful tool for exploring heart rhythm and its modification with both environmental and physiological influences.
Exploration of BisQ with different groups showed that BisQ is apparently higher among persons with epilepsy as well as among preterm infants.
With regard to the physiological meaning of obtained results, unlike other indices, BisQ was obtained as result of an attempt to find an association between bicoherence and other previously established nonlinear HRV measures.
As can be seen from its definition (10) Finally, the pair (1.05 Hz, 0.33 Hz) might be related to nonlinear interactions related to cardiorespiratory coupling. Respiratory signals usually display prominent spectral peaks at the first harmonic, at about 0.34 Hz. As it is well known, respiration is an important conditioner of heart rate, and nonlinearities in this important mechanism might be of special interest.
Thus, we assume the idea that the pair (1.05 Hz, 0.33 Hz) reflects nonlinear phenomena associated to cardiorespiratory coupling.
Worth of note, are the results obtained by Toledo et al. [19] with HRV bispectra after heart transplantation. They found few (about seven) significant phase coupled spots in the 0 to 1.5 Hz frequency range. Curiously, two of them were coincident with the frequency pairs and (1.19 Hz, 1.19 Hz), and (1.05 Hz, 0.33 Hz) included into BisQ. This could support the idea that the coupled frequencies that have been contributing to BisQ reflect mainly intrinsic heart regulatory mechanisms with a lesser influence from autonomic cardiovascular regulation.
Our results might point to an increase in BisQ in inter-ictal epilepsy. Increases in heart rate during an epileptic seizure is a well-known fact [30], but during inter-ictal periods, heart rate among epileptic patients is similar to that of healthy individuals. In our data base of epileptic patients, we found that inter-ictal heart rate ranges between 64 and 80 beats per minute (for a mean frequency of 1.20 Hz). This suggests that the observed changes are not associated to changes in heart rate in this group of patients.
At the same time, attempts to characterize the inter-ictal heart rhythm of persons with epilepsy is an ongoing task [27]. In this sense, BisQ, combined with other nonlinear HRV features could be a promising factor in epilepsy research and treatment.
On the other hand, our results deserve further clarification. As it has been reported, during epileptic seizures, neurovascular coupling is impaired [31,32]. However, these changes during inter-ictal periods do not appear in all types of epilepsies. Whether or not increased inter-ictal BisQ among focal drug-resistant patients with epilepsy might be a compensatory mechanism deserves a dedicated research.
In epilepsy, excessive synchronization has been reported as a hallmark for this condition. It is not excluded that brain synchronization might lead to over-synchrony in main brain rhythms. During, before, and after an epileptic seizure, cardiorespiratory coupling has been reported as reduced, but changes seem to depend on the type of epilepsy [23]. Since BisQ could combine both cardiac pacemaker mechanisms and cardiorespiratory coupling, a special study is required for assessing individual contributions corresponding to each frequency pair included into this index.
No differences in BisQ values were observed among persons with arrhythmia, a counter-intuitive result in view of impaired synchronization mechanisms among this group of persons [33,34]. However, the relatively low proportion of arrhythmic events in long-lasting ECG recordings could explain the absence of significant BisQ changes in this group of persons with arrhythmia.
Interpretation of results with preterm infants goes beyond the scope of this study. Heart rate and respiratory rate among newborns nearly double those of adults, and BisQ may not be interpreted in terms of corresponding frequencies as it was done with adult groups. It is not excluded that a tailored BisQ version for infants should be introduced when infants will be studied.
Biophysical mechanisms of cardiorespiratory coupling, though widely known, are scarcely clarified [35][36][37]. We hope that BisQ, an index that specifically includes frequencies putatively involved in this phenomenon, might contribute to cardiorespiratory coupling research.
Limitations of present research are mainly related to the low number of cases included into the study. The fact that in this exploratory study groups recorded in different conditions and even with different recording systems were compared, limits the clinical interpretation of our results: they can be interpreted as preliminary and guiding to further research. Thus, the validity of these results is merely exploratory to check for the plausibility of using BisQ with real data and to encourage properly designed study protocols. Factors as intra-individual variability and differentiation of epileptic patients from healthy individuals, as well relation to physical training will be addressed in future research from our group.

Conclusions
In this study, we proposed a bicoherence-derived index that exhibits a significant correlation with nonlinear heart rate variability indices. No differences were documented between normal adults of different age. Departures from normal values were obtained among a group of persons with epilepsy. Apparently, this index reflects nonlinear interactions at the level of sinus-and atrial-ventricular nodes, and most likely cardiorespiratory coupling. Hopefully, BisQ will allow further exploring cardiac nonlinear dynamics, complementing available HRV indices