Model For The Respiratory Modulation Of The Heart

ARTICLE IN PRESS

Physica A 355 (2005) 439–460 www.elsevier.com/locate/physa

Model for the respiratory modulation of the heart beat-to-beat time interval series Alberto Capurroa,, Luis Diambrab, C.P. Maltaa a

Instituto de Fı´sica, Universidade de Sa˜o Paulo, Rua do Mata˜o, Travessa R 187, CEP 05508-900 Sa˜o Paulo, SP, Brazil b Instituto de Cieˆncias Biome´dicas, Universidade de Sa˜o Paulo, Av. Lineu Prestes 1524, CEP 05508-900 Sa˜o Paulo, SP, Brazil Received 16 January 2004 Available online 18 April 2005

Abstract In this study we present a model for the respiratory modulation of the heart beat-to-beat interval series. The model consists of a set of differential equations used to simulate the membrane potential of a single rabbit sinoatrial node cell, excited with a periodic input signal with added correlated noise. This signal, which simulates the input from the autonomous nervous system to the sinoatrial node, was included in the pacemaker equations as a modulation of the iNaK current pump and the potassium current iK . We focus at modeling the heart beat-to-beat time interval series from normal subjects during meditation of the Kundalini Yoga and Chi techniques. The analysis of the experimental data indicates that while the embedding of pre-meditation and control cases have a roughly circular shape, it acquires a polygonal shape during meditation, triangular for the Kundalini Yoga data and quadrangular in the case of Chi data. The model was used to assess the waveshape of the respiratory signals needed to reproduce the trajectory of the experimental data in the phase space. The embedding of the Chi data could be reproduced using a periodic signal obtained by smoothing a square wave. In the case of Kundalini Yoga data, the embedding was reproduced with a periodic signal obtained by smoothing a triangular wave having a rising branch of longer duration than

Corresponding author. Tel.: +55 11 3091 6976; fax: +55 11 3091 6833.

E-mail address: [email protected] (A. Capurro). 0378-4371/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2004.04.135

ARTICLE IN PRESS 440

A. Capurro et al. / Physica A 355 (2005) 439–460

the decreasing branch. Our study provides an estimation of the respiratory signal using only the heart beat-to-beat time interval series. r 2005 Elsevier B.V. All rights reserved. Keywords: Respiratory modulation; Pacemaker; Heart rate variability; Rabbit SA node; Meditation

1. Introduction Cells within the sinoatrial (SA) node comprise the primary pacemaker site within the heart. These cells are characterized as having no true resting potential, but instead generate regular spontaneous action potentials. The membrane slowly depolarizes during the diastole of the cardiac contraction cycle leading to the onset of an action potential. The repolarization, that occurs after the spike, activates the so called slow diastolic depolarization again and initiates a new cycle. For a comprehensive review of the ionic currents underlying this membrane voltage oscillation see [1]. The intrinsic automaticity (spontaneous pacemaker activity) of the SA node determines a spiking rate of 100–110 action potentials per minute. This intrinsic rhythm is modulated by autonomic nerves, with parasympathetic (vagal) influence being dominant over sympathetic influence at rest. The ‘‘vagal tone’’ brings the resting heart rate down to 60–80 beats/min. Sympathetic activation increases the pacemaker rate with concomitant inhibition of vagal tone [2]. The acetylcholine (ACh) released by the vagal nerve binds to the muscarinic receptors at the post-synaptic membrane, and induces three main actions [1,3,4]: (i) through a reduction in the cyclic adenosine monophosphate (AMPc) production it inhibits the cationic current iNaK , thus decreasing the slope of the slow diastolic depolarization and slowing down the pacemaker rate without changing the membrane potential, (ii) it opens a G-protein coupled potassium channel that hyperpolarizes the cell and slows down the pacemaker rate, and (iii) it inhibits the long lasting component of the calcium current iCa . The first mechanism is opposite to the b adrenergic action of the sympathetic system that increases the AMPc production [5] and the slope of the slow diastolic depolarization. The sympathetic action also activates the iCa , but this current has little or no effect on the heart rate. Experimental evidence indicates that the main process that mediates the effect of low-concentration ACh actions on the heart rate is the inhibition of the cationic current iNaK [6]. The normal respiratory cycle is accompanied by changes in autonomic tone that modulates the heart rate. The activity of the vagal nerve endings increases during exhalation, and the activity of sympathetic fibers increases during inhalation, causing the ‘‘respiratory modulation’’ (RM) or ‘‘sinus arrhythmia’’, i.e. during inhalation the heartbeat intervals shorten and during exhalation they stretch. The oscillation in vagal action is responsible for most of the RM, because it is faster than the sympathetic action [2]. At typical respiratory frequencies (greater than

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

441

0.15 Hz), increase of the heart rate occurs simultaneously with the onset of inspiratory activity, while at frequencies less than 0.15 Hz, increase of the heart rate precedes inspiration [7]. The RM is reflected in the high-frequency power band of the heart rate variability (HRV), which is abolished by atropine infusion, implying that it is parasympathetically mediated [8]. It has been recently reported that the RM may improve the efficiency of pulmonary gas exchange in healthy humans [9]. The RM has larger amplitude in children and young subjects than in elders. In various pathological conditions the HRV is decreased, and the RM is partially or completely absent [8,10,11]. The chaotic structure found in normal subjects [12] tends to disappear in pathological states [13]. The heart rate dynamics was also studied during a variety of physiological conditions [14–17]. During Kundalini Yoga and Chi meditation practices in normal subjects, there have been reported extremely prominent heart rate oscillations correlated with slow breathing [18], with an amplitude significantly larger than the RM measured in the same individuals before meditation. In this study we present a model for the influence of respiration on the heart beatto-beat interval series, focusing on heart beat data from healthy subjects during meditation sessions of Kundalini Yoga and Chi techniques. We first show that the embedding of the experimental beat-to-beat interval data acquires a polygonal shape during meditation, and that the time course of respiration depends on the meditation technique. Then, we use the model to assess the waveshape of the RM needed to reproduce the trajectory of the experimental beat-to-beat interval series in the phase space.

