Intracardiac Origin of Heart Rate Variability, Pacemaker Funny Current and their Possible Association with Critical Illness

Heart rate variability (HRV) is an indirect estimator of autonomic modulation of heart rate and is considered a risk marker in critical illness, particularly in heart failure and severe sepsis. A reduced HRV has been found in critically ill patients and has been associated with neuro-autonomic uncoupling or decreased baroreflex sensitivity. However, results from human and animal experimental studies indicate that intracardiac mechanisms might also be responsible for interbeat fluctuations. These studies have demonstrated that different membrane channel proteins and especially the so-called ‘funny’ current (If), an hyperpolarization-activated, inward current that drives diastolic depolarization resulting in spontaneous activity in cardiac pacemaker cells, are altered during critical illness. Furthermore, membrane channels kinetics seem to have significant impact upon HRV, whose early decrease might reflect a cellular metabolic stress. In this review article we present research findings regarding intracardiac origin of HRV, at the cellular level and in both isolated sinoatrial node and whole ex vivo heart preparations. In addition, we will review results from various experimental studies that support the interrelation between If and HRV during endotoxemia. We suggest that reduced HRV during sepsis could also be associated with altered pacemaker cell membrane properties, due to ionic current remodeling.


INTRODUCTION
Healthy state exhibits some degree of stochastic variability in physiologic variables, such as heart rate (i.e. heart rate variability-HRV). This variability is a measure of complexity that accompanies healthy systems and has been suggested as responsible, for their greater adaptability and functionality related to pathologic systems [1]. Loss of such variability means loss of complexity that accompanies cardiovascular disease and critical illness [2,3], while it is associated with increased mortality rate after acute myocardial infarction [4], sepsis and multiple organ dysfunction syndrome (MODS) [3,5].
Studying physiological signals of critically ill patients, such as heart rate, can easily identify 'hidden' information concerning inherent dynamics and overall variability within time series [2]. Recognition that physiologic time series contain hidden information related to an extraordinary complexity that characterizes physiologic systems defies traditional mechanistic approaches based on conventional biostatistical methodologies and has fueled growing interest in applying techniques from statistical physics for the study of living organisms [1,2]. This effort has boosted research on heart rate and blood pressure dynamics through standardization of different signal processing techniques, frequency and duration of measurements and signal quality assessment, and has stimulated the development of more accurate diagnostic *Address correspondence to this author at the Vasilios E Papaioannou: Democritus University of Thrace, Alexandroupolis Medical School, Intensive Care Unit, Dragana 68100, Greece; Tel: 0030-6942551414; E-mail: vapapa@med.duth.gr and prognostic indices in cardiovascular diseases [6]. Furthermore, a number of international databases of heart rate signals have been developed with free access from different investigators, such as the Web Site Physionet (www.physionet.org) [6]. In conclusion, the combination of structural indices such as the left ventricular ejection fraction (LVEF) with autonomic function indices derived from heart rate variability analysis has been recently proposed as the stateof-the-art method for risk assessment among patients with acute myocardial infarction or severe congestive heart failure [7].
In this review article we will describe the basic methods for assessment of HRV, discuss potential mechanisms responsible for its generation at the intracardiac level and summarize the most recent and significant studies, from both basic and clinical research, which investigate the potential effects of cardiovascular diseases and sepsis upon HRV alterations. Finally, we will present research data concerning the association between critical illness, HRV and the 'funny' current (I f ) that is responsible, among others, for the generation of action potentials (APs) by pacemaker cells in the sinoatrial node (SAN).

MEASUREMENT OF HRV
Heart rate variability describes variations in both instantaneous heart rate and RR intervals. HRV is considered an indirect measure of autonomic regulation of cardiac activity, and can reflect the coupling between the autonomic nervous system (ANS) and the sinoatrial node [3]. In 1996, the Task Force of the European Society of Cardiology and the Northern American Society of Pacing and Electrophysiology pub-lished guidelines regarding standardization of nomenclature, specification of methods of measurement, definition of physiological and pathophysiological correlates, description of clinical applications and identification of different areas for future research [8].
The RR variations may be evaluated by two methods derived from: 1) time domain, and 2) frequency domain.

Time Domain Methods
Time domain methods determine heart rate or RR intervals in continuous electrocardiographic (ECG) recordings. Each QRS complex is detected and the normal-to-normal intervals (all intervals between adjacent QRS complexes) are calculated. Other time domain variables include the mean normal-to-normal interval, the mean heart rate or the difference between the longest and the shortest interval as well. There are also more complex statistical methods being used, particularly from heart rate signals being recorded for more than 24 hours. The simplest from these metrics is the standard deviation of the normal-to-normal intervals (SDNN), which is the square root of the variance. However, it should be emphasized that the SDNN variable becomes less accurate by shorter the monitoring periods since it depends on the duration of RR interval series and its duration must be standardized (e.g. 5 min or 24 hours) [8]. The most commonly used time domain methods are the square root of the mean squared differences of successive intervals (RMSSD), the number of interval differences of successive intervals greater than 50 ms (NN50) and the proportion derived from dividing NN50 by the total NN intervals (pNN50) [8].

Frequency Domain Methods
Akselrod et al. [9] introduced in 1981 power spectrum analysis of heart rate fluctuations in order to quantify beatto-beat cardiovascular control. Power spectrum density (PSD) analysis provides the basic information of how power (variance, msec 2 /Hz) distributes as a function of frequency. Spectral analysis of heart rate signals provides their power spectrum densityand displays in a plot the relative contribution (amplitude) of each frequency, after application of a Fast Fourier transformation (FFT) to the raw signal. This plot includes at least three frequency peaks. Fast frequency periodicities (high frequency, HF), in the range 0.15-0.4 Hz, are largely due to the influence of the respiratory phase on vagal tone. Low-frequency periodicities (LF), in the range of 0.04-0.15 Hz, are produced by baroreflex feedback loops, affected by both sympathetic and parasympathetic modulation of the heart. Very low frequency periodicities (VLF), i.e. less than 0.04 Hz, have been variously ascribed to modulation by chemoreception, thermoregulation and the influence of vasomotor activity, which is related, between others, to the renin-angiotensin-aldosterone system (RAS) [8][9][10]. The area under the power spectral curve in a particular frequency band (power) is considered to be a measure of heart rate variability at that frequency. The ratio LF/HF reflects sympathovagal balance whereas normalized units (nu) of both LF and HF (LF/total power and HF/total power, respectively) indicate heart rate variability in specific bands irrespectively of total variability of the whole signal [8]. In a double logarithmic plot of power versus frequency, their relation follows a straight line with a slope defined as β. This relation is known as the power law, whereas in normal subjects, β slope or exponent is close to -1 [8,11].

