Do meditators and non-meditators have different HRV dynamics?

Available online at www.sciencedirect.com

ScienceDirect Cognitive Systems Research 54 (2019) 21–36 www.elsevier.com/locate/cogsys

Do meditators and non-meditators have different HRV dynamics? Action editor: Luis H. Favela Atefeh Goshvarpour a, Ateke Goshvarpour b,⇑ a

Department of Biomedical Engineering, Faculty of Electrical Engineering, Sahand University of Technology, Tabriz, Iran b Department of Biomedical Engineering, Imam Reza International University, Mashhad, Razavi Khorasan, Iran Received 12 August 2018; received in revised form 22 October 2018; accepted 22 November 2018 Available online 27 November 2018

Abstract Purpose: The heart is a complex system and many researchers have been recently studying cardiac behavior using the theory of nonlinear dynamical systems. One of the most appealing tools for analyzing heart function is the heart rate variability (HRV) signal. This study aimed to elucidate the HRV dynamics of six distinct states: spontaneous normal breathing (SNB) and metronomic breathing (MB), as non-meditator groups, before Chinese Chi meditation (CCM), during CCM, before Kundalini yoga meditation (KYM), and during KYM, as meditator groups. Methods: The HRV data were obtained from the Physionet database. Lagged Poincare indices, Lyapunov exponent (LE), Lempel-Ziv complexity (LZ), and 4 types of entropy were calculated. Results: The results showed the greatest discrepancies in the lagged Poincare indices for KYM and CCM. In contrast, the least variations were achieved for MB. Compared to SNB, an enhancement in the log energy entropy and a reduction in the LZ and other entropies were concluded during KYM and CCM practices. In contrast, a reverse pattern was observed for MB. Using support vector machine, HRV dynamics were classified with average accuracies of 99.14 and 98.2% and average sensitivities of 99.87 and 99.57% for pre-KYM and during KYM, respectively. Conclusion: It was shown that the HRV dynamics were significantly different in meditators and non-meditators. Ó 2018 Elsevier B.V. All rights reserved.

Keywords: Nonlinear dynamics; Heart rate variability; Meditators; Non-meditators; Support vector machine

1. Introduction Abbreviations: ApEn, approximate entropy; CCM, Chinese Chi meditation; EEG, electroencephalogram; HF, high frequency; HRV, heart rate variability; kNN, k nearest neighbor; KYM, Kundalini yoga meditation; LE, Lyapunov exponent; LF, low frequency; LogEn, log energy entropy; LT, long term; LZ, Lempel-Ziv complexity; MB, metronomic breathing; PSD, power spectral density; SaEn, sample entropy; ShEn, Shannon entropy; SNB, spontaneous nocturnal breathing; ST, short term; SVM, support vector machine; VLF, very low frequency ⇑ Corresponding author at: Khayyam Boulevard, Khayyam 20, Number 54, 4th Floor, PO. BOX: 91876-83883, Mashhad, Khorasan Razavi Province, Iran. E-mail addresses: [email protected] (A. Goshvarpour), [email protected] (A. Goshvarpour). https://doi.org/10.1016/j.cogsys.2018.11.010 1389-0417/Ó 2018 Elsevier B.V. All rights reserved.

Changes in the heart’s beat to beat intervals are known as heart rate variability (HRV). HRV demonstrates autonomic regulations resulting from the function of both sympathetic and parasympathetic nerves. In pathological and psychological assessments, many researchers have been interested in analyzing HRV signals. HRV has been repeatedly used in the evaluation of heart diseases, such as arrhythmias and psychological states, such as emotion and concentration (Goshvarpour, Abbasi, & Goshvarpour, 2017b; Goshvarpour, Goshvarpour, &

22

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

Rahati, 2011; Phongsuphap & Pongsupap, 2011; Reed, Robertson, & Addison, 2005). Meditation belongs to a family of mental training practices, which aims to familiarize the practitioner with definite kinds of mental processes (Brefczynski-Lewis, Lutz, Schaefer, Levinson, & Davidson, 2007). Meditation causes reduced metabolic activity that provokes physical and mental relaxation and it is commonly suggested as a stress management technique. This technique is also suggested to improve psychological and emotional stabilities (Jevning, Wallace, & Beidebach, 1992). Since meditation (unlike medication) has no side effects, it is usually recommended for non-medicated management of some psychological conditions. So far, various meditation and relaxation protocols have been studied to get insight into the health benefits of these techniques (Peng et al., 1999). Up to now, numerous methods have been presented for analyzing and classifying HRV signals. Initial studies have focused on the statistical and spectral approaches. These methods are still used in many studies. By spectral analysis of HRV, the frequency band of the signal is divided into three fundamental parts: (1) very low frequency (VLF), (2) low frequency (LF), and (3) high frequency (HF) (Task Force of the European Society of Cardiology, 1996). The effects of thermoregulation and humoral factors, baroreflex-related activities, and the respiratory sinus arrhythmia can be studied through VLF, LF, and HF, respectively (Cerutti, Bianchi, & Mainardi, 1995; van Ravenswaaij-Arts, Kollee, Hopman, Stoelinga, & van Geijn, 1993). Several studies have shown that the heart is a complex system with nonlinear interactions. The continuous interaction of the sympathetic and vagus nerves in controlling HRV dynamics is an example of this nonlinear behavior (Bai et al., 2008). Therefore, many researchers believe that the conventional methods of HRV analysis cannot provide comprehensive information. Nonlinear time-series analysis is the alternative solution. These methods apply chaos theory and nonlinear dynamics to appreciate unpredictable behavior. Using the concept of state space or phase space reconstruction, calculation of some characteristic measures is possible, which allows us to guess the future behavior of the time-series. For example, two common measures are the Lyapunov exponent (LE) and entropies. They refer to unpredictability and irregularity because of their sensitive dependence on initial conditions. As already mentioned, the HRV signal has a chaotic nature. Previously, some nonlinear methods have been studied to reveal HRV dynamics during meditation. Nonlinear features applied in meditation studies are Hilbert transform (Peng et al., 1999), dynamical complexity (Li, Hu, Zhang, & Zhang, 2011), Poincare plots (Goshvarpour et al., 2011; Goshvarpour & Goshvarpour, 2015), higher order spectra (Goshvarpour & Goshvarpour, 2013), LE (Goshvarpour and Goshvarpour, 2012a; 2012b), recurrence plots (Goshvarpour & Goshvarpour, 2012c), fractal scaling

(Alvarez-Ramirez, Rodrı´guez, & Echeverrı´a, 2017), multifractal analysis (Song, Bian, & Ma, 2013), wavelet entropy (Gao et al., 2016), and approximate entropy (ApEn) (Goshvarpour & Goshvarpour, 2012a). From the perspective of physiological signal processing, valuable information has been achieved about the benefits and effects of meditation on the body. Table 1 summarizes the previous findings. Although these studies have emphasized the physiological changes during meditation, the comparative physiological effects of different meditation techniques have not been completely investigated. Additionally, still little information is available about the effect of meditation on cardiac dynamics. The innovations of this paper are summarized as follows: 1. The previous studies (Table 1) are limited to the analysis of one or two psychological modes (meditation technique). While we analyzed and compared the physiological effect of six different modes in this study. 2. In most of the previous studies, time- and frequencydomain techniques have been used for the analysis of biological signals, and only a limited number of them have used nonlinear approaches, which are based on the nature of the chaotic behavior of the biological signals. However, in this part of the literature, a limited number of features have been extracted. For example, lagged Poincare indices (lags of 1–6) were studied in Chi meditators in Goshvarpour and Goshvarpour (2015), LE and ShEn were used while subjects focused on breathing (Goshvarpour & Goshvarpour, 2012a). In this study, we applied 46 nonlinear features. These features were selected to provide global information about the reconstructed phase space structure of HRV signal (LE, entropies, and LZ) and to describe the shape details of HRV trajectories (Poincare Indicators). 3. Most studies that have examined the effects of meditation on physiological parameters are limited to extracting some features and little study has been done on the classification. This study also included the classification process. The goal of the classification step was to examine the separability of the HRV behaviors and characteristics in the different states. Totally, the main goals of this study are (1) to elucidate nonlinear dynamics of HRV signals using entropy (sample entropy (SaEn), ApEn, Shannon entropy (ShEn), and log energy entropy (LogEn)), Lempel-Ziv complexity (LZ), LE, and lagged Poincare indices; (2) to find the difference in chaotic degree of HRV signals between meditators practiced Chinese Chi meditation (CCM), Kundalini yoga meditation (KYM), and non-meditators, including metronomic breathing (MB) and spontaneous nocturnal breathing (SNB); (3) to classify nonlinear HRV dynamics in different states.

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