2. Methods 2.1. Meditation protocols Two specific meditation techniques were studied: (i) Chinese Chi (or Qigong) meditation (as taught by Xin Yan) and (ii) Kundalini Yoga meditation (as taught by Yogi Bhajan). In this subsection we reproduce a brief description of the meditation protocols followed [18]. The eight Chi meditators, 5 women and 3 men (age range 26–35 years), wore a Holter recorder for approximately 10 h during which time they went about their ordinary daily activities. After approximately 5 h of recording, each one of them practiced 1 h of meditation. Meditation beginning and ending time were delineated with event marks. During these sessions, the Chi meditators sat quietly, listening to the taped guidance of the master. The meditators were instructed to breath spontaneously while visualizing the opening and closing of a perfect lotus on the abdomen. The four Kundalini Yoga meditators, 2 women and 2 men (age range 20–52 years), wore a Holter monitor for approximately one and half hours, 15 min of baseline quiet breathing were recorded before the meditation period of 1 h. The

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

442

meditation protocol consisted of a sequence of breathing and chanting exercises, performed while seated in a cross-legged posture. The beginning and ending of the various meditation sub-phases were delineated with event marks. The Kundalini Yoga subjects were considered to be at an advanced level of meditation training. In the case of Chi meditation, the subjects were relative novices, having begun their meditation practice about 1–3 months before the recordings used for this study. The persons from both groups were in good general health conditions. 2.2. Data files We used time series, from physionet [19], of 4 subjects prior and during Kundalini Yoga meditation, 8 subjects prior and during Chi meditation, 4 healthy control subjects spontaneously breathing during sleeping hours, and 4 healthy subjects during supine metronomic breathing at 0.25 Hz. The heart beat-to-beat intervals were obtained from electrocardiogram (ECG) data by Peng et al. [18] and made available in physionet. The data files (in text format) consisted of two columns, the first containing the heart beat instant of time, and the second containing the instantaneous heartbeat rate. The sequence of beat-to-beat time intervals was generated using the data in the first column. From this sequence we selected segments of 300 intervals that had no trends or intervals of double length due to undetected heartbeats. From the selected segments of intervals we determined the autocorrelation function and a bidimensional embedding [20]. The value of the embedding delay was determined by the first zero of the autocorrelation function. 2.3. Model for the RM of the heart pacemaker 2.3.1. Pacemaker model Several models of the SA node pacemaking in the rabbit heart have been developed in the last years [21–25]. We selected the model in Endresen et al. [23] because it provides a realistic and simple representation of the membrane potential oscillations. This model, based on experimental in vitro studies of isolated rabbit SA node cells [1,26–29], was developed to fit experimental observations of (i) the action potential shape and (ii) the intracellular ionic concentrations of potassium, calcium and sodium. The model takes into account the ionic currents iK ; iCa ; iNa ; iNaK , and iNaCa . The membrane voltage v is given by the net charge difference between the intracellular and extracellular compartments as v¼

FV f½Ki þ 2½Cai þ ½Nai  ½Ke  2½Cae  ½Nae g , C

(1)

where ½ i ð½ e ) denotes the intracellular (extracellular) ionic concentration, C is the membrane capacitance, F is Faraday’s constant and V is the cell volume.

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

443

The intracellular ionic concentrations ½Ki , ½Cai and ½Nai evolve according to d 2iNaK  iK ½Ki ¼ , dt FV

(2)

d 2iNaCa  iCa ½Cai ¼ , dt 2FV

(3)

d iNa  3iNaK  3iNaCa ½Nai ¼ . dt FV The currents are given by   ev 1 ½Ke iK ¼ kK x sinh  ln , 2kT 2 ½Ki 

iCa

ev 1 ½Cae  ln ¼ ðkCa fd 1 þ kb;Ca Þ sinh kT 2 ½Cai

(4)

(5)  ,

 ev 1 ½Nae  ln iNa ¼ kNa hm1 sinh , 2kT 2 ½Nai   eðv  vATP Þ ½Ke 3 ½Nae þ ln iNaK ¼ kNaK tanh  ln , 2kT ½Ki 2 ½Nai   ev 3 ½Nae ½Cae  ln þ ln iNaCa ¼ kNaCa sinh , 2kT 2 ½Nai ½Cai

(6)



(7)

(8)

(9)

where kK , kCa , kNa , kNaK and kNaCa are the conductance parameters corresponding to each current (kb;Ca is the conductance of the background calcium current). vATP is the half activation potential of the NaK pump (ATP denotes cyclic adenosine triphosphate), k is the Bolztmann constant, e is the elementary charge, and T is the absolute temperature. The activation mechanism of calcium and sodium channels is very fast, so we use the steady state fraction of open channels (denoted by d 1 for the calcium channels, and m1 for the sodium channels) given by    1 2eðv  vd Þ d 1 ¼ 1 þ tanh , (10) 2 kT    1 2eðv  vm Þ m1 ¼ 1 þ tanh , (11) 2 kT where vd and vm are the half activation potentials of the calcium and sodium channels, respectively. The fraction of open potassium channels, x, evolves as        dx 1 2eðv  vx Þ 1 2eðv  vx Þ ¼ 1 þ tanh cosh  x , (12) dt tK kT 2 kT

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

444

where vx is the half activation potential of the potassium channel and tK is the relaxation time of the potassium channel. In addition to the activation mechanism, the calcium and sodium channels have an inactivation mechanism. The inactivation rates of the calcium and sodium channels are given by        2eðv  vf Þ 2eðv  vf Þ df 1 1 ¼ 1  tanh cosh  f , (13) dt tCa 2 kT kT        dh 1 2eðv  vh Þ 1 2eðv  vh Þ ¼ 1  tanh cosh  h , dt tNa kT 2 kT

(14)

where vf and vh are the half inactivation potentials of the calcium and sodium channels, respectively. tCa and tNa are the relaxation time of the inactivation process of the calcium and sodium channels, respectively. As the values of tK and vx are very close to the analogous magnitudes of tCa and vf , respectively, it is possible to relate the inactivation gating of calcium to the activation gating of potassium as x þ f ¼ 1. This allows to reduce the number of differential equations in the model [23]. Eqs. (12) and (13) are identical but have opposite dynamics. So, we have suppressed Eq. (13), and replaced f by 1  x in Eq. (6). 2.3.2. Pacemaker modulation The model for the rabbit SA node presented above [23] constitutes a plausible approximation to the human SA node, which is the source of the beat-to-beat interval data used in the present study. Our model for the RM of the human SA node was obtained by introducing a modulation in the sodium–potassium conductance parameter kNaK , and in the relaxation time of the potassium current tK in Eqs. (8) and (12) of the model for the rabbit SA node. This was accomplished by setting kNaK ðtÞ ¼ kNaK IðtÞ ,

(15)

and tK ðtÞ ¼

tK , 5:0IðtÞ

(16)