Extracardiac Origin of HRV
The LF component of HRV is probably the most contentious aspect with respect to cardiovascular variability. There are two opposing theories in the literature proposing different potential origins: 1) the central oscillator theory, and 2) the baroreflex feedback loop theory [12,13]. According to the first theory, it is believed that LF oscillations reflect sympathetic tone and are generated by the brain stem circuits. In cats, Montano et al. [12] analyzed the discharges of single sympathetic neurons located in the rostral ventrolateral medulla (RVLM) and caudal ventrolateral medulla (CVLM). They observed activity at 0.12 Hz, which was positively correlated with heart rate and blood pressure variability. As the above oscillations remained after sino-aortic and vagal resection, it was assumed that the central nervous system is able to generate such oscillations.
The second, more accepted theory is baroreflex feedback loop model [13], where a change in blood pressure is sensed by arterial baroreceptors, resulting in heart rate adjustment through the central nervous system and via both the fast vagal action and the slower sympathetic action. At the same time, baroreceptors induce a slow sympathetic withdrawal from the vessels. The delay in the sympathetic branch of the baroreflex in turn determines a new oscillation, which is sensed by the baroreflex and induces a new oscillation in heart rate. It has been also proposed that the LF oscillation arises from the interaction of slow sympathetic and fast vagal responses, where baroreflex buffering of the slow respiratory induced blood pressure oscillations results in resonant low frequency oscillations, due to the delay in the slow conducting sympathetic loop of the baroreflex [14].

Power Law and HRV
In terms of power spectrum density, the major component of HRV occurs at frequencies below 0.04 Hz, where its power spectrum exhibits a power law behavior.Fluctuations of a variable can be characterized by its probability density distribution. A way of estimating its characteristics is the construction of a histogram after normalization, so that the area under it will be equal to one. Often, this distribution N(x) of a variable x follows the power law form: N(x) = k*xd meaning that the relative frequency of a value x is proportional to x raised to the power of -d, whereas a constant multiplicative factor k (usually different from unity) must be incorporated in the equation. If we plot the logarithms of this relationship we have a linear equation: log (N) = log (k)d*log(x), whereas d is the negative slope of a straight line fit to N. This slope is the β slope or exponent. (Fig. 1A & B) [15,16].
Power law distribution behaves differently than Gaussian distributions. In power law distributions, tails are very long (long-tail distribution), representing the relative frequency of occurrence of large events. This means that the probability of large or rare events is much higher compared with a Gaussian. Power laws describe dynamics that have a similar pattern of change at different scales and they are called 'scale A B Fig. (1). Fast Fourier transformation (FFT) of a heart rate time series (A) and its power law distribution (B). The spectrum of the heart rate time series displays three peaks mainly at the lower frequencies (A) whereas the log transformation of the power of the signal gives rise to a more or less linear relation (B). If we apply log transformation of the frequency as well (not shown here), we will have a more smooth linear relation with a slope β. Data have been downloaded from the free-download website physionet (www.physionet.org).
invariant'. On the contrary, Gaussians are characterized by typical values, such as those corresponding to their peaks [16]. Moreover, the power law describes a time series whose pattern of variation is statistically similar regardless of its size. Magnifying or shrinking the scale of the signal reveals the same relationship, a property that has been called 'selfsimilarity' and is a fundamental characteristic of fractals [16,17]. These structures appear similar at different scales of magnification, meaning that in cases of different time series (e.g., heart rate) this property of self-similarity indicates the presence of long-term correlations across multiple temporal scales. Origin of power law behavior in healthy state is still not known, whereas decrease or more negative values of β slope have been described in different situations, such as myocardial infarction [11] and severe sepsis [5].

Intracardiac Origin of HRV
Except for extracardiac mechanisms, recent evidence from in vitro and ex vivo experimental studies suggests that HRV might also depend on factors intrinsic to cardiac tissue, whereas fluctuations in spontaneously beating rate of single cardiac cells , here defined as beating rate variability (BRV), seem to follow a power law behavior [18][19][20][21]. Moreover, different clinical studies in heart transplant recipients have found evidence for heart rate fluctuations originating from the heart itself [22,23]. Bernardi et al. [22] studied intrinsic mechanism regulating HRV in both transplanted and intact heart during exercise. They found that that at peak exercise a non-autonomic mechanism, probably intrinsic to the heart muscle, may determine heart rate fluctuations in synchrony with ventilation, in transplanted as well as in intact hearts. Hrushesky and colleagues [23] quantified respiratory sinus arrhythmia (RSA) and found that individuals with a transplanted heart had resting RSA values similar to healthy subjects.