23

Table 1 Previous achievements on meditation applying physiological parameters. Study

Meditation technique

Signal

Number of Methodology subjects

Main findings

Aftanas and Golocheikine (2001)

Sahaja Yoga

EEG

Theta and alpha synchronization during meditation

Takahashi et al. (2005)

Zen

11 ST- & 16 Spectral analysis LTmeditators 20 Spectral analysis

Phongsuphap, Pongsupap, Chandanamattha, and Lursinsap (2008) Sarkar and Barat (2008)

Samadhi

EEG & HRV & breath rate HRV 35 Power spectrum (PSD) meditators & 70 controls HRV 8 Scaling, complexity, and recurrence analysis

slow alpha1 and fast theta2 increased in frontal area during meditation

During meditation, the PSD of HRV shifts toward a particular frequency band to produce the resonant peak. Chi The complexity of the HRV increases and the long range correlation exhibited by the normal heart was destroyed during meditation. Goshvarpour et al. (2011) Chi HRV 8 Lagged Poincare As the lag increased, the width of Poincare plot increased during meditation. Li et al. (2011) Chi and KundaliniHRV 8 vs. 4 Dynamical complexity Dynamical complexity decreased in meditation Phongsuphap and Pongsupap Samadhi state and HRV 133 vs. 117 time - and frequency - based Classification accuracy: 94.8% (2011) non-Samadhi. HRV features and Fisher discriminant segments analysis Goshvarpour, Goshvarpour, Focused on EEG 11 Bispectrum More increments of phase coupling in Rahati, and Saadatian breathing meditators & occipital region (2012) 4 nonmeditator Goshvarpour and Focused on ECG 25 LE and entropy, Quadratic, The maximum classification accuracy was Goshvarpour (2015) breathing fisher, and kNN achieved using Quadratic with the rate of 92.31%. Goshvarpour and Chi and KundaliniHRV 8 vs. 4 Poincare plots, LE and Hurst Chaotic behavior of HRV was decreased in Goshvarpour (2012b) exponents both meditation techniques. Goshvarpour and Chi HRV 8 Recurrence plot Nonlinear interaction of variables decreased Goshvarpour (2012c) in meditation state Muralikrishnan, Isha Yoga HRV 14 Frequency and time domain Significant differences between meditators and Balakrishnan, meditators & indices controls in both frequency and time domain Balasubramanian, and 14 controls indices, with no difference in resting between Visnegarawla (2012) the groups. Peressutti, Martı´n-Gonza´lez, Zen HRV 20 Discrete wavelet transform, Different practice time of meditators reflect and Garcı´a-Manso (2012) Principal components analysis, stage-specific patterns of HRV. Symbolic dynamics Goshvarpour and Chi and KundaliniHRV 8 vs.4 Higher order spectra During both techniques phase-coupled Goshvarpour (2013) harmonics were moved to the higher frequencies. Goshvarpour, Shamsi, and Chi HRV 8 Spectral and time analysis Different ARMA model order before and Goshvarpour (2013) during meditation. Song et al. (2013) Chi HRV 8 Multi-fractal During meditation the HRV becomes regular and the degree of multifractality decreases. Goshvarpour and Chi HRV 8 Lagged Poincare Increase in SD1/SD2 and decrease in the area Goshvarpour (2015) of Poincare plots during meditation compared to that before meditation. Beutler, Beltrami, Boutellier, Yoga Respiration12 Tidal volume Yoga resulted different respiratory regulation and Spengler (2016) at rest and hypercapnia Gao et al. (2016) Mindfulness EEG & 11 Wavelet entropies Meditation increased the entrainment HRV between mind and body Jadhav, Manthalkar, and Focused on EEG 11 Spectral & Hjorth features & The positive affects has increased and the Joshi (2017) breathing kNN improvement of subject’s arousal state after 8 weeks of meditation. Alvarez-Ramirez et al. (2017) Chi and KundaliniHRV 8 vs.4 Fractal scaling Hurst exponent for meditators decreased, reflecting uncorrelated HRV dynamics Braboszcz, Cahn, Levy, Vipassana, EEG 20 vs. 27 vs., Spectral analysis All meditators showed higher parietoFernandez,and Delorme Himalayan Yoga 20 occipital gamma amplitude than control (2017) and Isha Shoonya subjects Note: EEG: Electroencephalogram; ST: Short term; LT: Long term; PSD: Power spectral density; kNN: k nearest neighbor.

24

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

2. Material and methods First, HRV signals were selected from the Physionet database (Peng et al., 1999). These signals belong to the following classes: SNB, MB, before CCM, during CCM, before KYM, and during KYM. Second, the data were segmented using a non-overlapped sliding window with the length of 180 samples. Then, the HRV segments were characterized by nonlinear features, including lagged Poincare measures, entropies, LZ, and LE. Finally, the normalized features were fed to the support vector machine (SVM) to automatically find the HRV differences between the classes. Fig. 1 describes these steps schematically. 2.1. Data selection In this study, the HRV signals of healthy meditators and non-meditators available in Physionet database (Peng et al., 1999) were analyzed. Two types of meditation, including the KYM (as trained by Yogi Bhajan) and CCM, were considered. In the former, four well-trained Kundalini yoga meditators were asked to wear a Holter monitor. Their ECG signals were recorded for about one and a half hours. Before starting the protocol, 15 min of ECG signals were acquired (pre-KYM). Then, the KYM protocol, including a sequence of breathing and chanting exercises, was done for about 1 h. The participants were asked to sit in a cross-legged posture during KYM. In the latter, HRV of 8 Chi meditators were recorded before and during CCM. All meditators were chosen from healthy graduated or post-doctoral students. Most of them started the CCM practice for about 1–3 months before the experiment. They were asked to wear a Holter monitor. In this situation, they performed their daily activities for about 10 h. At about 5 h into the ECG recording, each of them performed 1 h of CCM. During the experiment, the subjects did the task in a sitting position and listened to the taped guidance of the master. The participants were trained to breathe spontaneously while imagining the opening and closing of a perfect lotus in their stomach. To investigate the benefits of the meditation techniques (KYM and CCM) on cardiac functioning, the HRV signals of two control groups of non-meditating healthy subjects available in Physionet database (Peng et al., 1999) were also studied, including (1) a group of 11 healthy subjects which breathe spontaneously during sleep hours with dura-

tion of about 6 h, except for one subject (4.6 h) and (2) a group of 14 healthy volunteers during supine MB at a fixed rate of 0.25 Hz. The general health of all groups (meditators and nonmeditators) were comparable (Peng et al., 1999). The sampling rate was 128 Hz. More details about the groups are provided in Table 2. For initial evaluation of the data, the mean HRV signals were calculated for each class in a non-overlapped 180sample window length. The results are shown in Fig. 2. As the figure shows, the average HRV of two groups of meditation (CCM and KYM) was higher than that of SNB and MB. In addition, the HRV of the Kundalini’s meditators increased significantly during KYM compared to before KYM. 2.2. Feature extraction Nonlinear dynamics of a system can be examined by a trajectory in a multidimensional space, often mentioned as the state space or phase space. If the trajectory converges to a sub-space in the phase plane, this sub-space is entitled as the attractor of the system. In order to interpret the attractor in the phase plane, LE, entropy, and Poincare plot indices are commonly employed. LE and entropy measures can be considered as ‘dynamic’ indices of attractor complexity. Using the concepts of these measures, it is possible to get more detailed information about the chaotic behavior of HRV signals. There are many methods to investigate the chaotic behavior of a system. In this study, lagged Poincare plot indices, LE, entropies (ApEn, SaEn, ShEn, and LogEn), and LZ measures have been evaluated to distinguish the HRV signals of meditators and nonmeditators. Precisely, 4 Poincare measures, including SD1, SD2, SD1/SD2, and area were extracted from lagged Poincare plots (in 10 lags). Therefore, 40 features for lagged Poincare plots along with another 6 nonlinear indices were extracted. 2.3. Statistical analysis To specify significant differences between the extracted parameters, the Wilcoxon ranksum test (also known as Mann-Whitney U test) was performed with the null hypothesis that x and y vectors comprise independent points from a randomly selected value of one vector will

