 Review
 Open access
 Published:
Nonlinear analysis of heart rate variability signals in meditative state: a review and perspective
BioMedical Engineering OnLine volume 22, Article number: 35 (2023)
Abstract
Introduction
In recent times, an upsurge in the investigation related to the effects of meditation in reconditioning various cardiovascular and psychological disorders is seen. In majority of these studies, heart rate variability (HRV) signal is used, probably for its ease of acquisition and low cost. Although understanding the dynamical complexity of HRV is not an easy task, the advances in nonlinear analysis has significantly helped in analyzing the impact of meditation of heart regulations. In this review, we intend to present the various nonlinear approaches, scientific findings and their limitations to develop deeper insights to carry out further research on this topic.
Results
Literature have shown that research focus on nonlinear domain is mainly concentrated on assessing predictability, fractality, and entropybased dynamical complexity of HRV signal. Although there were some conflicting results, most of the studies observed a reduced dynamical complexity, reduced fractal dimension, and decimated longrange correlation behavior during meditation. However, techniques, such as multiscale entropy (MSE) and multifractal analysis (MFA) of HRV can be more effective in analyzing nonstationary HRV signal, which were hardly used in the existing research works on meditation.
Conclusions
After going through the literature, it is realized that there is a requirement of a more rigorous research to get consistent and new findings about the changes in HRV dynamics due to the practice of meditation. The lack of adequate standard open access database is a concern in drawing statistically reliable results. Albeit, data augmentation technique is an alternative option to deal with this problem, data from adequate number of subjects can be more effective. Multiscale entropy analysis is scantily employed in studying the effect of meditation, which probably need more attention along with multifractal analysis.
Methods
Scientific databases, namely PubMed, Google Scholar, Web of Science, Scopus were searched to obtain the literature on “HRV analysis during meditation by nonlinear methods”. Following an exclusion criteria, 26 articles were selected to carry out this scientific analysis.
Introduction
In the last decade, a growing investigation on the role of yoga and meditation, especially on psychological and cardiac health using heart rate variability (HRV) features, has been observed. Relevant research works claim that yoga/meditation not only improves mental calmness and boosts concentration level, but also keeps the cardiovascular system healthy [1,2,3,4]. It is also to be noted that HRV signal has received tremendous response from the biomedical research community probably due to its ease of acquisition, noninvasiveness, and direct association with the functioning of the autonomic nervous system (ANS) [5,6,7,8,9,10]. The ANS plays a very important role in controlling various physiological processes, like the cardiovascular regulation, respiration, thermoregulation, urinary system, digestion, etc., and thus maintaining the homeostasis. Over the last few decades, different works [10,11,12,13] have demonstrated that meditation has interesting impact on the activity of the ANS, which functions via its two divisions, namely, the parasympathetic nervous system (PNS) and the sympathetic nervous system (SNS). Under different actions and health conditions, the predominance of these two subsystems automatically changes. Activation of SNS and PNS causes the release of chemical messengers; epinephrine and norepinephrine in case of the former, while acetylcholine in the latter [5, 10]. These hormones control the firing rate of the sinoatrial (SA) node, responsible for maintaining the HR. Parasympathetic activation decreases HR and sympathetic activation increases HR. This leads to a variation in interbeat (RR) intervals, which is represented by the HRV signal. Changes in HR dynamics under pathological, normal sinus rhythm, and different meditative states are shown in Fig. 1.
Literature study shows that linear (time or frequency domain) HRV analysis methods are more commonly used to distinguish the meditative state and the premeditative state. Out of them, time domain methods essentially perform statistical and geometrical calculations using the RR intervals time series. Although they are popularly used for studying the effect of meditation [15,16,17], they have the limitation of providing only the average information about the signal, while frequency domain methods are used for evaluating the spectral power distribution of the HRV signal in three distinct bands. These bands are classified as very low frequency (VLF), low frequency (LF), and high frequency (HF) bands. It is to be noted that another band, named the ultra low frequency (ULF) band (0\(\)0.0033 Hz) is seldom used due to lack of reliable findings of its correlation with the ANS. The significance of VLF (0.0033\(\)0.04 Hz) power is its association with thermoregulation, renin–angiotensin system, and peripheral vasomotor activity [18, 19]. LF power (0.04\(\)0.15 Hz) reflects the stimulation of both the SNS and PNS with a greater influence from the former. On the other hand, HF power (0.15\(\)0.4 Hz) reflects the stimulation of PNS. Based on these measures, an important parameter: the ratio of LF power to HF power is computed, which reflects the sympathovagal balance between the SNS and PNS. The values of these parameters are used to analyze whether there is any impact of meditation of ANS [11, 20,21,22]. In majority of these frequency domain studies, the fast Fourier transform (FFT) is used to estimate the power spectral density (PSD) for studying the effect of meditation. However, considering the nonstationary behavior of HRV, few [21, 23] of them have resorted to waveletbased HRV analysis during meditation.
Even though linear methods of HRV analysis are being used in the detection of mind–body interactions, they are not able to perceive the underlying complex patterns in the HRV signal. Meditation causes a change in HRV dynamics arising from respiratory sinus arrhythmia (RSA) as well as the changes in autonomic activity because of concentrated mind. For the study of various phenomena, like self similarity, deterministic chaos, irregularity in patterns, nonlinear tools are assumed to be more effective. Different nonlinear techniques and tools are already popular in various fields including hydrodynamics, mechanical, and electrical engineering, etc. Subsequently, in the study of human physiological conditions also, they are found to be quite effective. The complexity and the dynamical behavior of HRV signals during meditation are mainly analyzed with the help of different nonlinear parameters, including the largest Lyapunov exponent (LLE), the correlation dimension (CD), and the nonlinearity score (NLS) [12, 24, 25]. Besides, several entropy measures, like the Shannon entropy (ShEn), the approximate entropy, the sample entropy, the correlation entropy, the symbolic entropy, the Renyi entropy, the permutation entropy, etc., are also used [13, 26, 27]. The presence/absence of longrange dependence has also been studied in some studies to examine whether any significant changes occur during meditation based on fractal/multifractal analysis [15, 28, 29] and visibility graph [28, 30, 31] techniques.
A thorough literature search on “HRV analysis for studying the effects of meditation” returns papers which are primarily focussing on linear analysis techniques and a massive vacuum is observed in the reviews related to nonlinear HRV analysis for studying the impacts of meditation. As numerous original works have demonstrated the nonlinear property of HRV signal, it is of paramount importance to have an analysis on the published papers related to the aforementioned subject. This led us to cumulate the findings of various nonlinear methods of HRV signal to analyze the impacts of meditation, followed by their critical analysis.
In order to provide a critical analysis on the works, we first mention different nonlinear parameters for characterization of HRV dynamics during and before the practice of meditation and yoga. Relevant theoretical backgrounds are briefed to extract specific knowledge and understanding of their usefulness in studying the dynamics of HRV signal. Next, we add our perspectives to aid in the quest for continuous and breakthrough findings on this multidisciplinary research domain.
The rest of the paper is organized as follows. "Results" section provides a details of the nonlinear HRV methods with their theoretical cores for studying the nonlinear behavior of HRV and the findings related to nonlinear HRV analysis under meditative intervention. "Discussion" section provides an overall summary of different approaches. The conclusions and a few future research directions are discussed in "Conclusion and future research trends" section. Finally, "Methods" section provides a detailed search strategy for the selection of relevant research articles analyzed in this study.
Results
Distribution of research articles on nonlinear HRV analysis in the last two decades
The literature search returns around 389odd articles on nonlinear HRV analysis for studying the cardiac conditions during meditation/yoga. Figure 2 shows yearwise numbers of publications on HRV studies on the above context, which clearly indicates that it is one of the highly researched areas by the scientific community and the trend is fairly on a rising scale. It is worth mentioning that in this study, we have considered Chi, Kundalini, mindfulness, Zen, slow/deep breathing, paced breathing, focused meditation as meditative interventions although there are copious forms of meditation. This is because our aim is to study various nonlinear methods useful for studying meditative interventions rather than to evaluate different meditative interventions. A pie diagram representation of the distribution of meditation techniques used in research articles is shown in Fig. 3.
Nonlinear methods for studying HRV dynamics
A brief account of the various nonlinear methods used for the study of HRV under meditative interventions; the databases, nonlinear parameters, and significant observations construed from them are detailed in Table 1.
Next, we elaborate the theoretical background and the mathematical formulation of some of the widely used techniques belonging to nonlinear regime; highlighting the parameters used for studying the dynamics of HRV as follows.
Application of Poincaré plots
Poincaré plot is popularly used to study the behavior of ANS by investigating the shortterm and longterm HRV [46, 47]. In this technique, each RR interval is plotted against its previous RR interval (i.e., \(RR_{n+1}\) vs. \(RR_{n}\) for a delay of 1), which forms a scatter plot as shown in Fig. 4. The shape of the plot gives a qualitative assessment, while the standard deviation parameters, SD1 and SD2 (dispersions along the minor and major axes of the fitted ellipse, respectively) could provide a quantitative assessment of the underlying dynamics of HRV. Particularly, SD1 corresponds to the standard deviation of instantaneous (shortterm) variability of RR intervals, which indicates the activity of PNS, while SD2 corresponds to the standard deviation of continuous longterm variability, which is an indicator of the activity of SNS [37, 47]. The SD1/SD2 ratio is the relative measure of shortterm and longterm variability of RR intervals.
A new approach of Poincaré analysis is also reported in the form of heart rate asymmetry (HRA) measure by Karmakar et al. [48], and later applied in [42] to assess its sensitivity towards delay. Mathematically, the HRA index, popularly known as Guzik’s index (GI), is given by [48]:
where \(C(P_{i}^{+})\) is the number of points (cyancolored dots) above the line of identity (LOI) as shown in Fig. 4, N is the number of RR intervals, \(D_{i}^{+}\) is the distance of ith point above the LOI, and \(D_i\) is the distance of ith point from the LOI given by:
It is reported that HRA index of meditators (Chi and KYM) are significantly different from nonmeditators especially at higher delays [42, 43]. Similarly, Rohila and Sharma [43] observe that there occurs a significant increase in GI index too during meditation, while Goshvarpour et al. [27, 35] have observed that SD1/SD2 ratio also increases during meditation, which connotes the dominant parasympathetic activity, while Dey et al. [49] have presented a 3D frequency delay plot, which shows that the effect of Chi and KYM on the ANS are not alike, rather their impacts are completely opposite. Goswami et al. [33] use a modified version of Poincaré plot, termed as the secondorder difference plot and observe that during meditation (Chi and KYM), the axis of the cluster rotates anticlockwise. Furthermore, correlation coefficient obtained from consecutive firstorder differences of RR intervals also increases during meditation. However, their results are not further validated by statistical significance tests. In order to compare the findings of the standard HRV features along with some of the commonly used nonlinear features (DFA scale exponent and PEn), against the findings of the existing works, we have demonstrated the simulated results in Table 2. It can be observed that only a few of the standard HRV features are able to distinguish the two states of mind–body states significantly. This portrays the importance of nonlinear analysis in distinguishing the HRV data of compared groups having subtle differences.
Nonlinear parameters in phase space
Application of phase space representation In case of nonlinear analysis of HRV, it is very common to transform the time series into the phase space form. The rationale behind the transformation is to capture all possible dynamical states of the 1D signal, which would not otherwise be possible from the magnitudes of a single variable. For that purpose, the signal is embedded into higher dimensions. The minimum embedding dimension (MED) provides the minimum number of variables required to represent the complete dynamics. But, in that process it is essential to ensure that the reconstructed attractor of different dimensions preserve the salient properties of the original attractor [24, 50]. The choice of proper embedding dimension and time delay play an important role in the efficient reconstruction of the attractor. Mathematically, the RR interval (HRV) time series of N data points, denoted by \(\text{RR}(i)=[\text{RR}(1), \text{RR}(2),.....,\text{RR}(N)]\) can be transformed into an m dimensional phase space with time delay, \(\tau \) by Taken’s embedding theorem [51] as follows:
where each time delay vector, Y(i) can be expressed as \(Y(i)=[\text{RR}(i), \text{RR}(i+\tau ), \text{RR}(i+2\tau ),...,\text{RR}(i+(m1)\tau )]\) and there are \(N_v=N(m1)\tau \) such time delay vectors.
For the estimation of MED, generally the methods based on the singular value decomposition (SVD) [53], the variation of the exponent of powerlaw based correlation integral [54], the count of false nearest neighbors (FNNs) based on Kennel’s [52] and Cao’s [50] approaches are commonly used. However, the method [50] has received wide acceptance among the chaos analysts because of the computational ease and reliability. Another advantage is that the estimation of FNNs is threshold independent. The attractor gets unfolded if the signal is embedded at higher dimension than MED, which can facilitate the study of stability, predictability, correlation behavior, and regularity. However, it is to be noted that the optimal timedelay value should be known prior to the estimation of MED. Average mutual information (AMI) and autocorrelation techniques are mostly used to determine the right delay at which these functions return the least value. Figures 5 and 6 show estimation of MED based on Cao’s and Kennel’s methods, respectively, for an HRV signal of meditators (PhysioNet). Here, E2(m) is given by \(\frac{E^*(m+1)}{E^*(m)}\), with \(E^*(m+1)\) and \(E^*(m)\) indicating the averages of the distances between the neighboring vectors at dimensions \(m+1\) and m, respectively, as given below:
where \(RR_{n(i,m)+m\tau }\) is the nearest neighbor of \(RR_{i+m\tau }\) in the mdimensional phase space. As reported in [50], E2(m) is specially tailored to distinguish deterministic signals from stochastic signals. For the test HRV signal, it attains the peak value at \(m=2\) and then saturates beyond \(m=3\). Saturation of E2(m) after \(m=3\), gives the value of MED, i.e., \(m+1=4\). As per Kennel’s method since the lowest percentage of FNN is found at \(m=4\) with threshold of 15% false neighbors, the MED = 4.
It can be observed from Fig. 7 that mutual information drops sharply and it remains negligible from \(\tau =4\) onward. With lower \(\tau \) values, the MED required to capture the attractor of a signal is higher than that with optimal \(\tau \) values. However, timedelay information is more relevant for continuous time signals, while for the analysis of discrete time signal (e.g., HRV), the timedelay of 1 is sufficient [50].
Phasespace transformation is found to be the basis for various types of nonlinear analyses, since it provides the scope to analyze a 1D time series under higher dimensions. This led to its use in studying the effects of meditation from the aspects of asymptotic stability [24, 27, 55], time correlation of phase trajectories [10, 27, 55], dynamical complexity [10, 14, 26, 37, 41, 45], etc. Nevertheless, the computation of MED from phase space representations using FNN and Cao’s techniques (shown above) have revealed that HRV signals under meditative and premeditative conditions can be best studied within MED=3. This indicates that even during meditative states, the HRV time series maintains the property of lower order deterministic chaos. Thus a general belief that the meditation may develop more regular and predictive cardiac dynamics is actually false, rather the cardiac dynamics still remain complex. It is worthwhile to mention that a complex heart dynamics indicates better adaptability of heart, whereas the monotonous heart dynamics enunciates pathological heart conditions.
Application of correlation dimension (CD) It gives the assessment of the complexity of a time series [54]. Since, for many complex signals, CD value is a fractional quantity, it is also known by fractal dimension (FD). Higher the magnitude, greater is the dynamical complexity. CDbased dimension calculation is having more adoption in biomedical applications as compared to the other techniques including, boxcounting dimension, Kolmogorov capacity dimension, and information dimension, since it needs relatively lesser data for its calculation. It can be obtained from a time series by finding correlations sum calculated from its data points as reported in [56]:
where
is the correlation integral of all the correlation functions, \(N_v\) is the number of states, and \({\mathcal {H}}\) is the Heaviside step function, which determines whether the distance between m dimensional vectors, \(Y_{m}(i)\) and \(Y_{m}(j)\) is within a sphere of radius, r or not. \(C_{m}(r)\) is found to have a powerlaw relationship with r, whose exponent gives the measure of CD [54]:
Raghavendra and Dutt [24] observed reduced MED and CD values during meditation (Chi and KYM) as shown in Fig. 8, which they considered to be due to a state of overwhelming calmness. Although CD is considered as an indicator of vagal activity, it also characterizes the sympathovagal balance [57, 58]. This demonstrates the depressed vagal activity during meditation. However, it is to be mentioned that estimation of CD requires large number of data, which makes its use dubious for the given dataset having smaller data length.
Application of Lyapunov exponent (LE) LEs indicate whether the underlying system of a time series is chaotic or not [59]. It determines the sensitivity of a system to the initial conditions by calculating the average rate of separation of infinitesimally close trajectories. A positive value of the largest LE (LLE) of a series means that the underlying system is chaotic, while negative LE describes the asymptotic stability (convergence of close trajectories). For a phase space representation of m dimensional dynamical series having N points, let the spectrum of LEs be \((\lambda _{1},\lambda _{2},...\lambda _{m})\). If the rate of divergence between any two trajectories is \(\delta Z(t)\) with initial separation vector \(\delta Z_{0}\), the LLE at a time t can be obtained by [24]:
Fig. 9 shows that LLE is higher during meditation [24, 27, 55], which is considered to be associated with enhanced alertness during meditation. The higher value of LLE during meditation may be considered as the manifestation of stronger neurocardiovascular couplings as well as the couplings between the ANS and other sensory systems [60, 61].
Application of recurrence plot and recurrence quantification analysis The recurrence plot (RP) is an array of dots, graphically represented by a squaresized (\(N \times N\)) image. A dot is placed at (i, j) whenever the delay vector \(Y_m( j)\) is sufficiently close to another vector \(Y_m(i)\). Therefore, the index (i, j) corresponds to the time information, and RP naturally describes the time correlation between two vectors. Its not restrained by stationarity, dimension size, and data length. Mathematically, the vectors, \(Y_m(i)\) and \(Y_m(j)\) are recurrent if they are within a cutoff distance (\(\varepsilon \)) [10, 62, 63], i.e.:
For clarity, recurrence plots for HRV signal corresponding to meditative and premeditative states are presented in Figs. 10 and 11 under a setting of embedding dimension, m=3 and time delay=1. It can be seen that during meditation, more prominent short lines parallel to the diagonal of the RP are found. This indicates that some states recur during meditation hinting a more periodic kind of oscillations. These periodic patterns may be due to the controlled breathing cycles during meditation. This illustrates the importance of respiratory data while analyzing the effects of meditation to arrive at a more profound decision. On the other hand, it is equally important to understand the subtle patterns in the RP, which may not be ascertained by visual inspection. Nonetheless, the correlations between different subsystems can be better analyzed by mathematical tools, such as recurrence quantification analysis (RQA) using the recurrence rate (REC), the determinism (DET), the maximal length of diagonal structures (\({L}_{\text{max}}\)), the Shannon entropy, and the trend based on the RP [63]. Later on, Marwan et al. [64] have derived few more quantifiers namely, laminarity, trapping time, and maximal length of vertical structures to extract further information from RP.
Goshvarpour and Goshvarpour [27, 55] observe that big black rectangular structures in the RP change to small rectangular patches when subjects sink in the meditative state from premeditative state. Sarkar and Barat [10] have observed that RQA parameters remain nearly constant at varying time scales during Chi meditation, whereas the parameters vary abruptly before meditation. On a nutshell, the information about the temporal correlations over certain time intervals and the information about the impact of RSA on baroreceptor sensitivity as well as the autonomic activity can be strongly analyzed by RP tools.
Application of entropy in phase space Entropybased parameters are also computed in phase space to reveal the regularity, complexity, and predictability of a system [26, 65, 66]. A brief overview of a few of them are given below:
Approximate entropy (ApEn) It provides the degree of regularity and unpredictability of the variations in a time series [67], which is given by:
where
Here, \(D_{i}^m(r)\) is almost the same as \(C_{i}^m(r)\) of Eq. 6, differing only for its consideration of distance measurement of time delay vectors with itself; \(\phi ^{m}(r)\) is the average of the correlation functions of all the time delay vectors with dimension m, and r is the tolerance limit.
Sample entropy (SampEn) Since approximate entropy is dependent on the data length and includes self matching of time delay vectors, it can give rise to biased result [67]. So it has been modified to derive the sample entropy (SampEn), which provides better realization of complexity of a time series. Its mathematical expression is given by [68]
where A and B are the number of time delay vectors whose Chebyshev distances are less than r in \(m+1\) and m dimensional phase space, respectively.
Permutation entropy (PEn) This entropy is suitable for identifying whether a temporal dynamics is originated from stochastic process or deterministic chaos, and follows the identical mathematical formulation as used in ShEn [66]. After embedding the time series into phase space, the unique ordinal patterns of the embedding vectors are designated as \(\pi \) forms. The maximum possible patterns configurable by the vectors of m dimension is given by m!. It is mathematical expression is given by:
Basescale entropy (BSEn) BSEn is another complexity metric suitable for short time series. For its computation, firstly the base scale (BS) is obtained from the phase space representation of the series and then symbolic sequences are formed on the basis of BS. Its mathematical calculation is given by [69]:
The symbolic sequences are coded in terms of symbols (say, a = 0, 1, 2, 3), which indicates about the amplitude resolution. In this case, for m dimensional phase space, there can be maximum of \(4^m\) (i.e., \( \{{\text{length(a)}}\}^m\)) different forms of \(\pi \). For each \(\pi \), the probability of occurrence is calculated to obtain BSE as given below:
Increment entropy (IncrEn) IncrEn evaluates the dynamical complexity of a time series by encoding each increment of it onto a word of two letters, where the first letter denotes the sign of change and the second letter denotes the amount of change [70]. Let these codewords be denoted by \(w_{n,k}\) for n=1, 2,...,\((Nm)\) and k = 1,..,m for an embedded increment series, HRV of length N; firstly, the successive differences are obtained and then the resultant increment series is embedded into m dimensional vectors, \(y_n\) for n = 1, 2,...,\((Nm)\). Each element of these vectors (\(y_n\)) are mapped onto a word \(w_{n,k}\) for k = 1,..,m having two letters, say \(s_{n,k}\) and \(q_{n,k}\) corresponding to sign and symbolic increment, respectively. These m words (\(w_{n,k}\)) in a vector, form a full word denoted by \(w_n\) [41]:
If \(t(w_{n})\) denotes the number of occurrences of nth unique word, the probability of occurrence of each unique word is given by:
The increment entropy, H(m) for a resolution level (r) is given by:
Deka and Deka [41] have observed reduced IncrEn value especially during KYM and marginally lower entropy during Chi meditation than that before meditation, though the results are not statistically significant. This indicates that there may be a reduced irregularity and dynamical complexity during meditation. However, many other markers have demonstrated more significant decrease in entropy values during meditation. For example, BSEn [26], SampEn and ApEn [27], PEn [38] return significantly lower value (\(p<0.05\)) for the HRV during meditation. However, based on multiscale entropy (MSE) analysis as introduced by Costa et al. [65], Sarkar and Barat [10] show that for higher scales, the SampEn values are quite higher during meditation. In [14], Deka and Deka have used multiscale PEn (MPE) to examine the dynamical complexity of temporal structures in HRV signals of meditators and observed consistent results with that of [10]. Consistent results are also obtained from a recently developed improved multiscale distribution entropy measure [45]. It provides an interesting fact that actually complexity of HRV signal increases during meditation.
Apart from these entropy markers based on Shannon entropy theory, generalized forms of entropies, such as Renyi entropy, Tsallis entropy, Kolmogorov–Sinai (KS) entropy as well as correlation measurebased correlation entropy are also used in practice.
Renyi entropy (ReEn) ReEn is the generalized form of ShEn, which is suitable for quantifying the complexity of multifractal time series. The mathematical formulation of ReEn of order q can be given as:
where n is the number of possible outcomes of a signal/series and \(q\ge 0\), but \(q\ne 1\).
Tsallis entropy (TsallisEn) TsallisEn is another generalized form of ShEn, which was introduced with the aim to provide a better analysis of especially, nonextensive systems, very common in natural phenomena. For an order q, TsallisEn can be given by [71]:
Kolmogorov–Sinai (KS) entropy KS entropy also quantifies the dynamic behavior of a time series, by determining the changes in entropy by partitioning the data at each iteration. To obtain KS entropy, the HRV time series is first mapped using measurepreserving transformation, \(T:\text{RR}\rightarrow {\text{RR}}\). The transformed RR series is then partitioned into several parts and each partition elements are considered as intervals \(I_i\). Then with the elements in partition at time n falling into an interval \(I_{i_1}\), while elements falling into interval \(I_{i_2}\) at time \(n+1\) and so on, block entropies of block size m can be obtained using the joint probability \(P_{i_1,i_2,..,i_m}\) function by [59]:
where \(h_q=\underset{P}{\sup }\ \underset{m\rightarrow \infty }{\lim }\frac{1}{m}[H_q(m,P_\epsilon )]\) and \(\underset{P}{\sup }\) means finding the supremum over all possible P, provided \(\epsilon \rightarrow 0\).
Correlation entropy (CorrEn) It gives the similarity measure between two random variables in the neighborhood of joint space. For this, the data are transformed into a higher dimensional Hilbert space and afterwards the correlation of the data in Hilbert space is obtained. If X, Y are two series of same dimension then CorrEn can be obtained by [37]:
where H(X) and H(Y) map the data X and Y to higher dimensional Hilbert domain. \(K_{\sigma }\) is a Kernel function, which could be a polynomial function, radial basis function (RBF) or sigmoid function, etc.
Goshvarpour and Goshvarpour [37] evaluate information theoretic descriptors, namely, the CorrEn and the Cauchy–Schwarz divergence (CSD) and observe very low CorrEn values during Kundalini yoga for smaller kernel sizes. It indicates higher relaxed state and lower SNS tone during Kundalini yoga meditation. However, as the kernel size increases, these values further decrease marginally. They further conclude that CSD fails to demonstrate the nonlinear similarities. Again, diffusion entropy analysis (DEA) as shown in Fig. 12 demonstrates regular repetitive oscillations during meditation, indicating the loss of longrange correlation. In [72], Porto et al. have found that although entropy markers, such as TsallisEn, ReEn could not provide significant difference between slow breathing and normal breathing conditions. They observed that dynamical complexity in case of slow breathing reduces and this leads to an reduced vagal control.
Fractal/multifractal analysis
Applications of detrended fluctuation technique Like many natural objects and signals, such as mountains, leaves of plants, DNA sequences, time series of network traffic, HRV time series also possesses the fractal behavior. Fractal object/process should have a fractional dimension and possess selfsimilarity behavior under different spatial/temporal scales. Fractal dimension (FD) is the most fundamental measure to check the fractality and hence the complexity of a signal, which is based on the degree of spacefilling of a curve constituting the signal. Even though CD, boxcounting dimension, Hausdorff dimension can provide the measure of FD, Hurst exponent and detrended fluctuation analysis (DFA) based FD measures are widely used [32, 73].
DFA method reveals the possession of selfsimilarity (fractal) behavior by eliminating any kind of trends. This trend unless removed, could mislead the interpretation of fractal analysis as it could add unnecessary signal components extraneous to the underlying dynamics of interest. The mathematical formulation of this method is given below [73]:
where \(\text{RR}_{cum}(m)=\sum _{i=1}^{m}[RR(i)\overline{RR(i)}]\), \(\overline{\text{RR}(i)}\) is the mean of RR interval series, \(\text{RR}_{cum}(m)\) is the cumulative sum, \(F_{rms}(n)\) is the rms fluctuation of the series at scale n, \(\text{RR}_{t_{n}}(m)\) is the trend (leastsquare fit line) of RR(m) at scale n, and \(\alpha \) is the scaling exponent which defines the selfsimilarity behavior of a signal. For \(0\le \alpha \le 1\), Hurst exponent, H = \(\alpha \) and for \(1<\alpha \le 2\), H = \(\alpha 1\). Here, H value reflects the correlation behavior of a series with 0 < H < 0.5 indicating antisymmetric correlation and 0.5 < H < 1 indicating longrange correlation. From the H value, FD can be obtained by: FD = 2H [74, 75]. It is to be noted that for the computation of asymptotic \(\alpha \), a very long time series is required (approximately 9000 data points), which is often unavailable as the physiological signals are concerned. This led to its evaluation for short time scales (4–16 data points), by shortterm DFA exponent (\(\alpha _s\)), which is also found to be effective in studying selfsimilarity property [27, 73]. In a work, Sarkar and Barat [10] observe that the longrange correlation of the HRV data is destroyed during Chi meditation as shown in Fig. 13, where there is a scaling crossover indicating the effect of meditation. In [12], Raghavendra and Dutt compute the mean and standard deviation of FD values for various scales (from 1 to 20) and it is found that FD values are significantly lesser during meditation ( \(p<\)0.05). During deep meditation, Guo et al. [39] have found that both shortterm and longterm scaling exponents have significantly increased. In another work, Kamath [15] observes decrement in HFD, which is interpreted as the increase of LFp owing to the shift of the RSA by parasympathetic mediation. AlvarezRamirez et al. [36] have observed a decrement in Hurst exponent (HE) during chi and KYM meditation from the rescaled range (R/S) analysis. From the time series, firstly overall rescaled range (R/S) is obtained. The constant of proportionality between the logarithm of R/S value and the logarithm of time scales gives the exponent known as HE. Reduced HE indicates uncorrelated HRV dynamics and destruction of longrange correlation due to meditation. In another study, Gronwald et al. [76] have observed that a reduced shortterm DFA exponent manifests the reduced dynamical complexity and withdrawal of organismic system to maintain homeostasis under the influence of central nervous system.
Moving forward, it was observed with time that HRV time series do not have the same fractal behavior all throughout its length, rather they possess multifractal behavior. In multifractal analysis method, the fractal behavior is analyzed under different scales locally. More mathematical details about MFDFA can be found in [77]. It is observed by Song et al. [78] that during Chi meditation, multifractal spectrum width reduces as compared to that in premeditative state. This indicates a reduced multifractality during meditation.
Application of visibility graph technique Visibility graph (VG) essentially converts a fractal time series into a scalefree graph. In this method, each data points are indicated by nodes and the nodes are connected by edges as shown in Fig. 14. If \(X_{p}\) and \(X_{q}\) are two data points at positions p and q, respectively, in a time series then they can be connected by an edge if it satisfies
where p, q and \(s\in Z^+\) and \(p<p+s<q\). The degree of a node is the number of edges connected with a node, and the degree distribution, P(k) is the fraction of the nodes with degree k. It is to be noted that as per the law of visibility graph, degree distribution follows powerlaw with the degree of a node, i.e., \(P(k) =k^{\lambda _{p}}\). The constant of proportionality, \(\lambda _{p}\) is called power of scalefreeness in visibility graph (PSVG) [28].
Based on VG, several complexity markers were derived, namely, medium articulation (\(\text {MA}_{VG}\)), graph index complexity (GIC), Visibility entropy (VGEn), offdiagonal complexity (ODC), etc. [79]. It was found that graph with medium number of edges is the most complex, whereas both sparse type and fully connected graphs are in fact less complex. \(\text {MA}_{VG}\) value attains the maximum value for such graphs with medium edges, which follows \(n_{VG}^{1.5}\), with \(n_{VG}\) being the number of nodes. Again, ODC provides the diversity in the value of node–node edges. More mathematical details about \(\text {MA}_{VG}\) and ODC can be found in [79]. However, since GIC and VGEn have been more frequently used in time series analysis including HRV, we provide below a brief mathematical details of them:
where D is the total number of unique degrees. On the other hand, GIC is determined from the eigenvalues of adjacency matrix of a graph. The largest eigenvalue (\(\lambda _{max}\)) for a graph lies with the range [\(2\cos (\frac{\pi }{n_{VG}+1}), (n_{VG}1)\)]. Mathematically, GIC is given by [40, 79]:
where \(c=\frac{\lambda _{max}2\cos (\frac{\pi }{n_{VG}+1})}{n_{VG}12\cos (\frac{\pi }{n_{VG}+1})}.\)
The PSVG depicts the amount of complexity and selfsimilarity of a dynamic time series, which is also linearly related to the HE. Higher the PSVG, more is the complexity. Bhaduri et al. [28] observe increase in both multifractal spectrum width and PSVG during Chi meditation. Unlike MFDFA technique, PSVG can be applied for small data size also. However their results contradict to that of [12, 29, 78]. On the other hand, Jiang et al. [30] observe that P(k) decreases as the degree of node k increases before meditation as shown in Fig. 15. Although during meditation it is found to be small for \(k<\)8, increases consistently to reach the peak value at \(k=\)11 indicating distinct change in dynamics. However, at large values of k, the P(k) distribution follows a powerlaw tail, which indicates that the longrange correlation of HRV time series is similar for both premeditative and meditative states. Again, from the complexity measure GIC, Nasrolahzadeh et al. [40] have found that during meditation, GIC values are much higher than that before meditation. They corroborate that the higher complexity associated with meditative state, reflects better adaptability of heart under different conditions (i.e., healthy ANS activity).
Application of higher order spectral (HOS) analysis HOS analysis provides the spectral analysis of higher order (\(>2\)) moments or cumulants of a signal. It is very effective in extracting hidden information from nonlinear and nonGaussian signal (e.g., HRV), which cannot be obtained from regular power spectrum. This is because regular power spectrum provides no information about the Gaussianity/nonGaussianity behavior of a signal. Besides, power spectrum (secondorder moment) is inadequate to provide the phase relations among spectral components [80,81,82].
If X(k) is an HRV time series, the nthorder spectrum \(S_n(f_1,f_2,...,f_{n1})\) can be obtained by taking the Fourier transform of the nthorder cumulant, \(c_n(t_1,t_2,...,t_{n1})\) as given below:
The spectrum obtained from 3rdorder cumulant, known by bispectrum is very popular among researchers to capture new information from HRV signal [34, 81,82,83]. The mathematical expression of bispectrum is given by:
where \(E[\cdot ]\) is the expectation operator, X(f) represents the Fourier transform of a signal X(k), and \(*\) indicates the complex conjugate operation. It is evident from Eq. (28) that bispectrum is capable of providing information about the phase coupling or spectral interactions of a signal at frequencies \(f_1\), \(f_2\), and \(f_1+f_2\). Similarly, trispectrum of a signal can provide details of the spectral interactions of \(f_1\), \(f_2\), \(f_3\), and \(f_1+f_2+f_3\) and given by:
In order to eliminate the impact of the variability of power spectrum across different signals on the salient phase coupling information, often bispectrum and trispectrum are normalized. The normalized versions of these two polyspectra are known by bicoherence and tricoherence, respectively. Their mathematical expressions are given by:
From the bispectrum analysis of HRV signals for meditative and premeditative states, Goshvarpour and Goshvarpour [34] have observed distinct differences in bispectrum amplitudes for the two states. They have observed an increase in mean bispectrum amplitude during meditation in case of KYM practitioners, while a significant decrease (\(p<0.5\)) in bispectrum amplitude during Chi meditation. However, their results are found to be containing bispectral components up to the range of \(\sim \) 15 Hz, which is found to be exceeding the spectral bandwidth of HRV signal. This may also happen due to the consideration of higher sampling frequency for the signals in the study. From, bispectral analysis of the HRV signals as shown in Figs. 16 and 17, it can be observed that oscillations are more periodic during meditation. However, the power concentrations in the spectral bands are more evident from bicoherence plots as portrayed in Figs. 18 and 19. It can be seen that before meditation, prominent spectral components are present in both 0.05\(\)0.15 Hz and 0.15\(\)0.45 Hz spectral ranges. However, during meditation, prominent spectral components are found only in 0.05\(\)0.1 Hz corresponding to LF band. It is to be mentioned that these two HOSbased plots are obtained by direct FFT technique. Multiple samples of length 256 are taken with an overlap of 50 samples, FFT length is chosen as 512 to estimate the bispectrum and bicoherence. In [84], Garcia et al. have corroborated that reduced HH power in a bispectrum analysis can be an indicator of alterations in vagal activity, which may drastically affect the psychological health by acting on the central autonomic network in response to psychological stress.
Application of waveletbased nonlinear analysis The nonstationarity in HRV signal could easily lead to misleading information if segmentationbased analysis is not conducted. Traditional time domain and spectral analysis methods consider the whole signal to be a stationary one. Besides, many nonlinear methods are suitable to stationary signals only. Under such circumstances, wavelet analysis is welldesigned to analyze such signals. It has the leverage of using a narrow wavelet for analyzing signals demanding high timeresolution and wider wavelet for those requiring highfrequency resolution [85]. A wavelet basically decomposes a signal into scalespecific waveforms by performing convolution operation with the dilated and translated wavelet function, \(\psi (t)\in L^2({\mathbb {R}})\) known as mother wavelet. If a is the scaling parameter and b is the translation parameter, the wavelet transform of a signal x(t) is given by:
Here, \(*\) is the complex conjugate operator. The scaling parameter can be varied to change the width of the analyzing wavelet function and the shift parameter, b facilitates the location specific analysis.
On the other hand, the decomposition of a signal x[n] using discrete wavelet transform (DWT) for a wavelet function (\(\psi [n]\)), is given by:
where \(\phi _{j,k}[n]=2^{j/2}\phi [2^{j}nk]\), \(\psi _{j,k}[n]=2^{j/2}\psi [2^{j}nk]\) are sampled versions of scaling and wavelet functions, respectively; \(j = 0,1,2,....J\) and \(k = 0,1,2,...,2^J1\) with \(j_0\)= 0 being the minimum scale, J the highest scale and N normally selected as integer power of 2. Decomposition of signals into different subbands with unique scales leads to detection of singularity under different scales, which makes it a proficient tool for multifractal analysis. This prompted researchers to analyze multifractality existing in various physical phenomena and contribute remarkable works [86, 87]. Besides, DWT and wavelet packet transform can yield various nonlinear metrics from their coefficients [88, 89].
Based on average wavelet coefficients (AWC), Papasimaki and Pallikari [32] have found that prominent peaks are obtained in the scales ranging from 8–32 beats during meditation. Their results indicates the presence of periodic behavior in LF region of HRV during meditation. Further, they observe a significant drop in Hurst exponent during meditation, which indicates destruction of longrange correlation behavior. In a study by Goshvarpour and Goshvarpour [90], daubechies wavelet ‘db4’ is applied to extract multiple features including the nonlinear parameters such as, skewness, kurtosis and entropy. They have achieved an accuracy of 99.61% in classifying the HRV of meditators (Chi and KYM) and nonmeditators using probabilistic neural network (PNN).
Discussion
We intend to present applications of various nonlinear dynamical approaches for studying HRV dynamics during meditation as reported in the literature over the last 20 years. It has been observed that due to the variation in meditation/relaxation and yoga techniques in terms of breathing procedures, duration and intensity, etc., their impact on ANS is different. To mention, slow breathing meditation and paced breathing (Kapalbhati) meditation adopt completely opposite procedure. Thanks to the PhysioNet that has provided standard datasets for studying the effects of meditation on a single platform [11,12,13, 15, 24, 27, 28, 30, 35, 37, 49]. However, we believe that availability of respiratory and blood pressure data could further enhance the depth of research in studying the impact of meditation. The rationale is that RSA also contributes spectral components at particular time scales which cause HRV modulation. Identically blood pressure also has impacts on HRV. It is to be noted that in case of classificationbased works, exclusivity is a concern as majority of the works are based on binary classifications. Multiclass screening of different groups of subjects (viz., subjects under stress, exercise, meditator, runner, healthy subjects under resting condition) is probably much desired in present scenario.
From the literature it can be inferred that researchers have amply studied the variations in longrange correlation behavior of HRV, dynamical complexity in terms of regularity, and predictability/stability of the HRV signal. It is also noticed that there are some degree of differences in the research outcomes of some works possibly due to the use different datasets or some other reasons. It is observed that during meditation, a majority of researchers demonstrate significant decrease in the complexity of HRV signal by computing nonlinear parameters (entropy, CD, MED, LLE, Poincaré indices, etc.) [12, 13, 24, 26, 27, 29, 33], while a few have observed an increase in dynamical complexity [15, 28, 39]. To examine the regularity and complexity of HRV signals, various forms of entropy parameters are proposed. There are profound interest among researchers to use entropy markers in evaluating HRV signals representing wellbeing and different ailments; a few have worked on to find the most suitable entropy measure for small data size [26, 27, 38]. However, the notion of multiscale analysis is not yet explored pervasively in the studies related to meditation except in [10, 14]. The multiscalebased entropy analysis has revealed that meditative state is rather a complex state unlike that revealed by the entropy analysis at scale 1 [13, 26, 27, 38]. Besides entropy, fractal/multifractal analysis is also performed by many researchers [10, 15, 28, 29, 91]. In these studies, results are found to be inconsistent from one another. On applying VG, multifractality, and FDI techniques to HRV signals for meditative states, a significant increase in the complexity is observed [11, 28, 92]. On a nutshell, it can be stated that nonlinear study demonstrates that HRV possesses more periodic and decimated longrange correlation behavior, low fractal dimension, but higher dynamical complexity during meditation as compared to the nonmeditative state. On the other hand, a prominent decrease in entropy, multifractality and dimensionality is observed in cardiac pathological conditions.
Moving on to other techniques, a new approach based on HRV auditory display is used in studying the impacts of meditation. In this technique, the HRV features (linear and nonlinear) are transformed to auditory signal pitch, timbres, etc., by a process known as sonification. On implementation of machine learningbased discrimination of meditative and premeditative state it has been noticed that only a few works [14, 27, 93, 94] employ machine learning techniques (pattern recognition, SVM, pNN, LVQ, etc., classifiers) to distinguish the meditative state from nonmeditative state, probably due to the small size of publicly available datasets. However data augmentation strategy [14] such as window slicing, concatenation based techniques can be very effective to enable machine learningbased classification works; whereas, in studies related to cardiac pathologies, comparatively larger (adequate) amount data are available in public domain.
Conclusion and future research trends
In this paper, we have carried out a brief review of works demonstrating scientific evidences of the impact of meditation and yoga on HRV dynamics under nonlinear domain. Nonlinear methods are chosen, since HRV is the output of nonlinear interactions of several subsystems. Recent studies revealed that meditation is a complex state in which larger entropy values are found at higher scales. This corroborates the previous findings of higher LLE during meditation. One of the major limitations of the research on meditation and yoga is the nonavailability of sufficient public labeled datasets for experimentation. Besides, given the fact that HRV signal often possess distinct temporal patterns over multiple scales, multiscale entropy/fractal analyses may divulge the hidden heartbeat dynamics in a more robust way. Some of the possible directions for future research in this topic could be:

(1)
No longterm analysis on the effects of different stages of meditation has been performed yet, using HRV signal before, during and after meditation, which may could help in tracking the transformation in ANS functionality over years of continuous practice.

(2)
Evaluation of variations in the dynamical complexity of HRV among experienced meditators, beginners, and control groups by taking sufficient experimental data. A study in this direction could provide more distinct and reliable outcomes.

(3)
To examine whether the changes in heart beat dynamics during meditation is due to controlled respiratory effort (RSA) or due to its notable effect on the ANS. This could be ascertained by using both the HRV and respiratory rate time series for analysis.
Methods
This study is based on the information provided by scientific reports, articles indexed by PubMed, Google Scholar, Web of Science, and Scopus with different keywords, like, “Nonlinear HRV meditation”, “HRV nonlinear analysis yoga”, “HRV and ANS”, “autonomic regulation and nonlinear HRV”.
After preliminary investigation based on the inclusion–exclusion criteria detailed in Fig. 20, we have shortlisted 26 research articles on nonlinear dynamical analysis of the HRV signal. Following this, we present analysis of the nonlinear methods that have gained much attention from the research community, while studying the HRV signal dynamics during meditation. Along with the theoretical background of various nonlinear methods of HRV analysis, we provide our perspectives supported by numerical simulations.
Availability of data and materials
Not applicable.
Abbreviations
 HRV:

Hear rate variability
 ANS:

Autonomic nervous system
 CHF:

Congestive heart failure
 CAD:

Coronary artery disease
 PNS:

Parasympathetic nervous system
 SNS:

Sympathetic nervous system
 SA:

Sinoatrial
 SDNN:

Standard deviation of NN intervals
 SDANN:

SD of average NN intervals
 RMSSD:

Root mean square of successive differences between NN interval
 HTI:

HRV triangular index
 TINN:

Triangular interpolation of NN interval histogram
 ANOVA:

Analysis of variance
 ANCOVA:

Analysis of covariance
 ULF:

Ultra low frequency
 VLF:

Very low frequency
 LF:

Low frequency
 HF:

High frequency
 FFT:

Fast Fourier transform
 PSD:

Power spectral density
 LLE:

Lyapunov exponent
 CD:

Correlation dimension
 NLS:

Nonlinearity score
 ShEn:

Shannon entropy
 CVD:

Cardiovascular disease
 GI:

Guzik’s index
 LOI:

Line of identity
 MSE:

Multiscale entropy
 HOS:

Higher order spectra
 ApEn:

Approximate entropy,
 CCTM:

Component central tendency measures
 DEA:

Diffused entropy analysis
 DFA:

Detrended fluctuation analysis
 DSJE:

Double symbolic joint entropy
 k:

Degree of VG
 LogEn:

Log energy entropy
 LZ:

Lempel–Ziv
 MED:

Minimum embedding dimension
 P(k):

Degree distribution
 PSVG:

Power of scalefreeness in VG
 SampEn:

Sample entropy
References
Hagins M, States R, Selfe T, Innes K. Effectiveness of yoga for hypertension: systematic review and metaanalysis. Evid Based Complement Alternat Med. 2013;2013: 649836.
Li AW, Goldsmith CA. The effects of yoga on anxiety and stress. Altern Med Rev. 2012;17:21–35.
Büssing A, Michalsen A, Khalsa SBS, Telles S, Sherman KJ. Effects of yoga on mental and physical health: a short summary of reviews. Evid Based Complement Alternat Med. 2012;3:1–7.
Innes KE, Bourguignon C, Taylor AG. Risk indices associated with the insulin resistance syndrome, cardiovascular disease, and possible protection with yoga: a systematic review. J Am Board Fam Pract. 2005;18:491–519.
Malik M. Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Eur Heart J. 1996;17:354–81.
Ivanov PC, Amaral LA, Goldberger AL, Havlin S, Rosenblum MG, Struzik ZR, Stanley HE. Multifractality in human heartbeat dynamics. Nature. 1999;399(6735):461–5.
Melillo P, Izzo R, Orrico A, Scala P, Attanasio M, Mirra M, Luca ND, Pecchia L. Automatic prediction of cardiovascular and cerebrovascular events using heart rate variability analysis. PLoS ONE. 2015;10(3):1–14.
Acharya UR, Faust O, Sree V, Swapna G, Martisa RJ, Kadri NA, Suri JS. Linear and nonlinear analysis of normal and cadaffected heart rate signals. Comput Meth Prog Bio. 2014;113:55–68.
Kamath MV, Fallen EL. Power spectral analysis of heart rate variability: a noninvasive signature of cardiac autonomic function. Crit Revs Biomed Eng. 1993;21:245–311.
Sarkar A, Barat P. Effect of meditation on scaling behavior and complexity of human heart rate variability. Fractals. 2008;16(3):199–208.
Peng CK, Mietus JE, Liu Y, Khalsa G, Douglas PS, Benson H, Goldberger AL. Exaggerated heart rate oscillations during two meditation techniques. Int J Cardiol. 1999;70(2):101–7.
Raghavendra BS, Dutt DN. Multiscale Fractal Dimension Technique for Characterization and Analysis of Biomedical Signals. In: 2011 Digital Signal Processing and Signal Processing Education Meeting (DSP/SPE), Sedona, USA, pp. 2011:370–374.
Goswami DP, Bhattacharya DK, Tibarewala DN. Analysis of heart rate variability in meditation using normalized Shannon entropy. J Int Acad Phys Sci. 2010;14(1):61–7.
Deka D, Deka B. Characterization of heart rate variability signal for distinction of meditative and premeditative states. Biomed Signal Process Control. 2021;66:1–13.
Kamath C. Analysis of heart rate variability signal during meditation using deterministicchaotic quantifiers. J Med Eng Technol. 2013;37(7):436–48.
Wang Y, Wei S, Zhang S, Zhanga Y, Zhaoa L, Liu C, Murray A. Comparison of timedomain, frequencydomain and nonlinear analysis for distinguishing congestive heart failure patients from normal sinus rhythm subjects. Biomed Signal Process Control. 2018;42:30–6.
Neves VR, Takahashi ACM, do SantosHiss MDB, Kiviniemi AM, Tulppo MP, de Moura SCG, Karsten M, BorghiSilva A, Porta A, Montano N, Catai AM. Linear and nonlinear analysis of heart rate variability in coronary disease. Clin Auton Res. 2012;22:175–83.
Cohen H, Benjamin J, Gevab AB, Matar MA, Kaplan Z, Kotler M. Autonomic dysregulation in panic disorder and in posttraumatic stress disorder: application of power spectrum analysis of heart rate variability at rest and in response to recollection of trauma or panic attacks. Psychiatry Res. 2000;96:1–13.
Acharya UR, Joseph KP, Kannathal N, Min LC, Suri JS. Heart rate variability. In: Acharya UR, Suri JS, Spaan JAE, Krishnan SM, editors. Advances in cardiac signal processing. Berlin, Heidelberg: Springer; 2007. p. 121–65.
Murata T, Takahashi T, Hamada T, Omori M, Kosaka H, Yoshida H, Wada Y. Individual trait anxiety levels characterizing the properties of Zen meditation. Neuropsychobiology. 2004;50:189–94.
Kheder G, Taleb R, Kachouri A, Messaoud MB, Samet M. Feature Extraction by Wavelet Transforms to Analyze the Heart Rate Variability During Two Meditation Technique. In: 6th WSEAS Int. Conf. on Circuits, Systems, Electronics, Control & Signal Processing, Cairo, Egypt, pp. 2007:374–378
Jovanov E. On spectral analysis of heart rate variability during very slow yogic breathing. In: Conf. Proc. IEEE Eng. Med. Biol. Soc., pp. 2005:2467–2470.
Peressutti C, MartinGonzalez JM, Garcia JM, Manso, Mesa D. Heart rate dynamics in different levels of Zen meditation. Int J Cardiol. 2009;145(1):142–6.
Raghavendra BS, Dutt DN. Nonlinear dynamical characterization of heart rate variability time series of meditation. Int J Biomed Biol Eng. 2011;5(9):429–40.
Goshvarpour A, Goshvarpour A. Chaotic behavior of heart rate signals during Chi and Kundalini meditation. Int J Image Graph Signal Process. 2012;2:23–9.
Li J, Hu J, Zhang Y, Wang J, Zhang X. Dynamical complexity changes during two forms of meditation. Physica A. 2011;390(12):2381–7.
Goshvarpour A, Goshvarpour A. Do meditators and nonmeditators have different HRV dynamics? Cogn Syst Res. 2019;54:21–36.
Bhaduri A, Ghosh D. Quantitative assessment of heart rate dynamics during meditation: an ECG based study with multifractality and visibility graph. Front Physiol. 2016;7(44):1–10.
Diosdado AM, Coyt GG, Uribe BMP. Oscillations in the Evaluation of Fractal Dimension of RR Intervals Time Series. In: 32nd Annual Int. Conf. of the IEEE EMBS Buenos Aires, Argentina, pp. 2010:4570–4573.
Jiang S, Bian C, Ning X, Ma QDY. Visibility graph analysis on heartbeat dynamics of meditation training. Appl Phys Lett. 2013;102(25):1–3.
Gao ZK, Cai Q, Yang YX, Dang WD. Timedependent limited penetrable visibility graph analysis of nonstationary time series. Physica A. 2017;476:43–8.
Papasimakis N, Pallikari F. Breakdown of LongRange Correlations in Heart Rate Fluctuations During Meditation
Goswami DP, Tibarewala DN, Bhattacharya DK. Analysis of heart rate variability signal in meditation using secondorder difference plot. J Appl Phys. 2011;109(114703):1–6.
Goshvarpour A, Goshvarpour A. Comparison of higher order spectra in heart rate signals during two techniques of meditation: Chi and kundalini meditation. Cogn Neurodyn. 2013;7:39–46.
Goshvarpour A, Goshvarpour A. Poincaré indices for analyzing meditative heart rate signals. Biomed J. 2015;38(3):229–34.
AlvarezRamirez J, Rodriguez E, Echeverria JC. Fractal scaling behavior of heart rate variability in response to meditation techniques. Chaos Soliton Fract. 2017;99:57–62.
Goshvarpour A, Goshvarpour A. A novel feature level fusion for heart rate variability classification using correntropy and CauchySchwarz divergence. J Med Syst. 2018;42(6):1–15.
Yao W, Zhang Y, Wang J. Quantitative analysis in nonlinear dynamic complexity detection of meditative heart beats. Physica A. 2018;512:1060–8.
Guo M, Guo X, Wang M, et al. Activation of Sympathetic Nervous System as a Biomarker for Deep Meditation. In: 9th International IEEE/EMBS Conf. on Neural Engineering (NER), 2019:546–549. IEEE, San Francisco, USA . https://doi.org/10.1109/NER.2019.8717006
Nasrolahzadeh M, Mohammadpoory Z, Haddadnia J. Analysis of heart rate signals during meditation using visibility graph complexity. Cogn Neurodyn. 2019;13:45–52.
Deka D, Deka B. Investigation on hrv signal dynamics for meditative intervention. In: et al., M.P. (ed.) Soft Computing: Theories and Applications. Advances in Intelligent Systems and Computing, 2020;1154:993–1005. Springer, Singapore
Goshvarpour A, Goshvarpour A. Asymmetry of lagged Poincaré plot in heart rate signals during meditation. J Tradit Complement Med. 2020;11(1):16–21.
Rohila A, Sharma A. Asymmetric spread of heart rate variability. Biomed Signal Process Control. 2020;60: 101985.
Goshvarpour A, Goshvarpour A. Verhulst map measures: new biomarkers for heart rate classification. Phys Eng Sci Med. 2022;45(2):513–23.
Deka B, Deka D. An improved multiscale distribution entropy for analyzing complexity of realworld signals. Chaos Soliton Fract. 2022;158:1–9.
Kamen PW, Krum H, Tonkin AM. Poincare plots of heart rate variability allows quantitative display of parasympathetic nervous activity in humans. Clin Sci. 1996;91:201–8.
Tulppo MP, Makikallio TH, Takala TES, Seppanen T, Huikuri HV. Quantitative beattobeat analysis of heart rate dynamics during exercise. Am J Physiol. 1996;271:244–52.
Karmakar CK, Khandoker AH, Palaniswami M. Investigating the changes in heart rate asymmetry (HRA) with perturbation of parasympathetic nervous system. Australas Phys Eng Sci Med. 2012;35:465–74.
Dey A, Bhattacharya DK, Tibarewala DN, Dey N, Ashour AS, Le DN, Gospodinova E, Gospodinov M. ChineseChi and Kundalini yoga meditations effects on the autonomic nervous system: comparative study. Int J Interact Multimed Artif Intell. 2016;3(7):87–95.
Cao L. Practical method for determining the minimum embedding dimension of a scalar time series. Physica D. 1997;110:43–50.
Takens F. Detecting strange attractors in turbulence. In: Rand D, Young L, editors. Dynamical systems and turbulence. Lecture notes in mathematics. Warwick: Springer; 1981. p. 366–81.
Kennel MB, Brown R, Abarbanel HDI. Determining embedding dimension for phasespace reconstruction using a geometrical construction. Phys Rev A. 1992;45(6):3403–11.
Broomhead DS, King GP. Extracting qualitative dynamics from experimental data. Physica D. 1986;20(2–3):217–36.
Grassberger P, Procaccia I. Characterization of strange attractors. Phys Rev Lett. 1983;50(5):346–9.
Goshvarpour A, Goshvarpour A. Recurrence plots of heart rate signals during meditation. I J Image Graph Signal Process. 2012;2:44–50.
Zurek S, Guzik P, Pawlak S, Kosmidera M, Piskorski J. On the relation between correlation dimension, approximate entropy and sample entropy parameters, and a fast algorithm for their calculation. Physica A. 2012;391:6601–10.
Bogaert C, Beckers F, Ramaekers D, Aubert AE. Analysis of heart rate variability with correlation dimension method in a normal population and in heart transplant patients. Auton Neurosci. 2001;90:142–7.
Beckers F, Verheyden B, Aubert AE. Aging and nonlinear heart rate control in a healthy population. Am J Physiol Heart Circ Physiol. 2006;290:2560–70.
Kantz H, Schreiber T. Nonlinear time series analysis. Cambridge: Cambridge University Press; 1997.
Valenza G, Allegrini P, Lanata A, Scilingo EP. Dominant Lyapunov exponent and approximate entropy in heart rate variability during emotional visual elicitation. Front Neuroeng. 2012;5:1–7.
Casties JF, Mottet D, Gallais DL. Nonlinear analyses of heart rate variability during heavy exercise and recovery in cyclists. Int J Sports Med. 2006;27(10):780–5.
Eckmann JP, Kamphorst SO, Ruell D. Recurrence plots of dynamical systems. Europhys Lett. 1987;4(91):973–7.
Jr Webber CL, Zbilut JP. Dynamical assessment of physiological systems and states using recurrence plot strategies. J Appl Physiol. 1994;76(2):965–73.
Marwan N, Wessel N, Meyerfeldt U, Schirdewan A, Kurths J. Recurrenceplotbased measures of complexity and their application to heartratevariability data. Phys Rev E. 2002;E66:026702–18.
Costa M, Goldberger AL, Peng CK. Multiscale entropy analysis of complex physiologic time series. Phys Rev Lett. 2002;88(6):068102–14.
Bandt C, Pompe B. Permutation entropy: a natural complexity measure for time series. Phys Rev Lett. 2002;88(17):174102–14.
Pincus SM. Approximate entropy as a measure of system complexity. Proc Natl Acad Sci. 1991;88:2297–301.
Richman JS, Moorman JR. Physiological timeseries analysis using approximate entropy and sample entropy. Am J Physiol Heart Circ Physiol. 2000;278(6):2039–49.
Li J, Ning X. Dynamical complexity detection in shortterm physiological series using basescale entropy. Phy Rev E. 2006;73:052902–14.
Liu X, Jiang A, Xu N, Xue J. Increment entropy as a measure of complexity for time series. Entropy. 2016;18(22):18010022–114.
Tsallis C. Possible generalization of BoltzmannGibbs statistics. J Stat Phys. 1988;52(1–2):479–87.
Porto AA, Tavares BS, Vidigal G, Garner DM, Raimundo RD, de Abreu LC, Bocalini D, Valenti VE. Nonlinear dynamics of heart rate during slow breathing exercise. Indian J Physiol Pharmacol. 2018;62(2):160–9.
Peng CK, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos. 1995;5(1):82–7.
Eke A, Herman P, Bassingthwaighte J, Raymond G, Percival DB, Cannon MIB, Ikrenyi C. Physiological time series: distinguishing fractal noises from motions. Pflugers Archiv. 2000;439:403–15.
Hardstone R, Poil SS, Schiavone G, Jansen R, Nikulin VV, Mansvelder HD, LinkenkaerHansen K. Detrended fluctuation analysis: a scalefree view on neuronal oscillations. Front Physiol. 2012;3:450.
Gronwald T, Ludyga S, Hoos O, Hottenrott K. Nonlinear dynamics of cardiac autonomic activity during cycling exercise with varied cadence. HMSC. 2018;60:225–33.
Kantelhardt JW, Zschiegner SA, KoscielnyBunde E, Havlind S, Bunde A, Stanley HE. Multifractal detrended fluctuation analysis of nonstationary time series. Physica A. 2002;316:87–114.
Song R, Bian C, Qianli MDY. Multifractal analysis of heartbeat dynamics during meditation training. Physica A. 2013;392:1858–62.
Kim J, Wilhelm T. What is a complex graph? Physica A. 2008;387:2637–52.
Nikias CL, Raghuveer MR. Bispectrum estimation: a digital signal processing framework. Proc IEEE. 1987;75:869–91.
Mohebbi M, Ghassemian H. Prediction of paroxysmal atrial fibrillation based on nonlinear analysis and spectrum and bispectrum features of the heart rate variability signal. Comput Meth Prog Bio. 2012;105:40–9.
Yu SN, Lee MY. Bispectral analysis and genetic algorithm for congestive heart failure recognition based on heart rate variability. Comput Biol Med. 2012;42:816–25.
Chua KC, Chandran V, Acharya UR, Lim CM. Cardiac state diagnosis using higher order spectra of heart rate variability. Journal Med Eng Technol. 2008;32(2):145–55.
Garcia RG, Valenza G, Tomaz CA, Barbieri R. Instantaneous bispectral analysis of heartbeat dynamics for the assessment of major depression. In: 2015 Computing in Cardiology Conference (CinC), 2015:781–784 .
Mallat S. A wavelet tour of signal processing: the sparse way. 3rd ed. San Diego: Academic Press; 1998.
Abry P, Veitch D. Wavelet analysis of longrangedependent traffic. IEEE Trans Inf Theory. 1998;44(1):2–15.
Muzy JF, Bacry E, Arneodo A. The multifractal formalism revisited with wavelets. Int J Bifurcat Chaos. 1994;4(2):245–302.
Acharya UR, Sudarshan VK, Ghista D, Lim WJE, Molinari F, Sankaranarayanan M. Computeraided diagnosis of diabetic subjects by heart rate variability signals using discrete wavelet transform method. Knowl Based Syst. 2015;81:56–64.
Avci E. Comparison of wavelet families for texture classification by using wavelet packet entropy adaptive network based fuzzy inference system. Appl Soft Comput. 2008;8(1):225–31.
Goshvarpour A, Goshvarpour A. Matching pursuit based indices for examining physiological differences of meditators and nonmeditators: an HRV study. Physica A. 2019;524:147–56.
Deka D, Deka B. Computeraided detection of congestive heart failure based on nonlinear hrv features. In: 2020 IEEE International Conference on Machine Learning and Applied Network Technologies, Hyderabad, 2020:1–6.
Gonzalez MAG, Castro JJR, FernandezChimeno M. A Novel Index Based on Fractional Calculus to Assess the Dynamics of Heart Rate Variability: Changes Due to Chi or Yoga Meditations. In: 2012 Computing in Cardiology, Krakow, Poland, 2012;39:925–928.
Phongsuphap S, Pongsupap Y, Chandanamattha P, Lursinsap C. Changes in heart rate variability during concentration meditation. Int J Cardiol. 2008;130:481–4.
Lee YH, Chen SCJ, Shiah YJ, Wang SF, Young MS, Hsu CY, Cheng GQJ, Lin CL. Supportvectormachinebased meditation experience evaluation using electroencephalography signals. J Med Biol Eng. 2014;34(6):589–97.
Acknowledgements
Authors would like to acknowledge AICTE, New Delhi, for providing financial support under the QIP scheme and Tezpur University for providing the necessary infrastructure to carry out this research.
Funding
Authors would like to acknowledge AICTE, New Delhi, for providing financial support under the QIP scheme.
Author information
Authors and Affiliations
Contributions
BD: design of the review work; BD and DD: methodology; BD: formal analysis; DD: simulation and investigation; DD: writing—original draft; BD: writing—review. Both authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
The authors declare that they consent for publication.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Deka, B., Deka, D. Nonlinear analysis of heart rate variability signals in meditative state: a review and perspective. BioMed Eng OnLine 22, 35 (2023). https://doi.org/10.1186/s12938023011003
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12938023011003