Cardiac Pacemaking Mechanism
Cardiac pacemaking is a basic physiological function required to match cardiac performance to metabolic demand. In the mammalian heart, it is accomplished by pacemaker cells in the sinoatrial node region, which is located in the right atrium. The SAN consists of specialized cells that exhibit spontaneous electrical activity, i.e., generating of recurring action potentials due to slow diastolic depolarization ( Fig. 2A). Diastolic depolarization is due to a small net inward current across the cell membrane ( (Fig. 2B), blue solid line), which is the result of the "membrane clock", a complex interaction of multiple time-and voltage-dependent inwardly and outwardly directed ion currents, including the 'funny' current, I f ( (Fig. 2C), red dashed line) [24,25]. I f is called 'funny' current because it is activated upon hyperpolarization and is carried by Na + and K + transport. (Fig. 3A) shows a typical example recorded in a SAN cell isolated from a human heart. Typically, I f becomes larger and activated more rapidly at increasingly negative potentials. (Fig.  3B) shows the average I step , i.e., the current measured at the end of the 2 sec hyperpolarizing voltage clamp step, and I tail , i.e. the currents measured by stepping back to the holding potential of -40 mV, of 3 human SAN cells. Voltagedependence of I f activation is analyzed by plotting normal- ized I tail amplitude against the preceding voltage steps, and is shown in (Fig. 3C). The average half-maximal activation voltage (V 1/2 ) and slope factor of the Boltzmann fit to the data (black dashed line) was -96.9±2.7 and -8.8±0.5 mV, respectively. Molecular characterization has demonstrated that I f is encoded by a family of 4 genes, termed Hyperpolarization-activated Cyclic Nucleotide (HCN) genes. I f observed in human SAN is likely carried by HCN4 channel proteins because they are the dominant expressed isoforms in human SAN [24][25][26]. I f is an important target of heart rate regulation from the autonomic nervous system [24]. In particular, β-adrenergic stimulation through its activation of adenylyl cyclase elevates cAMP, which then binds to the cytoplasmatic tails of I f channel and shifts the voltagedependency of activation to more positive potentials. (Fig.  3C) shows a schematic drawing of such an effect (green dashed line). As a consequence, there is an acceleration of spontaneous electrical activity through an increase of the rate of diastolic depolarization as schematically drawn in (Fig.  3D), green solid line. On the contrary, acetylcholine via muscarinic receptors inhibits adenylyl cyclase and shifts the I f voltage-dependency of activation to more negative poten-tials ( (Fig. 3C), red dashed line), thereby decreasing the diastolic depolarization rate (( Fig. 3D), red solid line) [24][25][26].
In the last decade, spontaneous Ca 2+ releases from the sarcoplasmic reticulum (SR) have been recognized as an additional mechanism for SAN cell pacemaker activity: the ''Ca 2+ clock'' [27]. Lakatta and colleagues [27] showed that spontaneous local subsarcolemmal Ca 2+ released from the sarcoplasmatic recticulum can result in an inward Na + /Ca + exchange current (I NCX ) that contributes to the diastolic depolarization rate and to the initiation of the AP.
The SAN displays a functional and anatomic inhomogeneity because SAN cells have different intrinsic cycle lengths (CLs) or interbeat intervals (IBIs), which correspond to the interval between consecutive AP upstrokes. Moreover, SAN cells differ in terms of responsiveness to ions, drugs, temperature or neuro-hormones [28]. These differences are related with regional anatomic and molecular heterogeneity of single cells within the node [29,30]. For instance, cells from the center of the SAN are small and have a more positive resting membrane potential, a slower diastolic depolarization phase, and slower AP upstroke, compared to larger cells that are found in the periphery. Multiple ion channels contribute to the differences in AP morphology of central and peripheral SAN cells (for review, see [26,29]). In addition, nodal cells in the center of the SAN are poorly electrically coupled due to low expression of connexins, which constitute the gap junctions. As a consequence, AP propagation in the center of the SAN is slow, making this region electrically insulated from the surrounding atrial muscle. However, expression of connexins Cx43 and Cx45 increases toward the periphery, thus improving electrical coupling [29,30]. This extensively distributed system of atrial pacemakers constitutes the pacemaker complex, which induces changes in heart rate and beat-to-beat cycle length due to frequent exchange of dominance between different pacemakers. Changes in any of the input signals may produce pacemaker shifts and alter prevailing beat-to-beat rate and variability. The same is true when dynamic competition between multiple pacemakers rather than multifocal origin of the impulse, or different areas of the leading pacemaker location site, produce changes in BRV [30].