Fig. 1. Proposed methodology.

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

25

Table 2 Subject characteristics in groups of meditators and non-meditators. Meditators

Num Age Mean Gender

Non-meditators

Kundalini yoga  before KYM  during KYM (Y1_pre-Y4_pre) (Y1_med-Y4_med)

Chinese Chi meditation  before CCM  during CCM (C1_pre-C8_pre) (C1_med-C8_med)

Spontaneous nocturnal breathing

Metronomic breathing

(N1-N11)

(M1-M14)

4 20–52 33 2 women, 2 men

8 26–35 29 5 women, 3 men

11 20–35 29 8 women, 3 men

14 20–35 25 9 women, 5 men

Note- Num: number of subjects; Age: Age range; Mean: Mean age.

Applying SVM classifier, a best separating hyper-plane is chosen for the HRV feature space. This hyper-plane results in the largest possible separation margin between two classes. Assume a data set (T) with i samples: T ¼ fðx1 ; y 1 Þ; ::::; ðxi ; y i Þg i ¼ 1; 2; 3; ::::; i xi 2 Rn ;

y i 2 f1g

ð2Þ

where the input and the class labels are xi and yi, respectively. The separating line has the following formulation: hw; xi þ b ¼ 0

ð3Þ

where w 2 R , hw; xi shows the inner product of w and x, and b represents the bias. To obtain the best separating line with the highest possible margin between classes, the following criteria should be satisfied: d

Fig. 2. The mean HRV values in SNB, MB, pre-CCM, during CCM, preKYM, and during KYM. The statistical test (Wilcoxon rank sum test) was performed to compare the significance of the difference between each pair of groups (15 different conditions); red arrows indicate that there is no significant difference between the two groups. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

be higher or smaller than a randomly selected point from a second vector. The result of the Wilcoxon ranksum test is mirrored in p-value. It is a non-parametric technique, which is applied to examine the similarity of two dependent data samples with an ordinal scale. In addition, it does not implicate the normal distribution assumption (unlike the ttest). In medical studies, when the use of conventional methods (parametric techniques) based on normal distribution does not fulfill or there is not any information about the distribution of the data, non-parametric techniques are implemented. To evaluate the substantive significance, the Cohen’s effect size (Cohen, 1988) was calculated. To compute this parameter, the difference between the group means is divided by the pool standard deviation.

ðhw; xi þ bÞ P 1

1 As kwk ¼ kwk is the interval between the closest data point and separating line, the margin is 2=k w k. hw; wi=2 should be minimized to maximize the margin. Therefore, instead of investigating the best separating line, the optimization criteria which is formulated as Eq. (5) should be fulfilled:

min 12 k w k

2

w;b

s:t: zi ðwt xi þ w0 Þ P 1

The classification of HRV signals between SNB and each of MB, pre-CCM, during CCM, pre-KYM, and during KYM states was implemented using SVM. Before feeding the extracted features (F) to the SVM, they were normalized (Norm_F) according to Eq. (1):
 F  F min NormF ¼ 2 1 ð1Þ F max  F min

8i

ð5Þ

By solving the above mentioned equation, the sample set is separated into two classes. Since samples are commonly not divided completely, the SVM employed kernel strategy which is formulated as follows: X f ðxÞ ¼ sgnð ai y i K ðxi :xÞ þ bÞ ð6Þ i¼1

2.4. Classification

ð4Þ

jwt xi þbj

  Pm where b ¼ y i  i¼1 yi ai xi :xj , ai depicts Lagrange multipliers, and Kðxi ; xj Þ is a chosen kernel function (Cortes & Vapnik, 1995). More details about SVM classifier can be found in (Cortes & Vapnik, 1995; Hu and Hwang, 2002). In this study, the ‘‘Gaussian radial basis function kernel” was carried out and the ‘‘Least square” method was applied to find the separating hyperplane. To evaluate the classification performances, accuracy (Acc) and sensitivity (Sen) were calculated in this study.

26

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

3. Results Several nonlinear features, including lagged Poincare measures, entropies, LZ, and LE were extracted from the HRV signals. Fig. 3 demonstrates the variability of the lagged Poincare measures in different states (top frames). Exactly, each box indicates the parameter variability in 10 lags (lag 1–10). To investigate the group differences, a specific lag (lag 10) was selected and the direct results of Poincare parameters in that lag are also visualized in Fig. 3 (bottom frames). The most variations in the lagged Poincare indices occur during KYM (Fig. 3, top frames). In contrast, the least variation is observed for MB. In addition, the results show that not only there are significant differences between the lagged Poincare indices in two groups of meditators (both before and during CCM and KYM) but significant differences are also observed between the HRV features of the meditators and SNB. Considering a particular lag (Fig. 3, bottom frames), the pattern of the Poincare parameters varies in different classes. As the figure shows, the highest amount of all Poincare measures is obtained for KYM. In addition, it can be easily understood that compared with non-meditator groups, the value of Poincare indices is higher for meditator subjects. Fig. 4 depicts the value of entropies, LZ, and LE in different states. In general, significant differences are perceived between the HRV measures of different states (Fig. 4). Compared to SNB, an increase in the logEn and a decrease in the LZ and other entropies are observed for two types of meditation. In contrast, an opposite pattern is perceived for MB. There is a drop in the logEn and a rise in the LZ and other entropies.

To compare the HRV measures of different states, the mean and standard deviation of the nonlinear characteristics have been provided in Table 3. Table 3 and Fig. 4 show that there are different patterns between nonlinear HRV features of different states. For example, compared to SNB, the mean value of ApEn decreases for CCM and KYM. A similar pattern is observed for SaEn and ShEn. Comparing MB with two classes of meditation (CCM and KYM) and SNB, there are different HRV dynamics. LE, LZ, ShEn, and SaEn are significantly higher during MB than those of the other five classes. Considering lagged Poincare measures, valuable information about dynamical attractors of the HRVs has been obtained for different states. As the lag increases, an increase in SD1 and a decrease in SD2 are observed for all states (Table 3). However, increasing the lags, the greatest variations are perceived for KYM, in which the highest increase in SD1 and the highest decrease in SD2 are detected. The second largest changes in the phase space are observed for CCM. Fluctuations in SD1 for lags of 1–10 are in the range of 2.44–6.43, while the corresponding SD2 changes from 4.8 to 4.33. In contrast, the lowest variations in the lagged Poincare measures are perceived for MB (Fig. 3 and Table 3). Subsequently, similar patterns are perceived for the other features extracted from the lagged Poincare indices (SD1/ SD2 and area). For example, considering the area, the data are less spread out in the phase space during MB, while data points are more dispersed for KYM. During the CCM and KYM practices, the lowest SD1/SD2 values in lag 1 are recognized compared to the other 4 groups. While, in lag 10, the highest SD1/SD2 values are obtained for these two classes. In other words, the greatest fluctuations in the SD1/SD2 are realized for CCM and KYM.

Fig. 3. Top: The variability of lagged Poincare measures in 10 different lags for all states. Bottom: lagged Poincare measures in lag 10 for all states. (a) SD1 and SD2; (b) SD1/SD2; (c) area.

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

27

Fig. 4. The Boxplot of nonlinear features for SNB, pre-CCM, during CCM, pre-KYM, during KYM states. (a) LE; (b) LZ; (c) LogEn; (d) ShEn; (e) SaEn; (f) ApEn. Red arrows indicate that there is no significant difference between the two groups. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To determine significant differences between the parameters, Wilcoxon rank sum test (Mann-Whitney U test) was performed. The p-values with their corresponding effect size (Cohen, 1988) are also listed in Table 3. Precisely, Table 3 shows the statistical results between SNB and other groups. For an easier and better interpretation of the results, the percentage of significant indices were computed. Table 4 briefly designates the results for the Poincare plot features. In addition, the non-significant differences between each pair of groups are shown with red arrows in Fig. 4. Totally, 189 out of 230 features show significant differences (p < 0.05) between SNB and the other classes. It can be concluded that there are significant differences between SNB and other groups for most of the nonlinear features. Considering all features, LZ, LogEn, ShEn, and

ApEn nominate as the best nonlinear features for showing the differences between the groups. However, the largest effect size is achieved for MB. In the next stage, dynamical HRV patterns in different states were classified using SVM. A binary SVM classifier was implemented to separate SNB from the other 5 groups. To this end, a 5-fold cross-validation approach was adopted. In this manner, data is separated into 5 equal portions: one for the test and the remaining for the training. The cross-validation was repeated for 50 times. The classification accuracies are shown in Fig. 5. The horizontal axis indicates the number of classification run (50 times), and the vertical axis represents the classification accuracies. As we can see, SVM gives the highest classification rates for pre-KYM and during KYM. In these cases, the mean accuracy rates are higher than 98%. For a closer look,