with IðtÞ being a time dependent input signal. By replacing kNaK by kNaK ðtÞ in Eq. (8), the current iNaK is modulated by IðtÞ. The current iK is modulated by IðtÞ in Eq. (5) by replacing tK by tK ðtÞ in Eq. (12). The reason for using the input signal multiplied by 5.0 in (16) is explained later in this section. The currents iNaK and iK where selected to be the target of regulation in the model because they mediate the main action of the autonomous nervous system on the heart rate [1]. As the inactivation gating of calcium channels was related to the activation gating of potassium channels as x þ f ¼ 1, we also included in the model a regulation of iCa , i.e. the relaxation time of the calcium inactivation is modulated in the same way as the relaxation time of the potassium activation. This last

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

445

modulation was included in the model because it is also described in the literature [1], although it may be not relevant for the heart rate. The input IðtÞ was formed by the sum of two time dependent signals (SðtÞ and NðtÞ) and the constant term A. SðtÞ is a periodic wave representing the influence of respiration. NðtÞ is a realization of correlated noise addressed to the variability of the breathing rhythm around a periodic pattern and to other influences coming through the autonomic innervation of the SA node. The constant term A is inversely related to the tonic activity of the colinergic input to the SA node (i.e. the vagal tone is represented by A). The sum SðtÞ þ NðtÞ evolves as a mirror image of the oscillatory component of the vagal action and parallel to the respiratory cycle. In this way, the inputs to the pacemaker model affect the main currents that are regulated by the action of the autonomous nervous system [1]. The manner in which the influences are introduced (Eqs. (15) and (16)) is arbitrary but physilologically plausible. In a previous study [30], we simulated the RM by injecting an external current in the membrane of a simpler model (Hodgkin–Huxley equations), without including direct modulations of the ionic currents. Indeed, the trajectory of the beatto-beat intervals in the phase space can be qualitatively reproduced by a general type of oscillator, provided the representation of the input signal in the output spike train is close to linear (see Section 3.1). With the present model, in which the input signals modulate specific ionic currents of the pacemaker, we could reproduce the results also quantitatively. The modulation of specific ionic currents is a necessary condition for this success, since we checked that the embedding shape of the data could not be reproduced when the input signals to the rabbit pacemaker model consisted of external current. In Section 3.1 it will be shown that, in order to reproduce the experimental results, the input signal has to be quasi-linearly represented in the output spike train. In order to achieve this condition kNaK and kK were set to 15.0 and 120.0 pA, respectively. The other parameter values were set as provided in Endresen et al. [23]. The values used for the parameters, initial conditions, observed constants, and fundamental physical constants are given in the Table 1. We integrated the system using a forward Euler algorithm (using C codes). At first the model was run excited only by the constant term A (with SðtÞ and NðtÞ set to zero) until a constant beat rate was reached. Then, the time dependent modulations were activated and the actual simulation started. Each realization lasted 106 sample steps, that corresponds to 200 s (sample step ¼ 0:2 ms). This was enough to obtain 200–300 simulated heart beat-to-beat intervals, using a spike detection algorithm explained at the end of this section. We checked that the time course followed by the membrane voltage and the ionic currents of our simulation are similar to those shown in Endresen et al. [23]. The Fig. 1A and B shows an example of the membrane potential oscillation and the currents iK ; iCa ; iNa ; iNaK , and iNaCa (the time dependent modulations were set to zero). The sense of the current is always positive (i.e. outward sense) for iNaK and iK . When the time dependent modulation SðtÞ þ NðtÞ is activated, we observe changes in these currents, which are related to the cycle of SðtÞ. The oscillation of iNaK in each action potential is shifted in the inward sense (i.e. to lower values) when the intervals

ARTICLE IN PRESS 446

A. Capurro et al. / Physica A 355 (2005) 439–460

Table 1 Values used for the parameters, initial conditions, observed constants, and fundamental physical constants Adjustable parameters

Physical constants

Name

Value

Name

Value

kCa kb;Ca kNa kNaCa kK kNaK

26.2 pA 0.01645 pA 112.7 pA 1400.0 pA 120:0 pA 15:0 pA

k e F R

1:38  1020 mJ K1 1:6  1019 C 96485:3 C mol1 8314:5 J kmol1 K1

Initial conditions

Observed constants Name

Value

T V

310.15 K 10  103 mm3

Name

Value

½Ke

5.4 mM

x0 h0 ½Ki0 ½Cai0 ½Nai0

0.1 0.008 130.66 mM 0.0006 mM 18.73 mM

½Cae ½Nae C vx ¼ vf vd vm vh vATP tK ¼ tCa ¼ tNa

2 mM 140 mM 47 pF 25.1 mV 6.6 mV 41.4 mV 91 mV 450 mV 200 ms

The asterisks mark the parameters that are different from the value provided in Endresen et al. (2000).

shorten to its minimum value (Fig. 1C and D), corresponding to the peak of SðtÞ. In the case of iK , the maximum value does not change with SðtÞ, but its duration increases at the valley of SðtÞ when the intervals reach their larger duration (not shown). The difference in the membrane potential from the peak to the valley of SðtÞ is only 2.5 mV (not shown) suggesting that the dominant mechanism of interval modulation is provided by the regulation of iNaK , as seems to be the case in the SA node cells [6,1]. In order to estimate the relative contribution of iNaK and iK to the interval modulation of our model, we focused at the maximum interval modulation amplitude induced by a typical RM oscillation of the Yoga data simulation. The oscillation amplitude reached 0.28 s when the two currents received the modulation, and was reduced to 0.24 s when we selectively eliminated the effect on iK . When we eliminated selectively the effect on iNaK , the amplitude was reduced to only 0.04 s (not shown). Similar results were obtained using a RM oscillation of the Chi data simulation. Thus, as occurs in the rabbit SA node [1], the main effect of interval modulation in our model is mediated by the iNaK , while the modulation of iK is needed to reproduce the heart rate of the different data sets. To set our model in this

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460 0.8

Beat-to-beat interval (s)

20

v (mV)

0 -20 -40 -60 -80

447

0

0.4

(A)

0.8

1.2

0.7

0.6

0.5

1.6

0

4

(C)

Time (s)

8

12

16

20

16

20

Time (s)

200

15

iK

iNaK (pA)

i (pA)

100

i NaK 0

i NaCa

iNa

13

11

-100

i Ca -200

(B)

0

0.4

0.8

Time (s)

1.2

9

1.6

(D)

0

4

8

12

Time (s)