Beating Irregularity of Single Pacemaker Cells and Monolayer Cultures of Cardiac Cells
Jose and Collison [31] defined 'inherent heart rate' as the heart rate under the simultaneous presence of β adrenoreceptor antagonist propranolol and the muscarinic receptor blocker atropine. Rosenblueth and Simeone [32] proposed a formula that links inherent heart rate (HR 0 ) with sympathovagal inputs upon sinoatrial node: HR = m*n* HR 0 , where the product m*n represents sympathetic and vagal tone, respectively. Despite its simplicity, the above relation demonstrates a continuous cross-talk between the two components of heart rate control, ANS and in situ SAN pacemaker activity.
Pacemaker activity has been found to exhibit continuous variability of IBI (BRV) in different in vitro experimental studies. Clay and DeHaan [28] observed random fluctuations of CLs in spontaneously firing cells isolated from embryonic chick heart. They hypothesized that such fluctuations arise from inherent variability of diastolic depolarization towards AP threshold, due to stochastic open-close kinetics of membrane ionic channels. In addition, they found that SAN cells' beating irregularity diminishes with increasing aggregate cell size, whereas the co-efficient of variation is approximately inversely proportional to the square root of the number of interconnected cells' [28].
Jongsma, et al. [33] came to the same conclusion in terms of BRV in neonatal rat atrial and ventricular cardiac cells, with an IBI coefficient of variation around 2%. They found that IBI variability diminished when cells were interconnected by gap junction channels. Jongsma, et al. noted: 'Whereas individual pacemaker cells discharge irregularly, denervated SAN exhibits a regular discharge pattern. The hypothesis of explanation is that small SAN cells contain a few channels that give rise to observed voltage fluctuations, whereas in the intact SAN the membrane of combined cells are connected in parallel by resistive coupling between the cells. Due to the huge increase in the total number of ionic channels, individual variations are averaged out and firing threshold is reached at the same time in every cycle' [33]. Another possible influence upon BRV derives from surrounding atrial cells. Using mathematical models, Joyner and van Capelle [34] found that electrotonic interactions may play a significant role in beating rate irregularity of the intact SAN and Wilders and colleagues [35] revealed an increase in IBI coefficient of variation of around 2.9% induced by 'atrial load'. In the latter study, it was also proposed that a gradation in coupling resistance from the center of SAN towards the atrium is necessary for impulse propagation [35].
Kucera et al. [18] studied BRV monolayers with spontaneously beating neonatal rat ventricular myocytes, devoid of extracardiac influences, using unipolar extracellular electrograms. The beat-rate time series were examined with Fast Fourier transformation and Hurst scaling exponents' calculation, a mathematical method used for assessing fractal properties of different time series. In 21 from the 22 monolayers studied, the authors found a power law behavior with a β exponent around -1.35 [18], Typically, healthy subjects have a value lower than -1, while after myocardial infarction and after cardiac transplantation it is between -1.2 and -2 [16]. Lower or more negative β slope values correspond to time series with a closer resemblance to Brownian motion (for which β = -2) [11]. The findings of Kucera and colleagues support the hypothesis that 'during cardiac diseases, a partial break-down of autonomic control may unmask components of HRV intrinsic to cardiac tissue'. In addition, Kucera et al. [18] found that β slope during long lasting experiments displays a high standard deviation, suggesting that autonomic modulation of BRV might stabilize β slope.
Another study of Kucera and coworkers [19] demonstrated that both the movement of the beat-to-beat focus and variability exhibited power law behavior. Stochastic fluctuations in trans-membrane potential currents and gating of Ca 2+ release channels were not responsible for such long term correlations. However, turnover of various ion channels reproduced fractal properties responsible for intrinsic beating variability [19].
Yokogawa and Harada [20] evaluated power law behavior of beat rate time series in isolated atrial and ventricular myocytes, in order to assess inherent fractal properties of cell membranes irrespective of cell network interactions that arise within a monolayer culture. Timing of spontaneous contractions was determined from a sequence of phase contrast microscopic images, based on changes in brightness of pixels. Using a method called detrended fluctuation analysis (DFA) for the assessment of autocorrelation of the contraction timings they observed similarity between atrial and ventricular myocytes in terms of long-term correlations of beat rate fluctuations. According to the authors [20], their findings describe a universality of such long-term properties of beat rate fluctuations, since they were found in cells from different chambers of the heart. In a closely related study [36], the same authors proposed that fluctuations in ATP concentrations, reflecting the metabolic state of the cell, might be the cause of long-range correlations characterizing BRV, through activation of ATP-regulated K + P currents.
However, they concluded that this mechanism is only valid for neonatal ventricular cells that exhibit spontaneous activity, and not for mature ventricular myocytes. SAN pacemaker cells have different characteristics, such as type of ionic current channels, signal transducing proteins etc, limiting generalization of results in different experimental settings [29].
In addition, Wilders and Jongsma [37] observed that in isolated SAN node cells, the auto-correlation function of IBIs drops to 0 immediately after a lag of 1 beat. This means that there are no long range correlations and that the power spectrum is flat. Thus, there is no power law behavior in the in vitro single cell model of Wilders and Jongsma, which is at odds with previous studies in which a power law behavior was observed [18][19][20]. It seems that differences between cells from different chambers of the heart or electrotonic interactions between interconnected cells might change their fractal properties, reflected in power-law behavior.

Pacemaker Cells' Response to Autonomic Inputs
SAN pacemaker activity can be viewed as a transduction mechanism between input signals, such as neural, endocrine and mechanical factors and output signal that is the cycle length of pacemaker discharge [26]. Rocchetti et al. [38] investigated the effects of muscarinic receptor activation upon CL variability in isolated SAN cells of rabbits using the patch clamp technique. In today's cardiac cellular electrophysiology, patch-clamp, as described in detail by Kornreich [39], is the common technique to record the electrical activity of single cells [24,25].
Using the patch clamp technique, Rocchetti et al. [38] found that SAN cells exposed to different acetylcholine (Ach) concentrations displayed a significant short-term variability of CL, estimated with its standard deviation, which was increased exponentially in a dose-dependent fashion. This implies a non-linear relation between input/output signals at the level of SAN, where a constant level of receptor input can affect both the level and the variability of the output pacemaker discharge. According to Zaza and Lombardi [40], the reduction in CL variability in different pathologic states could reflect both a reduced input from ANS and a blunted response of the SAN to autonomic influences. Finally, any change in diastolic depolarization rate through alterations in different ionic currents of pacemaker cells might modulate CL variability. One possible explanation could be the shift of leading pacemaker between regions within the SAN with different receptor densities and sensitivity, associated with changes in neural activation, temperature or extracellular ionic concentrations. In addition, LF periodicities in the concentration of second messenger substances that can affect BRV could be attributed to the existence of mutually antagonistic signaling systems (e.g kinases versus phosphatases), which are altered by neural and chemical inputs [40].