SNB mean ± std

28

Table 3 Mean and standard deviation (mean ± std) of the extracted features from HRV signals in different groups. Statistically significant indices between SNB and other groups are shown with: values < 103 and * p-values < 0.05. The effect size (ES) of the extracted HRV features between SNB and the other states are also inserted.

**

p-

ES

Before CCM mean ± std

ES

CCM mean ± std

ES

Before KYM mean ± std

ES

KYM mean ± std

ES

Entropy features ApEn 0.71 ± 0.09 SaEn 1.50 ± 0.46 ShEn (0.67 ± 0.23)  107 LogEn (1.50 ± 0.06)  103 LZ 3.88 ± 0.53 LE 0.86 ± 0.13

0.71 ± 0.04** 1.61 ± 0.26** (0.64 ± 0.18)  107** (1.50 ± 0.05)  103** 4.39 ± 0.32** 0.91 ± 0.07**

2.00 5.15 5.04 5.24 5.23 0.58

0.70 ± 0.07** 1.36 ± 0.42** (1.1 ± 0.34)  107** (1.58 ± 0.05)  103** 3.71 ± 0.33** 0.79 ± 0.11

1.82 2.73 2.80 2.96 2.94 0.09

0.68 ± 0.06** 1.24 ± 0.34** (1.01 ± 0.32)  107** (1.57 ± 0.05)  103** 3.73 ± 0.38** 0.86 ± 0.08**

1.77 2.42 2.86 3.01 2.96 0.45

0.74 ± 0.03* 1.50 ± 0.27 (0.61 ± 0.15)  107** (1.49 ± 0.05)  103** 3.97 ± 0.37* 0.83 ± 0.11

0.31 0.002 0.26 0.22 0.17 0.23

0.59 ± 0.03** 0.84 ± 0.10** (1.43 ± 0.27)  107 (1.62 ± 0.03)  103 3.88 ± 0.34** 0.88 ± 0.06

0.60 0.63 1.06 0.73 0.25 0.19

Lagged poincare indices 2.54 ± 1.21 SD1(1) SD1(2) 3.48 ± 1.74 SD1(3) 3.58 ± 1.80 SD1(4) 3.56 ± 1.69 3.89 ± 1.66 SD1(5) SD1(6) 4.18 ± 1.82 SD1(7) 4.27 ± 2.00 SD1(8) 4.31 ± 2.00 4.35 ± 1.93 SD1(9) SD1(10) 4.36 ± 1.95

2.30 ± 0.88** 3.24 ± 1.15** 3.17 ± 1.04** 2.88 ± 0.86** 3.10 ± 0.97** 3.43 ± 1.11** 3.44 ± 1.06** 3.23 ± 0.99** 3.13 ± 1.00** 3.24 ± 1.07**

4.89 4.94 4.98 5.02 5.00 4.98 5.01 5.01 4.99 4.98

2.65 ± 0.97** 3.31 ± 0.86** 3.73 ± 1.05** 4.01 ± 1.15** 4.35 ± 1.22** 4.61 ± 1.28** 4.79 ± 1.34** 4.88 ± 1.40** 4.94 ± 1.44** 5.02 ± 1.46**

2.75 2.85 2.83 2.83 2.83 2.84 2.83 2.83 2.82 2.82

2.44 ± 1.07** 3.47 ± 0.93** 4.34 ± 1.24** 4.97 ± 1.52** 5.51 ± 1.77** 5.93 ± 1.98** 6.19 ± 2.17** 6.31 ± 2.3** 6.43 ± 2.39** 6.43 ± 2.43**

2.65 2.85 2.86 2.85 2.84 2.83 2.81 2.79 2.79 2.78

2.16 ± 0.35 2.83 ± 0.45 2.96 ± 0.66 3.36 ± 0.81* 3.79 ± 0.77 3.82 ± 0.76 3.81 ± 0.89 3.99 ± 0.85 4.04 ± 0.81 3.98 ± 0.99

0.31 0.38 0.35 0.12 0.06 0.20 0.23 0.16 0.16 0.20

4.09 ± 2.29** 6.20 ± 2.50** 8.15 ± 2.78** 9.71 ± 3.00** 10.94 ± 3.17** 11.73 ± 3.23** 12.14 ± 3.22** 12.23 ± 3.28** 11.92 ± 3.26** 11.45 ± 3.45**

0.56 0.64 1.06 1.69 1.87 1.72 1.70 1.72 1.66 1.62

SD2(1) SD2(2) SD2(3) SD2(4) SD2(5) SD2(6) SD2(7) SD2(8) SD2(9) SD2(10)

6.77 ± 4.13 6.22 ± 4.13 6.22 ± 4.01 6.30 ± 3.96 6.08 ± 4.00 5.87 ± 3.95 5.81 ± 3.86 5.80 ± 3.82 5.79 ± 3.84 5.79 ± 3.81

4.86 ± 1.65** 4.24 ± 1.63** 4.35 ± 1.54** 4.59 ± 1.54** 4.43 ± 1.50** 4.17 ± 1.50** 4.15 ± 1.49** 4.33 ± 1.50** 4.43 ± 1.42** 4.33 ± 1.44**

4.96 4.89 4.94 4.97 4.96 4.95 4.93 4.95 4.99 4.97

9.11 ± 4.34** 8.87 ± 4.41** 8.69 ± 4.40** 8.54 ± 4.40** 8.36 ± 4.42** 8.21 ± 4.41** 8.10 ± 4.40** 8.05 ± 4.38** 8.01 ± 4.36** 7.97 ± 4.35**

2.63 2.61 2.60 2.59 2.57 2.56 2.56 2.56 2.56 2.55

8.47 ± 2.69** 8.14 ± 2.62** 7.70 ± 2.53** 7.31 ± 2.71** 6.88 ± 2.29** 6.50 ± 2.17** 6.23 ± 2.06** 6.07 ± 2.01** 5.93 ± 1.95** 5.92 ± 1.91**

2.84 2.84 2.83 2.83 2.82 2.82 2.83 2.82 2.83 2.83

6.51 ± 2.31 6.22 ± 2.35 6.17 ± 2.28 5.93 ± 2.31 5.64 ± 2.40 5.63 ± 2.36 5.66 ± 2.25 5.54 ± 2.27 5.51 ± 2.28 5.56 ± 2.16

0.06 0.001 0.01 0.09 0.11 0.06 0.037 0.07 0.07 0.06

14.47 ± 3.04** 13.69 ± 2.90** 12.59 ± 2.80** 11.38 ± 2.79** 10.17 ± 2.73** 9.26 ± 2.60** 8.80 ± 2.30** 8.73 ± 1.98** 9.11 ± 2.18** 9.57 ± 2.48**

0.93 0.93 0.84 0.63 0.52 0.50 0.44 0.42 0.49 0.52

SD12(1) SD12(2) SD12(3) SD12(4) SD12(5) SD12(6) SD12(7) SD12(8)

0.45 ± 0.21 0.70 ± 0.38 0.67 ± 0.30 0.64 ± 0.22 0.76 ± 0.27 0.84 ± 0.29 0.84 ± 0.29 0.84 ± 0.26

0.51 ± 0.23** 0.84 ± 0.37** 0.77 ± 0.25** 0.65 ± 0.16** 0.73 ± 0.21** 0.86 ± 0.24** 0.87 ± 0.24** 0.78 ± 0.21**

1.19 2.50 2.29 2.93 3.96 4.43 4.61 4.53

0.33 ± 0.14** 0.43 ± 0.16** 0.49 ± 0.17** 0.54 ± 0.19** 0.60 ± 0.21** 0.65 ± 0.22** 0.68 ± 0.23** 0.70 ± 0.22**

0.42 0.98 1.25 1.82 2.30 2.47 2.58 2.59

0.31 ± 0.18* 0.45 ± 0.13** 0.59 ± 0.16** 0.71 ± 0.21** 0.84 ± 0.24** 0.95 ± 0.27** 1.03 ± 0.31** 1.08 ± 0.37**

0.19 1.09 1.40 1.53 2.07 2.36 2.50 2.37