Fig. 1. (A) Membrane potential v vs. time (A ¼ 0:01, SðtÞ ¼ 0:0 and NðtÞ ¼ 0:0). (B) Ionic currents vs. time (corresponding to the membrane potential shown in the panel A): iNaK (full thick line), iK (full thin line), iCa (dashed line), iNaCa (dot-dashed line) and iNa (dotted line). (C) Beat-to-beat interval vs. time. The input modulation was the last cycle of IðtÞ, with the parameters and waveshape used for the simulation of the Kundalini Yoga data (Table 2), with the exceptions of SðtÞ amplitude and period, here set to 0.033 and 16.7 s, respectively. (D) Ionic current iNaK vs. time (corresponding to the beat-to-beat interval plot shown in the panel C). Note that the oscillation of iNaK is shifted in the direction of the zero value when the beatto-beat intervals reach their lower duration (higher beat rate), and vise-versa.

regime, the input to kNaK was IðtÞ, while the input to tK was 5:0IðtÞ (Eqs. (15) and (16)). The analysis of Eqs. (15) and (16) shows that the input signals are dimensionless. The methods to obtain these signals are presented in the following. In order to obtain NðtÞ, realizations lasting 106 sample steps of the following random differential equation (x is a white Gaussian noise), dN ¼ l1 N þ l2 x , dt

(17)

were low pass filtered with a Hanning window (half width ¼ 0:06 s). The equation was numerically solved using a forward Euler algorithm. The correlation level of

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

448

NðtÞ, controlled through the parameter l1 [31], was set to 0.04 (l2 was set to 1.0). The mean value of NðtÞ tends to zero, but in order to ensure an exact zero mean in realizations with a limited number of points, we added or subtracted an appropriate value so as to drive the signal mean value to zero. The standard deviation (SD) of NðtÞ was used as a measurement of the signal amplitude. It was initially set to 1.0, by multiplying each value by the appropriate factor. In this way, the final SD of the signal actually input into the model could be easily controlled by multiplying the time series by the SD value desired in each simulation. NðtÞ was produced with a code in C and fed into the pacemaker code. For the simulation of the Kundalini Yoga data, the periodic signal SðtÞ was obtained by smoothing a triangular wave of zero mean using a Hanning window (half width ¼ 6:0 s). For the simulation of the Chi data, SðtÞ was obtained by smoothing a square wave of zero mean, using a Hanning window (half width ¼ 3:0 s). For the simulation of pre-meditation and control data, we used a sine wave of zero mean. The period of SðtÞ was adjusted to each of the data sets to be modeled (Table 2). In the case of the triangular wave used to simulate the Kundalini Yoga data, the percentage of the period occupied by the rising branch of the wave was systematically changed in steps of 10%, to finally set it at 80%. After smoothing the triangular wave, this percentage changed to 60%. The amplitude of the three periodic waves used for SðtÞ was initially set to 1.0. In this way, the final amplitude of the signal injected into the model could be easily reached by multiplying SðtÞ by the value of amplitude desired in each simulation. The amplitudes of SðtÞ and NðtÞ (Table 2) were set, in successive tests, to reproduce (i) the shape of the interval plot and embedding of the beat-to-beat interval data, and (ii) the first autocorrelation coefficient of the beat-to-beat interval data equal or less than zero. Segments of the time dependent modulation SðtÞ þ NðtÞ used for the simulation of meditation are depicted in Fig. 7. SðtÞ was generated using a code in C and then fed into the code for the pacemaker. The model described above simulates the membrane potential vðtÞ of the pacemaker that triggers the heartbeat. The timing of the spikes was determined with an algorithm that detects the crossing of a threshold with positive derivative. The detecting level (threshold) was set at 20:0 mV. The absence of errors in the spike detection was checked graphically. Table 2 Parameter sets used in the numerical simulation of Kundalini Yoga data (row 2), Chi data (row 4) Case

SðtÞ period (s)

SðtÞ amplitude

NðtÞ amplitude

SðtÞ=NðtÞ

A

Pre-Yoga Yoga Pre-Chi Chi

3.3 15.0 4.0 18.2

0.006 0.03 0.008 0.028

0.01 0.018 0.018 0.018

0.6 1.67 0.44 1.56

0.1 0.01 0.09 0.06

The parameters used for the same individuals before meditation are given in rows 1 and 3, respectively.

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

449

3. Results In this section, we describe the experimental data, and make comparison with the model simulations. 3.1. Experimental data In Fig. 2 we display examples of the Kundalini Yoga and Chi data. As reported by Peng et al. [18], extremely prominent heart rate oscillations correlated with breathing appear during both forms of meditation. These oscillations have significantly larger amplitude and lower frequency during meditation than before meditation or in healthy control subjects [18]. The beat-to-beat intervals recorded before and during Kundalini Yoga meditation are shown in Fig. 2A and B, while the beat-to-beat intervals recorded before and during Chi meditation are shown in Fig. 2E and F. We see clearly that the beat-tobeat interval exhibits low-frequency oscillations with very large amplitude during both forms of meditation. The recordings performed in the same subjects before meditation do not display such large oscillations, although the RM can be recognized as low-amplitude oscillations of faster rate (Fig. 2A and E). In the 4 cases of Kundalini Yoga meditation analyzed, we observed (in addition to the oscillations) a pronounced decrease of the mean beat-to-beat interval (e.g., compare Fig. 2A with B). The comparison of the mean interval values before and during meditation with the U-test gives significant differences with a level of risk lower than 1%. The decrease of the mean beat-to-beat interval was small or did not occur in the 8 cases of Chi meditation (e.g., compare Fig. 2E with F) (there are no significant differences in the U-test). The oscillation of the heart beat interval duration is noteworthy in the autocorrelation function of the beat-to-beat series. The full lines in Fig. 3 show the autocorrelation of the cases depicted in Fig. 2. During meditation, the autocorrelation presents high amplitude slowly damped oscillations (full lines in Fig. 3B and D) due to the considerable regularity of the interval oscillations. Before meditation the RM can be recognized in the autocorrelation as low-amplitude oscillations of faster rate (full lines in Fig. 3A and C). The amplitude of the autocorrelation oscillations is larger in Yoga meditation than in Chi, indicating that the beat-to-beat interval oscillations are more regular in these last data sets (e.g., compare full lines in Fig. 3B and D). The larger regularity of the beat-to-beat interval oscillations during Yoga meditation holds in general for all data sets: the amplitude of the first autocorrelation oscillation is significantly larger for Yoga (U-test at 5% of risk). The mean values SD are 0:85 0:22 for Yoga and 0:33 0:24 for Chi. The amplitude of the beat-to-beat interval oscillations, estimated through the SD of the complete beat-to-beat interval series, is not significantly different in both types of meditation. The shape of the embedding changes significantly during meditation in the cases depicted in Fig. 2. The points lie in a roughly circular shape before meditation (Fig. 2C and G). During Kundalini Yoga meditation the points of the embedding lie

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