From Cellular Dynamics to Whole Heart Preparations
There are a limited number of studies in the literature investigating origins of HRV in isolated hearts. In 1999, Langer and colleagues [21] were the first to quantify the ef-fects of afterload, preload, and temperature changes on IBI variability in intact isolated hearts of Sprague-Dawley rats. Using C 90 , the central 90% range of the beating intervals during 10 mins periods, as a time domain method of assessment of HRV and frequency domain methods, they found that changes in end-diastolic and aortic pressures had no effect, where decrease in temperature from 37 0 C down to 27 0 C increased C 90 about sevenfold. The beating intervals of isolated SAN cells were found to be normally distributed with significant variability. On the contrary, isolated rabbit right atrium and more pronouncedly, the intact rat hearts exhibited significantly smaller IBI fluctuations [21]. According to Langer and colleagues [21], at low heart rates, as typically occur during hypothermia, the steepness of the diastolic depolarization is reduced, leading to a variation of the trigger point that corresponds to the onset of AP.
The SAN includes clusters of pacemaker cells depolarizing independently with a multi-centric origin of a depolarization wave [29]. Since SAN cells exhibit the phenomenon of mutual synchronization, simulation studies from Michaels [41] demonstrated that a reduction in gap junction densities and conductivity may cause a decrease in intercellular connections, leading to a more prolonged and heterogeneous propagation of mutual synchronization. This decrease in the effective size of pacemaker cell clusters is responsible for the increase in IBI variability, estimated with C 90 , during hypothermia. In addition, individual BRV could be associated with different size of pacemaker clusters, implying that every factor that diminishes inter-node conductivity, such as low temperature, or accelerates conduction velocity, such as catecholamine application, can produce an increase or decrease in BRV, respectively [21]. Another finding from the study of Langer was that the SAN is the sole source of IBI variability, since neither preload or afterload conditions nor impulse conduction variation, which has been associated with the fractal structure of His-Purkinje network by Goldberger [42], contributed to the observed CLs variations.
More recently, Frey et al. [43] analyzed HRV in adult rabbit hearts, perfused for twenty minutes in a Langendorff set-up. An inherent HRV was found using metrics from the frequency domain (total power as an index of overall variability) in isolated denervated hearts, whereas both left atrial and ventricular chamber filling did not affect HRV. In another experimental study, Janousek et al. [44] compared HRV in ex vivo and in vivo whole heart preparations of rabbits and found a significant decrease in all frequency domain metrics of HRV in the first, compared to the second group. In addition, LF/HF ratio was five times higher in isolated heart than in in vivo hearts.
Monfredi and colleagues [abstract] tried to unmask possible mechanisms of HRV in isolated Langendorff perfused hearts of rabbits. They found that HRV, estimated with SDNN and power spectrum density, existed in denervated cardiac preparations and was significantly higher in isolated pacemaker cells than in whole hearts. Lack of effect of both atropine and propranolol that were infused in isolated hearts suggests the presence of HRV unrelated with possible residual effect of inherent autonomic ganglia. Moreover, HRV in terms of total power was decreased after administration of drugs which affects intracellular Ca 2+ homeostasis, such as ryanodine, implying that Ca 2+ release from the SR could be a possible mechanism of HRV generation. Finally, application of Cs + , frequently used in experimental laboratories as I f inhibitor [24,25], increased total power of heart rate signals. Table 1 summarizes major findings from studies discussed so far, regarding the presence of inherent beat-rate variability with power law behavior in different experimental in vitro and ex vivo studies.

Clinical Implications of Low HRV
The first large prospective population study that proved the significant prognostic value of low HRV after an acute myocardial infarction was the Autonomic Tone and Reflexes After Myocardial Infarction Study (ATRAMI), and included 1284 patients with a recent (<28 days) myocardial infarction [4]. A 24 h Holter recording was done to quantify HRV (using SDNN values) and ventricular arrhythmias. Low values of HRV (SDNN<70 ms) carried a significant multivariate risk of cardiac mortality. Furthermore, the association of low SDNN with left ventricular ejection fraction <35% carried a relative risk of 6.7, compared with patients with LVEF >35%. In the Framingham Heart Study [45], HRV time (SDNN) and frequency domain measures were computed in 736 patients and correlated with all-cause mortality, during 4 years of follow-up. The authors concluded that HRV offers prognostic information for mortality independent of that provided by traditional risk factors.
In the Zutphen study [46], 885 middle-aged (40-60 years old) and elderly Dutch men (aged 65-85) were followed from 1960 until 1990, whereas SDNN was determined from the resting 12-lead ECG. It was shown that low HRV is predictive of mortality from all causes, indicating that it can be used as an index of compromised health in the general popu-lation. It seems that the predictive value of low SDNN is independent of other factors, such as depressed left ventricular ejection fraction and presence of late potentials.
It is supposed that the change in the geometry of a beating heart due to necrosis may abnormally increase the firing of sympathetic afferent fibers by mechanical distortion of their sensory endings [47,48]. This excitation attenuates the vagal activity to the SAN. After acute myocardial infarction, HRV is reduced for a period of few weeks, while it is maximally, but not fully, recovered after 6 to 12 months [49].
In addition, retrospective ECG data analysis from 127 patients included in the Veterans Affairs' Survival Trial of Antiarrhythmic Therapy in Congestive Heart Failure (CHF) [50] demonstrated that CHF patients with SDNN <65.3 ms had a significantly increased risk of sudden death. Moreover, this study demonstrated that every 10 ms increase in SDNN conferred a 20% decrease in risk of mortality.
Furthermore, it has been demonstrated that subclinical inflammation and the concentration of inflammatory markers, such as cytokines, correlate strongly with cardiovascular mortality and morbidity in both healthy subjects and in those with known coronary artery disease (CAD) [51,52]. In this respect, a limited number of clinical studies investigating a possible association between ANS outflow and various inflammatory indices in patients with heart diseases have appeared in the literature [53][54][55].
Lanza and colleagues [53] assessed HRV and measured C-reactive protein (CRP) serum levels within 24 hours of admission in 531 patients with unstable angina pectoris. They found a small but statistically significant negative correlation between CRP levels and all HRV metrics derived from both time and frequency analysis, with the highest r coefficient with SDNN and VLF. After categorizing patients into 4 subgroups according to CRP quartile levels, significantly lower HRV (low SDNN) values were found in the upper CRP quartile. The subsequent multivariate analysis Table 1. revealed that SDNN and VLF were the most significant predictors of increasing CRP, whereas CRP itself was proven to be a strong predictor of impaired ANS activity as well. Psychari et al. [54] also reported a strong inverse association between CRP and several HRV indices (SDNN, HF and LF) in patients after acute myocardial infarction after adjustment for left ventricular function.