0.37 ± 0.11 0.51 ± 0.16* 0.52 ± 0.14* 0.62 ± 0.19* 0.77 ± 0.27 0.76 ± 0.24* 0.73 ± 0.18 0.79 ± 0.20

0.37 0.52 0.51 0.08 0.03 0.26 0.38 0.19

0.27 ± 0.09 0.44 ± 0.12* 0.65 ± 0.17** 0.87 ± 0.26** 1.13 ± 0.37** 1.34 ± 0.44** 1.42 ± 0.36** 1.41 ± 0.27**

0.39 0.46 0.23 0.91 0.87 0.60 0.96 1.14

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

MB mean ± std

1.07 1.23 1.43 1.51 1.44 1.31 1.25 1.22 1.21 1.19 204.52 ± 163.69** 284.74 ± 178.27** 338.95 ± 190.33** 361.92 ± 195.30** 362.82 ± 195.21** 354.50 ± 190.25** 349.62 ± 181.34** 348.87 ± 170.53** 353.53 ± 173.34** 353.36 ± 176.23**

29

Table 4 The percentage of significant differences (p < 0.05) between each pair of groups for Poincare plot features in 10 lags.

SNB MB pre-CCM CCM pre-KYM KYM

SNB

MB

pre-CCM

CCM

pre-KYM

100 100 100 15 97.5

90 95 42.5 97.5

60 75 97.5

72.5 90

95

KYM

Note – SD12: SD1/SD2; ES: effect size; ES  0:14: small effect; ES = 0:15–0.39: moderate effect; ES  0:40: large effect, emphasized by bold numbers.

44.86 ± 18.35 56.72 ± 24.55 59.41 ± 27.35 64.82 ± 30.50 69.33 ± 33.99 70.45 ± 35.27 71.11 ± 35.63 72.99 ± 37.68 73.83 ± 39.21 73.87 ± 39.68 2.32 2.44 2.40 2.38 2.36 2.34 2.33 2.31 2.30 2.30 61.39 ± 59.03 77.19 ± 77.06 81.51 ± 81.91 84.56 ± 85.67 89.13 ± 91.08 93.25 ± 97.19 95.08 ± 100.95 96.48 ± 103.48 97.80 ± 105.95 98.82 ± 108.73 Area(1) Area(2) Area(3) Area(4) Area(5) Area(6) Area(7) Area(8) Area(9) Area(10)

36.94 ± 21.92** 44.90 ± 26.00** 46.04 ± 26.35** 44.53 ± 25.16** 46.38 ± 27.24** 48.07 ± 28.04** 47.94 ± 27.61** 47.20 ± 27.52** 47.07 ± 27.84** 47.29 ± 27.81**

4.50 4.53 4.54 4.56 4.52 4.52 4.54 4.52 4.51 4.51

79.64 ± 58.95** 95.86 ± 61.79** 107.03 ± 71.94** 114.12 ± 78.81** 121.17 ± 85.36** 126.86 ± 91.62** 130.49 ± 95.50** 132.85 ± 99.30** 134.60 ± 101.46** 135.96 ± 103.07**

66.55 ± 35.30** 93.37 ± 48.10** 111.38 ± 60.54** 121.01 ± 68.26** 126.78 ± 71.78** 128.48 ± 72.17** 127.80 ± 71.06** 125.92 ± 69.47** 125.14 ± 68.84** 124.97 ± 69.16**

2.62 2.64 2.61 2.58 2.58 2.59 2.59 2.60 2.60 2.60

0.28 0.27 0.27 0.23 0.22 0.24 0.24 0.23 0.23 0.23

0.61 0.75 0.81 ± 0.22* 0.76 ± 0.17 2.60 2.61 0.87 ± 0.27 0.86 ± 0.24 SD12(9) SD12(10)

0.72 ± 0.13** 0.78 ± 0.19**

4.51 4.55

0.71 ± 0.22** 0.72 ± 0.22**

1.13 ± 0.41** 1.13 ± 0.42**

2.33 2.36

0.22 0.42

1.33 ± 0.31* 1.24 ± 0.39**

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

Table 5 provides the mean classification accuracies and sensitivities. The average classification accuracies of 99.14, 98.2, 96.76, 92.03, and 91.46% were achieved for pre-KYM, KYM, MB, CCM, and pre-CCM, respectively. The corresponding sensitivities are 99.87, 99.57, 99.04, 97.77, and 97.38%. 4. Discussions In this experiment, the nonlinear analysis was used to examine HRV dynamics in different states. To this end, HRV times-series of the Physionet database were applied. In an attempt to explain the HRV dynamics during different states, including CCM (pre and during the performance), KYM (pre and during the performance), and MB as compared to SNB, lagged Poincare measures, entropies, LZ, and LE were calculated. The results showed that not only HRV signals have different chaotic behaviors in SNB, MB, pre-CCM, during CCM, pre-KYM, and during KYM, but there is an obvious discrepancy between the HRV signals of meditators and non-meditators. Using SVM, pre-KYM, KYM, MB, CCM, and pre-CCM were discriminated from SNB with the average accuracy rates of 99.14, 98.2, 96.76, 92.03, and 91.46%, respectively. 4.1. Feature extraction The results showed that mean HRV is not significantly different in the conditions studied, except for pre-KYM (Fig. 2). May be this question arises that what is changing to stabilize the hemodynamic. Since HRV demonstrates characteristics of complex systems (Goldberger and Giles, 2006), traditional linear based methods such as the mean are not able to show HRV changes which are associated with the nonlinear interaction of physiological variables, such as autonomic arms (Castiglioni & Parati, 2011; Silva, Silva, Salgado, & Fazan, 2017). In this context, we studied different nonlinear based approaches. Investigation of the lagged Poincare indices in different lags (Fig. 3) showed that there are significant changes in the geometry of the HRV phase space. Therefore, it can be concluded that the chaotic behavior of the HRV signals varies in different states. This finding is in agreement with the previous conclusions on the meditative physiological signals, in which it was shown that increasing the lag leads

30

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

Fig. 5. Classification accuracies of nonlinear dynamical feature using SVM.

Table 5 Mean and standard deviation (mean ± std) of the classification accuracies and sensitivities for 50 times run of a 5-fold cross validation based SVM.

Accuracy Sensitivity

MB

Before CCM

CCM

Before KYM

KYM

96.76 ± 0.61 99.04 ± 0.53

91.46 ± 1.03 97.38 ± 0.94

92.03 ± 1.10 97.77 ± 0.89

99.14 ± 0.22 99.87 ± 0.22

98.20 ± 0.35 99.57 ± 0.39

to a change in Poincare plot’s shape during meditation (Goshvarpour et al., 2011; Goshvarpour & Goshvarpour, 2015). Poincare plots of different states have dissimilar patterns. The center of the plot is moved to the higher values for KYM (Fig. A.1), which represents the higher mean of the signal during KYM (in accordance with the results of Fig. 2). The most discrepancies in the lagged Poincare measures were observed for KYM and the least variations were observed for MB. The discrepancy of the points in lag 1 designates the correlation of the consecutive points in the time-series and is related to the amplitude of the signal. These results are consistent with the findings of the Peng et al. (1999), where it was shown that the amplitude of the HRV in KYM was significantly higher than the other states and the lower HRV amplitude was assigned to MB. CCM has the second most variations in the lagged Poincare measures. Therefore, it can be concluded that the Poincare indices are significantly different in meditators compared to nonmeditators (Fig. 3). Additionally, a lag variation causes the alteration in the shape of the Poincare plots (Fig. 3). According to the results, by increasing the lags, an increase in SD1 and a decrease in SD2 were perceived for all states (Table 3). Consequently, a cigar-shaped Poincare plot, with a higher correlation between HRV intervals, turned towards a round cloud of points in higher lags, which highlights a typical lack of correlation between HRV points. This result is consistent with the previously reported findings on HRV Poincare plots