450

1.2

Beat-to-beat interval (s)

Beat-to-beat interval (s)

1.2

1

0.8

0.6

0.4

0

40

(A)

80

120

Beat

0.4

Beat-to-beat interval [i+7] (s)

Beat-to-beat interval [i+5] (s)

0.8

0.6

0.4 0.4

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

80

120

1

0.8

0.6

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

1.2

Beat-to-beat interval (s)

Beat-to-beat interval (s)

1

0.8

0.6

0

40

80

120

Beat

0.8

0.6

0.4

0

40

80

120

Beat 1.2

1

0.8

0.6

0.4 0.4

1

(F)

Beat-to-beat interval [i+6] (s)

Beat-to-beat interval [i+6] (s)

40 Beat

0.4 0.4

1.2

(G)

0

(D)

1.2

0.4

0.6

1.2

1

(E)

0.8

(B)

1.2

(C)

1

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

1

0.8

0.6

0.4 0.4

(H)

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

451

1

1

0.5

0.5

Autocorrelation

Autocorrelation

0.75

0.25 0 -0.25

0

-0.5

-0.5 -0.75

0

1

2

(A)

3

4

5

-1

6

0

5

10

(B)

Interval order

15

20

25

30

20

25

30

Interval order 1

1

0.5

0.5

Autocorrelation

Autocorrelation

0.75

0.25 0 -0.25

0

-0.5

-0.5 -0.75

(C)

0

1

2

3

Interval order

4

5

1-

6

(D)

0

5

10

15

Interval order

Fig. 3. (A) Full lines: autocorrelation functions of the beat-to-beat interval series shown in Fig. 2. Dashed lines: 99% confidence intervals of the autocorrelation functions corresponding to the simulations depicted in each panel (N ¼ 50 model realizations). The beat-to-beat intervals of single realizations of these four simulations are shown in Fig. 5. (A) Before Yoga meditation. (B) During Yoga meditation. (C) Before Chi meditation. (D) During Chi meditation.

on a triangular shape (Fig. 2D), and during Chi meditation the shape is quadrangular (Fig. 2H). The shape of the embedding is related to the shape of the beat-to-beat interval plot. The triangular embedding is associated with the beat-to-beat interval oscillation having sharp peaks (Fig. 2B), while the quadrangular embedding is related to beat-to-beat interval oscillation having peaks and valleys of similar width (Fig. 2F). A remarkable property of the beat-to-beat interval data during meditation is the presence of noise induced indentations both in the peaks and valleys of the interval Fig. 2. Beat-to-beat interval (s) vs. beat, (A) before and (B) during Kundalini Yoga meditation. (C) Embedding of the time series shown in the panel A (300 intervals with delay ¼ 5). (D) Embedding of the time series shown in the panel B (300 intervals with delay ¼ 7). Beat-to-beat interval (s) vs. beat, (E) before and (F) during Chi meditation. (G) Embedding of the time series shown in the panel E (300 intervals with delay ¼ 6). (H) Embedding of the time series shown in the panel F (300 intervals with delay ¼ 6).

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

452

oscillations (Fig. 2B and F). In other words, the edges of the embeddings (triangular or quadrangular) have similar width, i.e. they are similarly noisy (Fig. 2D and H). This indicates that the SA node operates in a regime in which the relationship between the input modulations from the autonomous nervous system and the output beat-to-beat interval series is approximately linear, i.e. the saturation level described in Goldberger et al. [32] is not reached. This property has to be taken into account by the model in order to properly reproduce the embedding shape of the data. In the next subsection, we shall assess through numerical simulations the shape of the input signals needed to reproduce the results presented above. The simulations are based in the cases shown in the Fig. 2, because these segments have a very regular waveshape in the beat-to-beat interval plot. Regarding the recordings from the other individuals, in 3 out of 4 cases of Kundalini Yoga meditation analyzed, the embedding and the beat-to-beat interval plots are similar to the case depicted in the Fig. 2, for example Fig. 4A and B. In the remaining case the embedding is more irregular although it also evokes a triangular shape (not shown). In the case of Chi meditation we found larger differences between the different subjects. In 3 out of 8

1.2 Beat-to-beat interval [i+8] (s)