Summary of Different Experimental In Vitro and Ex Vivo Studies Concerning the Presence of Beat-rate Variability (BRV) and Power-law Behavior of Isolated Cells and Denervated Hearts
Malave et al. [55] examined HRV in relation to circulating levels of tumor necrosis factor-α (TNF-α), TNF receptors and norepinephrine in 10 controls, 15 patients with mild CHF and 14 subjects with moderate heart failure. They observed a significant inverse linear correlation between increased levels of all biomarkers and SDNN, LF and HF power among CHF patients [50]. In addition, the LF power was more closely correlated with circulating levels of TNF-α than was the HF component, whereas multiple linear regression analysis showed that TNF-α was a stronger predictor of reduced HRV than was the circulating levels of norepinephrine [55]. Therefore it was concluded that over-expression of TNF-α and subsequent loss of βadrenergic responsiveness contributes to the decrease in HRV observed in heart failure. According to recent experimental studies, TNF-α might inhibit β-adrenergic signal transduction through either activation of Gi proteins or impairment of activation of Gs proteins [56], something that could be viewed as an adaptive mechanism in the early stages of CHF through protection of cardiac myocytes from the deleterious actions of catecholamines. However, in the more advanced stages of the disease, this mechanism could become maladaptive, leading to a reduction in cardiac output [57].
Similar results in terms of alterations in HRV during septic shock and multiple organ dysfunction syndrome have been reported from different research groups [2,3,5]. For instance, using spectral analysis of HRV and blood pressure of critically ill patients, Goldstein et al. [5] was able to show that increased total variability and LF power were associated with recovery and survival, whereas a decrease in total power, LF/HF and LF power correlated with severity of illness and mortality in septic patients, 48 hours after being admitted to the Intensive Care Unit (ICU). This loss of variability of heart rate signals has been attributed to a 'decomplexification', which is associated with a defective communication between different organs, due to ANS dysfunction and parallels severity of disease [5,58]. Moreover, early HRV alterations could reflect different cytokine patterns as their local production can be significantly increased by denervation of an organ [59]. In this respect, Tateishi et al. [60] investigated the relationships between HRV and interleukin (IL) 6 upon admission in a cohort of 45 septic patients and they found that IL6 exhibited significant negative correlations with both LF and HF power values. These finding indicate an association between low HRV indices and hyper-cytokinemia, suggesting that there is a continuous cross-talk between ANS and immune-regulated inflammatory response.