(Goshvarpour et al., 2011; Goshvarpour & Goshvarpour, 2015). The results of Alvarez-Ramirez et al. (2017) showed that Hurst exponents reduced during meditation, which designates a reduction in the correlation between data points during meditation. SD1 illustrates the instantaneous beat-to-beat variability of a time-series and reveals the parasympathetic activity. In contrast, SD2 shows the overall beat-to-beat variability of the HRV signal. The area is under the vagal influence and the SD1/SD2 ratio of the Poincare plots is an index of sympathetic-parasympathetic balance (Goshvarpour et al., 2011; Kamen, Krum, & Tonkin, 1996; Tulppo, Makikallio, Takala, Seppanen, & Huikuri, 1996; Tulppo, Makikallio, Seppanen, Laukkanen, & Huikuri, 1998). Therefore, it can be concluded that KYM has the most effect on the function of sympathetic and parasympathetic nerves. In general, these results can be extended to the other extracted features. Examination of all features indicates dissimilar patterns between nonlinear HRV features regarding the influence of different states. For example, compared to the other states, significant decreases in the entropy values (ApEn, SaEn, and ShEn) were perceived for KYM. Therefore, an average irregularity of the signal decreases during KYM. Khandoker, Jelinek, Moritani, and Palaniswami (2010) reported that the lower entropy measures cause withdrawing the parasympathetic efferent pathway activities. Therefore, it can be concluded that KYM results in the withdrawal of the parasympathetic efferent pathway activities.

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

Although numerous nonlinear techniques have been employed in the past decades to evaluate HRV dynamics, interpretation of them in terms of parasympathetic and sympathetic activity has not been entirely described. Bolea, Pueyo, Laguna, and Bailo´n (2014) studied some linear and nonlinear HRV features during ANS blockades to evaluate the relation between parasympathetic/sympathetic activities and HRV indices. It was shown that parasympathetic blockade (by atropine administration) resulted in a significant decrease in nonlinear indices (SaEn and ApEn) (Bolea et al., 2014). In contrast, inhibiting sympathetic activity via propranolol caused a significant increase in ApEn. Their results showed that in the control condition, a natural increment in the sympathetic modulation (switching from the supine position to standing position) produces a significant decrement in the nonlinear indices (Bolea et al., 2014). In another study (Silva et al., 2017), selective blockade of cardiac autonomic control was proposed to examine the role of sympathetic and vagal function on nonlinear dynamics of HRV. They studied multiscale entropy (MSE) over short (1–5, MSE1–5) and long (6–30, MSE6 –

31

30) scales, where MSE at scale 1 is equivalent to SaEn. Their results showed that sympathetic and parasympathetic control differently affect entropy. The greatest effect of parasympathetic control of the heart to maintain the complexity of RR intervals is at low scales. Therefore, it can be concluded that our results on SaEn and ApEn (Fig. 4(e) and (f)) showed a higher parasympathetic activity for MB. While meditation (CCM and KYM) caused higher sympathetic activities. Higher LE in MB group indicates the more chaotic behavior of HRV signals (Table 3). In addition, LE indices were higher during CCM and KYM compared to preCCM and pre-KYM, respectively. However, these values were lower than that of the MB and SNB. The LZ increased during MB and decreased during CCM. These results confirm the previous findings on HRV dynamics during meditation (Goshvarpour & Goshvarpour, 2012c; Li et al., 2011; Song et al., 2013). Li et al. (2011) performed a study on CCM and KYM using the same database. They concluded that the dynamical complexity of HRV decreases in meditation. Song et al.

Fig. A.1. Poincare plots of different HRV states for lag 1 and lag 10. (a) SNB; (b) MB; (c) pre-CCM; (d) during CCM; (e) pre-KYM; (f) during KYM.

32

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

Fig. A.1. (continued)

(2013) showed that during CCM, the HRV becomes regular and the degree of multifractality decreases. By analyzing recurrence plots of HRV signals (Goshvarpour & Goshvarpour, 2012c), it was shown that nonlinear interactions of variables decrease in a meditation state. The results of the statistical test revealed that the total number of significant features (p < 0.05) was 189 out of 230 features (Table 3). The effect size analysis also demonstrated that 78.7% of all features (181 features) has a large effect (ES  0.4). These results put emphasis on the proficiency of the proposed feature extraction approach in discriminating between each state and SNB. However, some nonlinear features performed better in differentiating the states (Table 3). For instance, significant differences (p < 0.05) between all groups were observed for LZ, LogEn, ShEn, and ApEn. In addition, the results showed that there is a significant difference between each pair of classes for most nonlinear features. Totally, the results emphasized that there are significant differences in cardiac dynamics of SNB, MB, pre-CCM, during CCM, pre-KYM, and during KYM. Since there are significant differences between the nonlinear HRV

indices in pre-meditation techniques (CCM and KYM) and SNB, it can be concluded that performing meditation in a period of time can affect the base HRV dynamics. Therefore, the HRV patterns of meditators differ from the HRV patterns of non-meditators. 4.2. Classification In this section, we discuss the classification results. Discrimination of HRV dynamics was performed using a binary SVM classifier. The highest recognition rate was obtained for KYM. HRV features of pre-KYM and during KYM were classified with the average accuracy rates of 99.14 and 98.2%, and the average sensitivities of 99.87 and 99.57%, respectively. The results of MB classification were also promising. The average accuracy was 96.76% for MB. HRV dynamics were categorized with accuracy rates of 92.03 and 91.46% for CCM and pre-CCM, correspondingly. As far as the authors know, this is the first study that attempts to classify HRV signals of 2 meditation groups and MB from SNB. Previously, Phongsuphap and Pongsupap (2011) proposed

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

a system to classify HRV characteristics of Samadhi group and non-Samadhi. To this end, time - and frequency based features were extracted. Using Fisher discriminant analysis, the accuracy of 94.8% was obtained. Another study attempted to distinguish between HRV dynamics of CCM and pre-CCM (Goshvarpour & Goshvarpour, 2012a). To this end, LE and entropy were calculated and fed into the Quadratic, fisher, and kNN classifiers. The maximum accuracy of 92.31% was reported using Quadratic classifier. 4.3. Limitations and future works This experiment attempted to investigate the dynamics of the HRV signals of meditators and non-meditators using nonlinear analysis. There are significant differences between the parameters of different states. However, to confirm the findings, more data should be considered in future works. In addition, different meditation techniques should be analyzed in future. The KYM subjects were not age-matched with the other groups. In addition, some differences in the gender distribution could have further affected the results. These influences cannot be excluded. The effect of individual differences, such as gender- and age- differences, should be considered in future studies. Previous studies have shown that data length has significant impact on the results (Goshvarpour, Abbasi, & Goshvarpour, 2017a; Jiang et al., 2017) and some nonlinear measures like LE strictly demands a big amount of data. The length of some available signals in the Physionet dataset was limited, and we were not able to examine this parameter on our own results. Consequently, by collecting a new dataset, the effect of the data length should be considered in the future. In this study, a classification scheme was presented to discriminate between distinct ‘‘non-pathological” states. The objective of the classification step was to examine the separability of the HRV behaviors and characteristics in the different states. This goal can be used in medical and therapeutic applications. For example, as we know, meditation practice is used as a complementary medicine to reduce stress (Hoge et al., 2018). On the other hand, stress is associated with physiological changes such as a change in heart function (Lischke et al., 2018). Therefore, using the classification framework proposed in this study, in future work a procedure should be provided for managing stress using meditation techniques. It can help a meditator during a training session with proper feedback to bring the level of HRV to the appropriate level. 5. Conclusions In this study, a nonlinear-based feature extraction technique was proposed to classify HRV time-series. To this

33

end, lagged Poincare measures, entropies, LZ, and LE were calculated. In addition, statistical test and effect sizes were applied to examine significant differences between the states. Finally, the SVM classifier was employed for HRV classification. HRV recordings of Physionet dataset were used. The results showed that nonlinear measures had a significant difference between non-meditator and meditator groups. Precisely, there was a significant difference between nonlinear indices; this difference was not only between SNB and the other states, but also between different pairs of groups, and thus could provide further insight into underlying HRV dynamics. Additionally, the proposed framework was able to successfully classify the groups, with the mean accuracy rates of >96.76% for the pre-KYM, KYM, and MB groups. Compliance with ethical standards Funding: ‘‘This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors”. Ethical approval: ‘‘This article examined the HRV signals of Physionet dataset (Peng et al., 1999), which is freely available in the public domain. This article does not contain any studies with human participants performed by any of the authors.” Informed consent: Informed consent was obtained from all individual participants included in the study (Peng et al., 1999). Appendix A Poincare plot measures. Consider HRV time-series as X 1 ; X 2 ;    ; X N . Poincare plot is a 2D demonstration of the phase space of the HRV data. In this circumstance, 2 consecutive HRV data points are plotted ðX i ; X iþ1 Þ. To offer a quantifiable representation of data attractor, an ellipse shape is projected on the Poincare plot. In this study, the following indices were calculated: SD1 and SD2 were obtained from the minor and the major axis of the fitted ellipse on the Poincare plot, respectively. To calculate these indices, the basic statistical measures, i.e. standard deviation of RR interval (SDRR) and standard deviation of the successive difference of RR interval (SDSD), are used. These indices can be formulated by Eqs. (A.1) and (A.2). 1 SD21 ¼ SDSD2 ¼ cRR ð0Þ  cRR ð1Þ 2 1 SD22 ¼ 2SDRR2  SDSD2 2  2