Beat-to-beat interval [i+9] (s

1.2

1

0.8

0.6

0.4 0.4

(A)

0.6 0.8 1 Beat-to-beat interval [i] (s)

0.6

(B)

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

0.6

1.2

1.2 Beat-to-beat interval [i+4] (s)

Beat-to-beat interval [i+6] (s)

(C)

0.8

0.4 0.4

1.2

1.2

1

0.8

0.6

0.4 0.4

1

0.6

0.8

1

Beat-to-beat interval [i] (s)

1

0.8

0.6

0.4 0.4

1.2

(D)

0.8

1

Beat-to-beat interval [i] (s)

Fig. 4. Embedding of the beat-to-beat intervals of different individuals during Yoga meditation, (A) 1121 intervals, (B) 906 intervals. Embedding of the beat-to-beat intervals of different individuals during Chi meditation, (C) 4027 intervals, (D) 3759 intervals.

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

453

cases the embedding shape and the beat-to-beat interval plot are similar to the case shown in the Fig. 2H and F, for example Fig. 4C. In the remaining 5 subjects, the embedding is oval or has irregular shape (for example Fig. 4D), and the oscillation of the beat-to-beat interval plot is less regular than in Fig. 2F. In some of these cases the waveshape of the oscillation resembles that before meditation although the amplitude is higher and the frequency lower. In general the pre-meditation recordings present roughly circular shapes as shown in Fig. 2C and G, but in many cases they are more irregular as also happens with the control subjects during sleeping hours, and with the subjects during metronomic breathing (not shown). The larger variability within different subjects observed in the case of Chi meditation may be related to the shorter training period of the practitioners used for recording the data. The reason being that a larger training period enables the practitioner to obtain a more regular breathing rhythm during meditation. In Fig. 4 we show examples of embeddings calculated from the complete recording time in different individuals during meditation. 3.2. Simulation of the beat-to-beat interval series Using the model for the RM of the SA node (Eqs. (1)–(17)), we generated numerically the beat-to-beat interval series during and before both types of meditation. In successive tests, we selected the appropriate parameter sets for obtaining good agreement with the experimental data. These parameters are displayed in the Table 2. From this table we can draw the following conclusions. The period and the amplitude of SðtÞ are larger for meditation data. The ratio SðtÞ=NðtÞ (given in the 5th column of the Table 2) is larger during meditation. The parameter A, set so as to fit the mean rate of each recording, resulted larger during meditation, and resulted larger for Yoga than for Chi meditation. The parameters obtained for the control subjects are similar to those for the pre-meditation situation (not shown). 3.2.1. Simulation of Kundalini Yoga data The Kundalini Yoga data could be reproduced using a signal SðtÞ obtained by smoothing a triangular wave. The rising branch of this signal lasted 60% of the period (Fig. 7A). In a simulation using the parameter set given in the second row of the Table 2, we obtained a beat-to-beat interval oscillation with sharp peaks (Fig. 5B), and a triangular embedding (Fig. 5D), just like the experimental data (Fig. 2B and D). The pre-meditation data could be reproduced using a sinusoidal SðtÞ, the points of the embedding falling roughly inside a circle (Fig. 5C). The waveform of the periodic signal was not crucial in this case, because the amplitude of SðtÞ for pre-meditation is much smaller than in the case of meditation (see parameters in the Table 2). In these conditions, SðtÞ cannot be distinguished from the ‘‘background noise’’ represented by NðtÞ. The key factor to obtain a polygonal embedding in the simulations is the amplitude ratio of the periodic and aperiodic modulations, SðtÞ=NðtÞ (5th column of the Table 2). The embedding of the simulation before meditation becomes triangular

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

454

1.2

Beat-to-beat interval (s)

Beat-to-beat interval (s)

1.2

1

0.8

0.6

0.4

0

40

(A)

80

0.8

0.6

0.4

120

Beat

1

1

0.8

0.6

0.4 0.4

(C)

0.6 0.8 1 Beat-to-beat interval [i] (s)

100

120

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

1.2

Beat-to-beat interval (s)

Beat-to-beat interval (s)

80

0.6

0.4 0.4

0.8

0.6

0

40

80

0.8

0.6

0

(F)

Beat

40

80

120

Beat

Beat-to-beat-interval [i+6] (s)

1.2

1

0.8

0.6

0.4 0.4

1

0.4

120

1.2 Beat-to-beat interval [i+6] (s)

60 Beat

0.8

(D)

1

(G)

40

1

1.2

1.2

0.4

20

1.2 Beat-to-beat interval [i+7] (s)

Beat-to-beat interval [i+5] (s)

1.2

(E)

0

(B)

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

1

0.8

0.6

0.4 0.4

(H)

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

1

0.8

0.6

20% 0.4 0.4

(A)

1.2 Beat-to-beat interval [i+7] (s)

Beat-to-beat interval [i+7] (s)

1.2

455

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

(B)

1

0.8

0.6

50% 0.4 0.4

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

Beat-to-beat interval [i+7] (s)

1.2

1

0.8

0.6

80% 0.4 0.4

(C)

0.6 0.8 1 Beat-to-beat interval [i] (s)

1.2

Fig. 6. Embedding of the beat-to-beat interval series simulated using periodic signals with different percentages of the period occupied by the rising branch of the wave. The general parameters were as in the second row of Table 2 with the exceptions of SðtÞ amplitude and period, here set to 0.033 and 16.7 s, respectively. The rising branch of SðtÞ occupied 20% of the period in (A), 50% of the period in (B), and 80% of the period in (C). These percentages are printed in the respective panels.

if the amplitude of the periodic signal SðtÞ is increased, the other parameters being kept fixed. Similarly, the polygonal embedding obtained in the simulations of meditation progressively looses definition as the ratio SðtÞ=NðtÞ decreases. The spatial orientation of the embedding varied with the percentage of the period occupied by the rising branch of SðtÞ. This feature can be better appreciated in simulations in which the triangular wave SðtÞ was not smoothed with the Hanning window (Fig. 6). Using rising branches ranging from 10% to 30% we obtained a

Fig. 5. Numerical simulation of the beat-to-beat interval (s) vs. beat, (A) before and (B) during Kundalini Yoga meditation. (C) Embedding of the time series shown in the panel A (240 intervals with delay ¼ 5). (D) Embedding of the time series shown in the panel B (240 intervals with delay ¼ 7). Numerical simulation of the beat-to-beat interval (s) vs. beat, (E) before and (F) during Chi meditation. (G) Embedding of the time series shown in the panel E (240 intervals with delay ¼ 6). (H) Embedding of the time series shown in the panel F (240 intervals with delay ¼ 6).

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

456

triangular embedding (Fig. 6A), but its orientation is different from the Yoga data (Fig. 2D). For rising branches of 50% the embedding tends to a quadrangular shape (Fig. 6B) because the left edge of the triangle shown in the Fig. 6A bends forming a new angle. The orientation of this quadrangular shape (Fig. 6B) is different from the embedding of the Chi data (Fig. 2H). For rising branches ranging from 70% to 90% of the period, the embedding becomes triangular again, and this time the orientation is as observed in the Yoga data (Fig. 6C). Thus, the orientation of the embedding is related to the percentage of the period occupied by the rising branch of SðtÞ. We chose the percentage that better reproduced the embedding shape of the data. This percentage was 80% without smoothing the triangular wave (Fig. 6C), and changed to 60% after the smoothing (Fig. 5D). The shape of the simulated embedding matched the data embedding only when we used the smoothed SðtÞ. To check the criterion of the embedding orientation, we counted the number of beats in the increasing and decreasing branches of the beat-to-beat interval oscillations of the experimental and numerical data. In both cases, the percentage of the decreasing branch of the interval oscillation (that corresponds to the rising branch of SðtÞ) has a good agreement with the value selected according to the embedding orientation. The results described above suggest that the respiratory signal of the Kundalini Yoga meditation may look like the noisy periodic wave shown in Fig. 7A, in the sense that the wave has sharp peaks, and that the inhalation lasts longer than the exhalation.

S(t)+N(t)

3.2.2. Simulation of Chi data The Chi data could be reproduced using a signal SðtÞ obtained by smoothing a square wave (Fig. 7B). In the simulation using the parameter set given in the fourth

0.1

0.1

0.05

0.05

0

0

-0.05

-0.05

-0.1

(A)

0

20

40

60

t (s)

80

100

-0.1

(B)

20

40

60

80

100

t (s)

Fig. 7. Single realization of the time dependent input modulation SðtÞ þ NðtÞ obtained with the parameter set provided in the Table 2 used for the numerical simulation of Kundalini Yoga (A) and Chi (B) meditation data. The ordinate axis shows the amplitude of the modulation on iNaK (dimensionless). The same modulations were used for iK with an amplitude 5.0 times larger.

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

457

row of the Table 2, we obtained a beat-to-beat interval oscillation with wide peaks and valleys (Fig. 5F), and a quadrangular embedding (Fig. 5H), just like the experimental data (Fig. 2F and H). As the embedding orientation of the model coincided with that of the experimental data, we concluded that in this case the rising and decreasing branches of the signal SðtÞ should have the same duration. To confirm this, we checked that both branches of the beat-to-beat interval oscillation had similar number of intervals, both in the experimental and numerical data. These results suggest that the respiratory signal of Chi meditation may look like the noisy periodic wave shown in Fig. 7B, in the sense that it has peaks and valleys of approximatively same width. The pre-meditation data could be reproduced using a sine, as for the pre-Yoga data. The points of the embedding fall in a roughly circular shape (Fig. 5G). The data sets recorded during Chi meditation in which the embedding shape was not square but oval (e.g., Fig. 4D) could be reproduced using a sine of larger amplitude and period than for the case before meditation (not shown). 3.2.3. Statistical validation of the numerical simulations The simulations presented in Fig. 5 correspond to single realizations of the model. As each realization is performed with a different seed of noise (Eq. (1)), we calculated the averaged autocorrelation function of 50 model realizations with its 99% confidence intervals, in order to check if the autocorrelation of the beat-to-beat interval data is contained within these limits. In Fig. 3 the first oscillation of the autocorrelation functions estimated from the four data segments that provided the parameter sets for each simulation (i.e. before and during Yoga meditation, and before and during Chi meditation) are plotted with full line. Note that in the four cases the full lines are within the dashed lines that represent the confidence intervals of the simulations. The only exception is the case of Chi meditation for autocorrelation coefficients of order larger than 16. In this case the autocorrelation oscillation of the simulation is damped more slowly than in the data.

4. Discussion The present work is based on the paper by Peng et al. [18] that first reported the occurrence of high-amplitude heart rate oscillations during Chi and Kundalini Yoga meditation practices. These authors applied spectral analysis and a technique based on Hilbert transform to quantify the heart rate oscillations. They found that the amplitude of the oscillations and its Fourier and Hilbert powers are significantly greater during meditation than before meditation or in control subjects. Measuring fluctuations in the mean cardiac electrical axis that accompany breathing [18,33,34] derived respiratory signals from body surface ECGs and showed that the heart rate oscillations are correlated with respiration at the same rate. Our contribution starts with the analysis of data segments using a bidimensional embedding [30], and continues with the development of a biophysical model to

ARTICLE IN PRESS 458

A. Capurro et al. / Physica A 355 (2005) 439–460

reproduce the RM of the beat-to-beat interval data. In this model, we assumed that the respiratory input from the autonomous nervous system to the SA node mediates the fluctuations between heart beat intervals in short time scales [8], as the duration of the selected data segments in which our simulations are based (200–300 s). The origin of the rhythmic activity of the colinergic input to the SA node is a very complex issue, but it seems clear that the respiratory movements cause parallel variations in the vagal tone that affect the discharging rate of the SA node [2]. This is supported by early works showing that the RM is suppressed by cervical vagotomy [35]. In addition, the high-frequency power band of the HRV (0.15–0.4 Hz), which contains the RM, is suppressed by parasympathetical blocking with atropine infusion [36]. Thus, in our model we represent the high-frequency power band of the HRV with a periodic wave with added correlated noise. The lowfrequency band (0.04–0.15 Hz), due to baroreceptor and thermoregulatory influences, is represented only generically by the added correlated noise. As the power spectrum of NðtÞ covers both the high- and low-frequency bands of the HRV, this signal is addressed to aperiodic influences affecting the RM and the lowfrequency power band. We ignored other components of heart rate regulation, contained in the very low power (0.0033–0.04 Hz) and ultra low power (less than 0.0033 Hz) frequency bands of the HRV, which are excluded due to the use of short lasting data segments without trends. The very low-frequency band is due to the renin-angiotensin and peripheral vasomotor systems, while the ultra low-frequency band is due to neuroendocrine and circadian rhythms [8]. In the following we discuss some specific points of our results. 4.1. The RM during meditation The first claim that we can make from our analysis of the data confirms the results shown in Peng et al. [18]: during both types of meditation the RM on HRV becomes very large, and predominates over aperiodic components. Before meditation, although noteworthy in the autocorrelation function, the RM is hidden in a background of aperiodic variability. The large difference in the value of the parameter SðtÞ and in the ratio SðtÞ=NðtÞ used in the simulation of meditation and pre-meditation data (Table 2) is a consequence of the very high amplitude of the oscillations that appear during meditation. It is important to stress the similarity of the parameters for both meditation techniques. This is in accordance with the suggestion that meditation protocols from different cultures may share certain types of autonomic exercises mediated in part by specialized breathing maneuvers [18]. 4.2. Sustained rate changes and resting vagal tone The parameter A allowed us to reproduce in the model the sustained changes observed in the heart rate. This is remarkable because, using simpler oscillators we could reproduce qualitatively the waveshape of the meditation data [30], but not the sustained changes in rate. With the present model, we are able to reproduce

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

459

quantitatively the heart rate as well as the amplitude and waveshape of the variability found in different data sets. The parameter A represents the existence of a resting vagal tone that continuously slows down the frequency of the SA node, more precisely the vagal tone is represented by A. The pronounced decrease in the mean value of the beat-to-beat interval observed during Kundalini Yoga meditation may be related to a decrease in the tonic vagal action, i.e. this meditation type has the larger value of A (Table 2). The concomitant large increase of the RM amplitude suggests that the oscillatory component of the vagal action (i.e. the mirror image of SðtÞ þ NðtÞ) may be greatly enhanced. Thus, there may exist a dissociation in the regulation of the oscillatory and tonic components of the vagal action. 4.3. Differences in the time course of respiration The use of the embedding allowed us to focus on certain differences between the time course of respiration followed in both forms of meditation. The shape of the embedding is triangular for Kundalini Yoga data, and quadrangular for Chi data. These shapes are related to the beat-to-beat interval plots, but in such plots the differences are not easily detected due to the noisy components of the variability. In the embedding, although the trajectories vary, they are limited to geometric shapes easy to recognize. These shapes can be considered as hallmarks of each meditation type. If we knew the relationship between SðtÞ þ NðtÞ and the respiratory movements, the time course of respiration could be inferred from the variability of the heart beatto-beat interval series. We assumed that this relationship is linear, because this seems to be the case for frequencies in the range of breathing [7]. With this condition, the respiration during Kundalini Yoga and Chi meditations should have the waveshapes shown in Fig. 7. Respiratory signals may be derived from body surface ECGs by measuring fluctuations in the mean cardiac electrical axis that accompany breathing [33,34]. These ECG derived estimations can be produced using the same ECG voltage files from which the heart beat interval series where obtained. There exists the possibility of contrasting the respiratory signals shown in Fig. 7 with these ECG derived estimations of the respiratory movements. These comparisons will be carried out in a future study, when these ECG voltage file are made available in physionet [19].

Acknowledgements AC acknowledges financial support from FAPESP (Fundac- a˜o de Amparo a` Pesquisa do Estado de Sa˜o Paulo), process number 00/14357-8. CPM acknowledges partial financial support from CNPq and FAPESP. AC is researcher of PEDECIBA (Uruguay) and thanks to the EU Advanced Course in Computational Neuroscience 2001 (Trieste/Italy).

ARTICLE IN PRESS A. Capurro et al. / Physica A 355 (2005) 439–460

460

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

[20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]

D. DiFrancesco, Ann. Rev. Physiol. 55 (1993) 455. R.M. Berne, M.N. Levy (Eds.), Physiology, Mosby-Year Book, St. Louis, MO, 1993 (Chapter 25). D. DiFrancesco, C. Tromba, J. Physiol. 405 (1988) 477. D. DiFrancesco, C. Tromba, J. Physiol. 405 (1988) 493. D. DiFrancesco, P. Tortora, Nature 351 (1991) 145. D. DiFrancesco, P. Ducouret, R.B. Robinson, Science 243 (1989) 669. J.P. Saul, R.D. Berger, M.H. Chen, R.J. Cohen, Am. J. Physiol. Heart Circ. Physiol. 256 (1989) H153. P.K. Stein, R.E. Kleiger, Annu. Rev. Med. 50 (1999) 249. N.D. Giardino, R.W. Glenny, S. Borson, L. Chan, Am. J. Physiol. Heart Circ. Physiol. 284 (2003) H1585. A.L. Goldberger, L.A.N. Amaral, J.M. Hausdorff, P.Ch. Ivanov, C.K. Peng, H.E. Stanley, Proc. Natl. Acad. Sci. USA 99 (2002) 2466. M. Costa, A.L. Goldberger, C.K. Peng, Phys. Rev. Lett. 89 (2002) 068102. A. Babloyantz, A. Destexhe, Biol. Cybern. 58 (1988) 203. C.S. Poon, C.K. Merrill, Nature 389 (1997) 492. Y. Ichimaru, K.P. Clark, J. Ringler, W.J. Weiss, Comput. Cardiol. 17 (1990) 657. M. Sakakibara, S. Takeuchi, J. Hayano, Psychophysiology 31 (1994) 223. L.A. Lipsitz, F. Hashimoto, L.P. Lubowsky, J.E. Mietus, G.B. Moody, O. Appenzeller, A.L. Goldberger, Br. Heart J. 74 (1995) 390. P.F. Migeotte, G.K. Prisk, M. Paiva, Am. J. Physiol. Heart Circ. Physiol. 284 (2003) H1995. C.K. Peng, J.E. Mietus, Yanhui Liu, K. Gurucharan, P.S. Douglas, H. Benson, A.L. Goldberger, Int. J. Cardiol. 70 (1999) 101. A.L. Goldberger, L.A.N. Amaral, L. Glass, J.M. Hausdorff, P.Ch. Ivanov, R.G. Mark, J.E. Mietus, G.B. Moody, C.K. Peng, H.E. Stanley, PhysioBank, physioToolkit and physionet: Components of a new research resource for complex physiologic signals, Circulation 101 (2000) 23 http:// circ.ahajournals.org/cgi/content/full/101/23/e215. M. Brennan, M. Palaniswami, P. Kamen, Am. J. Physiol. Heart Circ. Physiol. 283 (2002) H1873. H.F. Brown, J. Kimura, D. Noble, S.J. Noble, A. Taupignon, Proc. R. Soc. London Ser. B 222 (1984) 329. S.S. Demir, J.W. Clark, C.R. Murphey, W.R. Giles, J. Physiol. Cell physiol. 266 (1994) C832. L.P. Endresen, K. Hall, J.S. Hoye, J. Myrheim, Eur. Biophys. J. 29 (2000) 90. H. Zhang, A.V. Holden, I. Kodama, H. Honjo, M. Lei, T. Vagues, M.R. Boyett, Am. J. Physiol. Heart Circ. Physiol. 279 (2000) H397. Y. Kurata, I. Hisatome, S. Imanishi, T. Shibamoto, Am. J. Physiol. Heart Circ. Physiol. 283 (2002) H2074. T. Shibasaki, J. Physiol. (London) 387 (1987) 227. J. Guo, K. Ono, A. Noma, J. Physiol. (London) 483 (1995) 1. A. Noma, Jpn Heart J. 37 (1996) 673. A. Zaza, M. Micheletti, A. Brioschi, M. Rocchetti, J. Physiol. (London) 505 (1997) 677. A. Capurro, L. Diambra, C.P. Malta, Physica A 327 (2003) 168. A. Capurro, K. Pakdaman, T. Nomura, S. Sato, Phys. Rev. E 58 (1998) 4820. J.J. Goldberger, S. Challapalli, R. Tung, M.A. Parker, A.H. Kadish, Circulation 103 (2001) 1977. G.B. Moody, R.G. Mark, A. Zoccola, S. Mantero, Comput. Cardiol. 12 (1985) 113. G.B. Moody, R.G. Mark, M.A. Bump, J.S. Weinstein, A.D. Berman, J.E. Mietus, A.L. Goldberger, Comput. Cardiol. 13 (1986) 507. M. De Burgh Daly, in: R.M. Berne, N. Sperelakis (Eds.), Handbook of Physiology, 1979, p. 529 (Chapter 16). R. Pomeranz, J. Macaulay, M.A. Caudill, I. Kutz, D. Adam, D. Gordon, K.M. Kilborn, A.C. Barger, D.C. Shannon, et al., Am. J. Physiol. Heart Circ. Physiol. 284 (1985) H151.