HRV and I f in Critical Illness: Lessons from Cardiac Cellular Electrophysiology
The reasons for reduced HRV, estimated with time or frequency domain methods, which has been found in dif-ferent pathologic states, such as myocardial infarction [4] and severe sepsis or septic shock [5], have been debated and two theories have been developed. The first theory focuses on reduction of vagal tone and has been introduced by Akselrod et al. [9]. The second theory developed by Goldberger and colleagues [61] states that normal physiology has fractal-like properties with high levels of complexity that explain phenomena such as HRV. Severe disease reflects a 'decomplexification' mostly attributed to uncoupling between different restorative mechanisms [58]. Accumulating evidence from both in vitro and ex vivo experiments and human studies, with cardiac transplant recipients with hearts devoid from autonomic nervous inputs, support a potential third mechanism proposed recently by Griffin et al. [62], which is associated with an intracardiac origin of HRV. According to this hypothesis, SAN cells with an extreme heterogeneity in electrophysiological properties and intercellular connections of SAN cells can be viewed as an amplifier of various input signals [40]. During severe sepsis or in cardiovascular diseases, an unfavorable metabolic milieu could affect ion channel gating properties or membrane receptor densities, with significant impact upon level and variability of pacemaker activity. In addition, a possible reduced responsiveness of SAN cells to external stimuli could also negatively affect HRV [62].
In this respect, Janousek and colleagues [63] studied ex vivo heart preparations from rabbits, and using different protocols of ischemia and reperfusion, they computed thirty-five indices of HRV. They concluded that ischemia decreases significantly two metrics that seem to reflect intracardiac mechanisms of HRV, the high frequency component and SD2, which is the standard deviation of the Poincaré plot (RR n+1 versus RRn) along the line of the identity. The last method is used as a visual tool for variability assessment in different time-series. In a study of isolated SAN cells from control rabbits and rabbits with volume and pressure overload-induced heart failure, Verkerk et al. [64] found increased intrinsic CL and decreased diastolic depolarization rate, associated with a reduced density of I f . This is supported by findings in dogs, where heart failure results in down regulation of HCN4 and HCN2 expression in the SAN [29].
Administration of lipopolysacchraride (LPS, endotoxin derived from the cell wall of Gram-negative bacteria) has been found to decrease HRV [5,15]. In an animal study of experimental endotoxemia, Fairchild et al. [65] demonstrated a strong inverse correlation between SDNN and total power of RR time series and peak concentrations of different cytokines, 3-9 hours post-LPS. The same results were found after administration of recombinant TNF-α. Mechanisms responsible for decrease in HRV could be related with effects of LPS and/or cytokines on various ion channels.
Zorn-Pauly and colleagues [66] studied the effects of LPS on I f in human myocytes isolated from atrial appendages of patients undergoing open-heart surgery. After incubation with 10 µgr/mL of LPS, for 6 to 10 hours they found that the voltage-dependency of I f activation was shifted to more negative potentials. Consequently, I f was significantly lower in the physiological voltage range for SAN pace-maker activity. Using a mathematical computer model of a SAN cell, the effects of LPS on I f resulted in a longer CL. Using the mathematical SAN model cells, they further showed a reduced response of IBI fluctuations to ANS stimuli after LPS-induced I f impairment. A reduced responsiveness to external autonomic inputs may affect HRV. According to the authors 'altered I f gating during endotoxemia may originate from conformational changes, induced by direct binding of LPS to the channel and/or changes in membrane properties' [66]. The results of Zorn-Pauly et al. agree with the findings of Khaykin and colleagues [67], who showed that 10 mg of zatebradine, an I f blocker, reduced all frequency components of HRV in patients without structural heart disease. In addition, Joannides et al. [68] found that Ivabradien, a new selective I f current blocker at the level of SAN, significantly decreased LF/HF ratio in healthy volunteers during tilt and exercise.
However, in another experimental study [69], an increased HRV estimated with SDNN after administration of zatebradine was observed in conscious rats, three days after coronary artery ligation or sham-operation. Such results could support the hypothesis of Noble and DiFransesco [70] regarding the stabilizing effect of I f upon pacemaker frequency, 'against changes induced by altering other conductances whatever the magnitude of its contribution to the depolarizing current'. Reasons for such discrepancies could be differences in study design (early versus late administration, amount of dose, pre-existing pathology, or differences between species). In addition, I f blockers can affect HRV in different ways, either through a modulation of heart rate, effects upon baroreflex sensitivity in in vivo experimental settings through central action [71] or depending on levels of stress [69].
Recently, the BEAUTI f UL trial (MorBidity-mortality EvalUation of the I f inhibitor Ivabradien in patients with coronary disease and left ventricUlar dysfunction) including 10.917 patients with CAD, LVEF < 40% and heart rate ≥ 60 beats/min [72], has shown that Ivabradien in patients with a heart rate of >70 beats/min decreased admission to the hospital for fatal and non-fatal myocardial infarction or coronary revascularization [73]. Moreover, the SHIFT trial (Systolic heart failure treatment with I f inhibitor Ivabradien trial) including patients with stable symptomatic heart failure, LVEF ≤ 35% and heart rate ≥ 70 beats/min, showed a significant reduction in deaths from heart failure after heart rate reduction with Ivabradien [74].
Effects of LPS upon heart rate in septic patients seem to involve rather inherent heart rate than ANS influences at the level of SAN [75]. In different studies it has been shown that tachycardia in sepsis develops slower and lasts longer compared to reduction in blood pressure, or even continues after adequate volume resuscitation. In addition, septic tachycardia is apparent even after a decline in the plasma levels of catecholamines, indicating a no baroreceptor mediated mechanism [76,77]. Finally, sepsis seems to increase sympathetic tone at the level of peripheral nerve endings independent of increases in firing heart rate [78].For these reasons, it has been suggested that septic tachycardia is due to an increase in in vitro inherent heart rate [79].
This non-autonomic origin of increased heart rate in sepsis is paradoxically associated with desensitization of beta-receptors [80], whereas analysis of HRV in septic animals by Goldstein et al. [81] has found a decrease or absence of sympathetic input upon the heart, estimated with LF/HF and LFnu, despite a decrease in overall HRV. However, these experiments do not elucidate whether loss of HRV is related to an endotoxin effect at the level of ANS output, baroreflex sensitivity or the pacemaker cell itself [79]. In an experimental study of single intraperitoneal LPS administration of 5 mg/Kg in rats, Wearden [82] showed that inherent heart rate increased at approximately 4 hours after endotoxin injection, whereas LFnu and LF/HF ratio did not change significantly. In addition, total power of HRV began to decrease at 2 hours, reaching its minimal values 7 hours later. According to Wearden [82], loss of HRV estimated with total power of frequency spectrum, reflects baroreceptor mediated mechanisms of heart rate variability, whereas increased inherent heart rate might be due to an unknown cellular mechanism. Moreover, this mechanism could initiate an LPS-induced septic tachycardia, however, catecholamines are required to maintain this effect. The fact that tachycardia in the absence of significant sympathetic input upon the heart, as it was estimated with HRV analysis, paralleled a decrease in contractility, could reflect a maladaptative mechanism responsible for an increase of oxygen delivery to the tissues with less autonomic stimulation [82].
Time changes of HRV and inherent heart rate during sepsis might reflect study designs. In different in vivo experimental septic models that include different doses or routes of administration of LPS (i.e., intravenous bolus, cecal ligation or puncture models), it has been shown that there is a varying time lag after injection of endotoxin before endotoxemia develops [83]. This may be due to different time of peaks in the concentration of various cytokines [84]. However, different experimental conditions, i.e., animal's weight, temperature, anesthetic agents, residual effects of catecholamines or even doses of LPS can bias such findings significantly [86][87].
The finding that I f could also mediate HRV changes during endotoxemia [66], might explain findings from other studies concerning both septic tachycardia and loss of overall HRV [75,81]. Although, the combination of a markedly attenuated vagal tone that is found in severe sepsis [79], a massive release of endogenous catecholamines or their exogenous administration could override in vivo the bradycardic effects of LPS, endotoxin can also sensitize the HCN channels for sympathetic stimulation, thereby increasing heart rate [66]. As a consequence, a combination of increased heart rate and narrowed HRV indicate autonomic dysfunction and poorer prognosis in patients with multiple organ dysfunction syndrome , who display similar findings in terms of pathophysiologic mechanisms, with those suffering from chronic heart failure [88].
The above association is based on the presence of sympathetic overactivity, autonomic dysfunction, inappropriately increased heart rate, insulin resistance and in some cases, cardiomyopathy with reduced cardiac contractility [88,89]. As with heart failure, a down regulation of beta receptors is also found in septic patients, limiting the possible use of beta blockers [80]. In addition, increased heart rate in these patients has been associated with fatal outcome [90]. However, in cardiac patients sympathetic activity dominates over vagal tone where in MODS patients both branches of ANS are attenuated [79]. For these reasons, it has been hypothesized that during critical illness except for ANS impairment, a defective signal transduction at the level of pacemaker cells could also account for observed differences between cardiac and septic patients. I f could be a possible candidate due to its reduced amplitudes during endotoxemia. In this respect, a possible beneficial effect of I f inhibitors merits further investigation and could shed more light into the role of I f in septic cardiomyopathy. A randomized, controlled phase-2 trial -the MODI f Y trialis underway (www. clinicaltrials. gov, NCT 01186783), trying to explore a possible benefit from the use of Ivabradien in 70 critically ill patients with multiple organ dysfunction syndrome, with heart rate ≥ 90 beat/min and contraindication to beta-blockers [91].
There are a few more studies in the literature, trying to investigate a possible impact of LPS upon HRV, at the cellular level. In monolayers of spontaneously contracting rat neonatal ventricular cardiomyocytes, Schmidt and colleagues [92] showed that application of 1 µgr/mL of LPS significantly reduced SDNN, RMSSD and pNN50. Since the authors did not study real SAN cells, they hypothesized that endotoxemia induces an ionic remodeling, based on already proven negative effects of LPS upon different ion channels, such as a decrease in density and activation properties of L-type Ca 2+ channels and down regulation of β-adrenergic receptors [79,89].
Recently, Wondergem et al. [93] confirmed the previous results in immortalized HL-I cardiomyocytes. Application of 1 µgr/mL of LPS resulted in a rapid reduction in Ca 2+ oscillations and total intracellular Ca 2+ concentration. In addition, I F was inhibited at very negative potentials. The effects occurred rapidly after LPS application (within minutes), which argue against the contribution of release of second messengers, such as cytokines and cAMP.
One more study that was performed in Human Embryonic Kidney 293 (HEK293) cell lines, confirmed negative effects of LPS on I f properties. Klöckner et al [94] transfected HEK cells with cDNA encoding the HCN2 channel proteins, for comparison reasons with Zorn-Pauly findings [66], since these HCN transcripts are predominantly expressed in human right atrium myocytes [95]. LPS was found to induce NF-kB activation through Toll-likereceptors-4 (TLR-4), leading to IL-6 production within 1-2 hours that reached its maximum after more than 24 hours. However, decrease in activation of I f with reduced steepness of diastolic depolarization seemed to occur within approximately 8 seconds. As a result, it was shown that the polysaccharide part of LPS 'O-chain' is required for immediate reduction of HCN channel activity. Thus, endotoxin seems to affect with different structures of its molecule the activation of the host immune system and the modification of HCN channel electrophysiological properties. According to Klöckner et al [94], their results parallel reports from the anesthesia literature, where both propofol and halothane induced a dose-dependent inhibition of HCN channel activity and an associated decrease in heart rate variability [96,97]. However, possible differences of biophysical properties of various HCN transcripts could bias previous results, since these isoforms have been found to coassemble in varying ratios, thus modulating SAN responsiveness to different stimuli [98]. Table 2 presents some of the studies investigating the role of funny current in cardiovascular diseases and sepsis, in both experimental and clinical settings.