¼ cRR ð0Þ þ cRR ð1Þ  2RR

ðA:1Þ

ðA:2Þ

where cRR(0) is the autocorrelation function for lag-0 and cRR(1) is the autocorrelation function for lag-1 RR inter

val. In addition, RR shows the mean of RR intervals.

34

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

Accordingly, the SD1/SD2 ratio and the area of the plot (p  SD1  SD2) were also computed. Totally, four indices, including SD1, SD2, SD1/SD2, and area were constructed from the HRV data in 6 states, including SNB, MB, preCCM, during CCM, pre-KYM, and during KYM in 10 different lags (lag1-10). Fig. A.1 provides a comparison between Poincare plots of the states for lag1 and lag10. In this figure, to plot these Poincare maps, 540 points of the data were selected for different scenarios, equivalent to three 180-sample windows. The colorization has been used to show the density of the points. The most repeated data were shown in red and the least frequent data were displayed in blue. As we know, the dynamics of chaotic systems are confined to the attractor, but the same states are never repeated. Therefore, Fig. A.1 confirms the chaotic behavior of the signals. Fig. A.1 shows that the KYM recordings have a larger variability. In this case, the plot center has shifted towards higher numbers. The Poincare plot center represents the mean of the signal. Fig. 2 graphically shows that the average HRV in KYM is greater than the other conditions. On the other hand, in the Poincare plot with lag 1, the dispersion of the points refers to the correlation of the successive data points and is related to the data amplitude. Previously (Peng et al., 1999), it was shown that the HRV amplitude in KYM is significantly higher than the other states. Appendix B. Lyapunov exponent Some features have been proposed to discriminate deterministic chaos from random or periodic behavior. A Chaotic system can be distinguished by its sensitive dependence on initial conditions. In these systems, small variations in the state variables will cause great differences in the behavior at some future points. In other words, in chaotic systems, two adjacent points in the trajectory at time 0 diverge widely from each other at time t. Consider two points at time 0 and at time t, which are adjacent at t = 0. Distances between them in the ith direction is shown by k dxi ð0Þ k and k dxi ðtÞ k, respectively. By quantifying the average separation rate (ki), LE is defined as Eq. (B.1). kdxi ðtÞk kdxi ð0Þk

¼ 2ki t ðt ! 1Þ

kdxi ðtÞk ki ¼ lim 1t log2 kdx i ð0Þk

ðB:1Þ

t!1

A positive LE shows sensitive dependence on initial conditions; therefore, it can demonstrate the deterministic chaos and the loss of predictability (Ruelle, 1979). Negative and zero LEs show a common fixed point and stable attractor, respectively. In this study, LE of HRV signals have been calculated for SNB, MB, pre-CCM, during CCM, pre-KYM, and during KYM conditions.

Appendix C. Entropy measures To investigate the irregularity and complexity of a signal, entropy is applied. To this end, different entropy measures have been introduced (Pincus, 1991). In this study, ShEn, LogEn, SaEn, and ApEn were calculated. Consider HRV time-series fX 1 ; X 2 ;    ; X N g with the probability mass function of P ðX Þ. ShEn and LogEn measures are formulated as Eqs. (C.1) and (C.2), respectively: n X P ðXiÞlog2 P ðXiÞ ðC:1Þ ShEn ¼  i¼1

LogEn ¼ 

n X

ðlog2 P ðXiÞÞ2

ðC:2Þ

i¼1

To calculate ApEn, two parameters, including m and r must be selected. Next, we need to compute the correlation integral C m ðrÞ. ApEnðm; r; LÞ ¼

Lm 1 X 1 logC imþ1 ðrÞ  L  m i¼1 Lmþ1



Lm X

logC mi ðrÞ

ðC:3Þ

i¼1

where r and m are the effective filter and the pattern length, respectively. Higher values of ApEn designate high irregularity, while lower values specify a more regular time-series. For SaEn, again two parameters of m and r must be chosen. In this study, m = 2 and r = 15% of the standard deviation of each time-series. By considering Xm data with the length of m, and the distance function of d ½X m ðiÞ; X m ðjÞ , we can count the number of vector pairs of length m and m + 1 which is symbolized as B and A. Finally, SaEn has the following equation: SaEn ¼ log

A B

ðC:4Þ

Appendix D. Complexity In this study, LZ (Lempel & Ziv, 1976) was applied to the HRV data. For calculating this measure, firstly, numerical data is transformed into a symbolic form. For example, binary data is composed of sequences of 0 and 1. A limitedlength sequence (S) is created as follows: S ¼ s1 s2 . . . sn ; 8si 2 A. The set of all subsequences of S is considered as the ‘vocabulary’ or ‘word’ (v(S)) of the S. Suppose S and W are two sequences. In addition, the concatenation of S and W is SW and when the last character of SW is deleted, the sequence is presented as SWp. Suppose we have a sequence of n-lengths. First, LZ is 1, S ¼ s1 ; W ¼ s2 ; and SW p ¼ s1 . To define a general statement, let S ¼ s1 s2 . . . sr and W ¼ srþ1 . If W 2 vðSW pÞ, then

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

srþ1 is a subsequence of S. Therefore, the new character does not add any new information to the S. Then, S remained unchanged and a new character (srþ2 ) is added to the W. If W 2 vðSW pÞ, LZ is increased by one. Simultaneously, S ¼ SW ¼ s1 s2 . . . sr srþ1 . . . srþi and a new character is added to the W (srþiþ1 ). The procedure is repeated until W holds the final character of S. Finally, LZ is defined as the number of different subsequences of S, which represents the amount of the new pattern resulting from the increase in the length of S. Generally, its normalized form is used in articles. Briefly, after transforming the numerical data into a symbolic data, the symbolic data sequence parsed to create distinct words. Encoding them with the length of L(n), the LZ is defined as (D.1): LZ ¼

LðnÞ n

ðD:1Þ

Appendix E. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.cogsys.2018.11.010. References Aftanas, L. I., & Golocheikine, S. A. (2001). Human anterior and frontal midline theta and lower alpha reflect emotionally positive state and internalized attention: High-resolution EEG investigation of meditation. Neuroscience Letters, 310, 57–60. Alvarez-Ramirez, J., Rodrı´guez, E., & Echeverrı´a, J. C. (2017). Fractal scaling behavior of heart rate variability in response to meditation techniques. Chaos, Solitons & Fractals, 99, 57–62. Bai, Y., Siu, K. L., Ashraf, S., Faes, L., Nollo, G., & Chon, K. H. (2008). Nonlinear coupling is absent in acute myocardial patients but not healthy subjects. American Journal of Physiology-Heart and Circulatory Physiology, 295, H578–H586. Beutler, E., Beltrami, F. G., Boutellier, U., & Spengler, C. M. (2016). Effect of regular yoga practice on respiratory regulation and exercise performance. PLoS One, 11, e0153159. Bolea, J., Pueyo, E., Laguna, P., & Bailo´n, R. (2014). Non-linear HRV indices under autonomic nervous system blockade. In Conf. proc. IEEE. eng. med. biol. soc., Chicago, IL (pp. 3252–3255). Braboszcz, C., Cahn, B. R., Levy, J., Fernandez, M., & Delorme, A. (2017). Increased gamma brainwave amplitude compared to control in three different meditation traditions. PLoS One, 12, e0170647. Brefczynski-Lewis, J. A., Lutz, A., Schaefer, H. S., Levinson, D. B., & Davidson, R. J. (2007). Neural correlates of attentional expertise in long-term meditation practitioners. PNAS, 104, 11483–11488. Castiglioni, P., & Parati, G. (2011). Present trends and future directions in the analysis of cardiovascular variability. Journal of Hypertension, 29, 1285–1288. Cerutti, S., Bianchi, A. M., & Mainardi, L. T. (1995). In Spectral analysis of the heart rate variability signal Heart Rate Variability (pp. 63–74). Armonk NY: Futura Publishing Company. Cohen, J. (1988). Statistical power analysis for the behavioral sciences. Hillsdale, NY: Lawrence Erlbaum Associates. Cortes, C., & Vapnik, V. (1995). Support-vector networks. Machine Learning, 20, 273–297. Gao, J., Fan, J., Wu, B. W., Zhang, Z., Chang, C., Hung, Y. S., ... Sik, H. H. (2016). Entrainment of chaotic activities in brain and heart during MBSR mindfulness training. Neuroscience Letters, 616, 218–223.