CONCLUSIONS AND FUTURE SUGGESTIONS
Contradictory results in terms of HRV changes during severe sepsis and MODS have been found, mainly due to differences in experimental design, (i.e. differences in species studied, isolated SAN cells versus intact SAN preparations and/or intact hearts, etc). In addition, sepsis and MODS are considered extremely complex and available experimental models fail to represent accurately human diseases. Although in cardiovascular medicine different experimental and clinical studies have managed to confirm the value of HRV measurement in the clinical setting, critical illness still remains an open field for future investigation, concerning implementation of HRV quantitative analysis and evaluation of its potential for risk stratification or even building different prognostic models. In addition, understanding of the complex mechanisms responsible for its generation requires in vitro or ex vivo models, representative of different diseases, in order to develop future therapies for cardiovascular and septic patients, targeting increased heart rate.
Existing literature supports the presence of intracardiac mechanisms of heart rate variability at the level of either membrane channel kinetics or intracellular transduction signaling. However, and despite evidence from numerous studies that demonstrate negative inotropic effects of LPS during sepsis [99], there are only a limited number of reports investigating its chronotropic effects, particularly in association with inherent heart rate and loss of beat-to-beat rate variability in animal septic models.
The complexity of HRV and the existing debate regarding its origin puts another obstacle for understanding underlying pathophysiology and assessing its value as a monitoring tool, in terms of both prognosis and treatment effect, in different groups of patients. Many details remain to be discovered, such as effects of LPS on conductivity between SAN pacemaker cells, since recent findings support a reduced expression of various gap junction proteins in ventricular cells during severe sepsis [100]. In addition, association of sepsis with pacemaker shifts, ionic remodeling of different currents responsible for spontaneous diastolic depolarization or membrane potential propagation within the SAN, could be some interesting areas for future research, in relation with reduced HRV as an index of cellular metabolic stress. Moreover, the use of simulation studies with mathematical models that are more valid for larger mammal's SAN electrophysiological properties might be of significant value in the near future.
Finally, standardization of experimental protocols and methods used for the estimation of HRV and the use of open-source data and knowledge bases, such as the website physionet (www. physionet.org) [6], with different analytic tools that can be adopted for free by independent investigators, is an ongoing effort that could enhance understanding and implementation of this technology in both basic and clinical research.

CONFLICT OF INTEREST
The authors confirm that this article content has no conflicts of interest.  Verkerk et al. [64] SAN cells from control rabbits versus rabbits with volume and pressure overload-induced heart failure Increased intrinsic CL and decreased diastolic depolarization rate, associated with a reduced density of If.

ACKNOWLEDGEMENTS
Zorn-Pauly and colleagues [66] Effects of LPS on If in human myocytes isolated from atrial appendages of patients undergoing open-heart surgery.
Voltage-dependency of If activation was shifted to more negative potentials.
Reduced response of IBI fluctuations to ANS stimuli after LPS-induced If impairment Joannides et al. [68] Healthy volunteers during tilt and exercise. Ivabradien, a new selective If current blocker at the level of SAN, significantly decreased LF/HF ratio BEAUTIfUL trial [72,73] 10.917 patients with CAD, LVEF < 40% and heart rate ≥ 60 beats/min Ivabradien in patients with a heart rate of >70 beats/min decreased admission to the hospital for fatal and non-fatal myocardial infarction MODIf Y trial [91] Effects of Ivabradien in 70 critically ill patients with multiple organ dysfunction syndrome, with heart rate ≥ 90 beat/min and contradiction to beta-blockers (still ongoing) Schmidt and colleagues [92] Impact of LPS upon HRV (SDNN, RMSSD), in monolayers of spontaneously contracting rat neonatal ventricular cardiomyocytes 1 µgr/mL of LPS significantly reduced SDNN, RMSSD and pNN50 Klöckner et al [94] Transfected