35

Goldberger, A. L. (2006). Giles f. Filley lecture. Complex systems. Proceedings of the American Thoracic Society, 3, 467–471. Goshvarpour, A., Abbasi, A., & Goshvarpour, A. (2017a). Do men and women have different ECG responses to sad pictures? Biomedical Signal Processing and Control, 38, 67–73. Goshvarpour, A., Abbasi, A., & Goshvarpour, A. (2017b). Fusion of heart rate variability and pulse rate variability for emotion recognition using lagged poincare plots. Australasian Physical and Engineering Science in Medicine, 40, 617–629. Goshvarpour, A., Goshvarpour, A., Rahati, S., & Saadatian, V. (2012). Bispectrum estimation of electroencephalogram signals during meditation. Iranian Journal of Psychiatry and Behavioral Sciences, 6, 48–54. Goshvarpour, A., Goshvarpour, A., & Rahati, S. (2011). Analysis of lagged Poincare plots in heart rate signals during meditation. Digital Signal Processing, 21, 208–214. Goshvarpour, A., & Goshvarpour, A. (2012a). Chaotic behavior of heart rate signals during Chi and Kundalini meditation. International Journal of Image, Graphics and Signal Processing, 2, 23–29. Goshvarpour, A., & Goshvarpour, A. (2012b). Classification of heart rate signals during meditation using Lyapunov exponents and entropy. International Journal of Intelligent Systems and Applications, 2, 35–41. Goshvarpour, A., & Goshvarpour, A. (2013). Comparison of higher order spectra in heart rate signals during two techniques of meditation: Chi and Kundalini meditation. Cognitive Neurodynamics, 7, 39–46. Goshvarpour, A., & Goshvarpour, A. (2015). Poincare indices for analyzing meditative heart rate signals. Biomedical Journal, 38, 229–234. Goshvarpour, A., & Goshvarpour, A. (2012c). Recurrence plots of heart rate signals during meditation. International Journal of Image, Graphics and Signal Processing, 2, 44–50. Goshvarpour, A., Shamsi, M., & Goshvarpour, A. (2013). Spectral and time based assessment of meditative heart rate signals. International Journal of Image, Graphics and Signal Processing, 5, 1–10. Hoge, E. A., Bui, E., Palitz, S. A., Schwarz, N. R., Owens, M. E., Johnston, J. M., ... Simon, N. M. (2018). The effect of mindfulness meditation training on biological acute stress responses in generalized anxiety disorder. Psychiatry Research, 262, 328–332. Hu, Y. H., & Hwang, J. N. (Eds.). (2002). Handbook of neural network signal processing. Electrical engineering and applied signal processing (Series). New York: CRC Press. Jadhav, N., Manthalkar, R., & Joshi, Y. (2017). Effect of meditation on emotional response: An EEG-based study. Biomedical Signal Processing and Control, 34, 101–113. Jevning, R., Wallace, R. K., & Beidebach, M. (1992). The physiology of meditation: A review. A wakeful hypometabolic integrated. Neuroscience & Biobehavioral Reviews, 16, 415–424. Jiang, M., Mieronkoski, R., Rahmani, A. M., Hagelberg, N., Salantera, S., & Liljeberg, P. (2017). Ultra-short-term analysis of heart rate variability for real-time acute pain monitoring with wearable electronics. IEEE int. conf. bioinformatics. biomed. (BIBM). 13–16 Nov. 2017, Kansas City, MO, USA. . Kamen, P. W., Krum, H., & Tonkin, A. M. (1996). Poincare plot of heart rate variability allows quantitative display of parasympathetic nervous activity. Clinical Science (London), 91, 201–208. Khandoker, A. H., Jelinek, H. F., Moritani, T., & Palaniswami, M. (2010). Association of cardiac autonomic neuropathy with alteration of sympatho-vagal balance through heart rate variability analysis. Medical Engineering & Physics, 32, 161–167. Lempel, A., & Ziv, J. (1976). On the complexity of finite sequences. IEEE Transactions on Information Theory, 22, 75–81. Li, J., Hu, J., Zhang, Y., & Zhang, X. (2011). Dynamical complexity changes during two forms of meditation. Physica A, 390, 2381–2387. Lischke, A., Jacksteit, R., Mau-Moeller, A., Pahnke, R., Hamm, A. O., & Weippert, M. (2018). Heart rate variability is associated with psychosocial stress in distinct social domains. Journal of Psychosomatic Research, 106, 56–61. Muralikrishnan, K., Balakrishnan, B., Balasubramanian, K., & Visnegarawla, F. (2012). Measurement of the effect of Isha Yoga on cardiac

36

A. Goshvarpour, A. Goshvarpour / Cognitive Systems Research 54 (2019) 21–36

autonomic nervous system using short–term heart rate variability. Journal of Ayurveda and Integrative Medicine, 3, 91–96. Peng, C.-K., Mietus, J. E., Liu, Y., Khalsa, G., Douglas, P. S., Benson, H., & Goldberger, A. L. (1999). Exaggerated heart rate oscillations during two meditation techniques. International Journal of Cardiology, 70, 101–107. Peressutti, C., Martı´n-Gonza´lez, J. M., & Garcı´a-Manso, J. M. (2012). Does mindfulness meditation shift the cardiac autonomic nervous system to a highly orderly operational state? International Journal of Cardiology, 154, 210–212. Phongsuphap, S., Pongsupap, Y., Chandanamattha, P., & Lursinsap, C. (2008). Changes in heart rate variability during concentration meditation. International Journal of Cardiology, 130, 481–484. Phongsuphap, S., & Pongsupap, Y. (2011). Analysis of heart rate variability during meditation by a pattern recognition method. Computers in Cardiology, 38, 197–200. Pincus, S. M. (1991). Approximate entropy as a measure of system complexity. Proceedings of the National academy of Sciences of the United States of America, 88, 2297. Reed, M. J., Robertson, C. E., & Addison, P. S. (2005). Heart rate variability measurements and the prediction of ventricular arrhythmias. QJM: An International Journal of Medicine., 98, 87–95. Ruelle, D. (1979). Sensitive dependence on initial condition and turbulent behavior of dynamical systems. Annals of the New York Academy of Sciences, 316, 408–416. Sarkar, A., & Barat, P. (2008). Effect of meditation on scaling behavior and complexity of human heart rate variability. Fractals, 16, 199–208.

Silva, L. E., Silva, C. A., Salgado, H. C., & Fazan, R. Jr., (2017). The role of sympathetic and vagal cardiac control on complexity of heart rate dynamics. American Journal of Physiology-Heart and Circulatory Physiology, 312, H469–H477. Song, R., Bian, C., & Ma, Q. D. Y. (2013). Multifractal analysis of heartbeat dynamics during meditation training. Physica A, 392, 1858–1862. Takahashi, T., Murata, T., Hamada, T., Omori, M., Kosaka, H., Kikuchi, M., ... Wada, Y. (2005). Changes in EEG and autonomic nervous activity during meditation and their association with personality traits. International Journal of Psychophysiology, 55, 199–207. Task Force of the European Society of Cardiology, North American Society of Pacing and Electrophysiology (1996). Heart rate variability: Standards of measurement, physiological interpretation, and clinical use. European Heart Journal, 17, 354–381. Tulppo, M. P., Makikallio, T. H., Seppanen, T., Laukkanen, R. T., & Huikuri, H. V. (1998). Vagal modulation of heart rate during exercise: Effects of age and physical fitness. American Journal of Physiology, 274, H424–H429. Tulppo, M. P., Makikallio, T. H., Takala, T. E., Seppanen, T., & Huikuri, H. V. (1996). Quantitative beat to beat analysis of heart rate dynamics during exercise. American Journal of Physiology, 271, H244–H252. van Ravenswaaij-Arts, C. M. A., Kollee, L. A. A., Hopman, J. C. W., Stoelinga, G. B. A., & van Geijn, H. P. (1993). Heart rate variability. Annals of Internal Medicine, 118, 4–47.