This Article Abstract Full Text (PDF) Supplement All Versions of this Article: biophysj.106.091892v1 92/10/3734    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Ponard, J. G. C. Articles by Kucera, J. P. Search for Related Content PubMed PubMed Citation Articles by Ponard, J. G. C. Articles by Kucera, J. P.

Mechanisms of Intrinsic Beating Variability in Cardiac Cell Cultures and Model Pacemaker Networks

Department of Physiology, University of Bern, Bühlplatz 5 CH-3012 Bern, Switzerland

Correspondence: Address reprint requests to Jan P. Kucera, MD, Dept. of Physiology, University of Bern, Bühlplatz 5 CH-3012 Bern, Switzerland. Tel.: 41-31-631-87-59; Fax: 41-31-631-46-11; E-mail: kucera{at}pyl.unibe.ch.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Heart rate variability (HRV) exhibits fluctuations characterized by a power law behavior of its power spectrum. The interpretation of this nonlinear HRV behavior, resulting from interactions between extracardiac regulatory mechanisms, could be clinically useful. However, the involvement of intrinsic variations of pacemaker rate in HRV has scarcely been investigated. We examined beating variability in spontaneously active incubating cultures of neonatal rat ventricular myocytes using microelectrode arrays. In networks of mathematical model pacemaker cells, we evaluated the variability induced by the stochastic gating of transmembrane currents and of calcium release channels and by the dynamic turnover of ion channels. In the cultures, spontaneous activity originated from a mobile focus. Both the beat-to-beat movement of the focus and beat rate variability exhibited a power law behavior. In the model networks, stochastic fluctuations in transmembrane currents and stochastic gating of calcium release channels did not reproduce the spatiotemporal patterns observed in vitro. In contrast, long-term correlations produced by the turnover of ion channels induced variability patterns with a power law behavior similar to those observed experimentally. Therefore, phenomena leading to long-term correlated variations in pacemaker cellular function may, in conjunction with extracardiac regulatory mechanisms, contribute to the nonlinear characteristics of HRV.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
The normal heartbeat is initiated in the sinoatrial node (SAN), a pacemaking structure that generates action potentials which propagate across the atria and ventricles and trigger their rhythmic contractions. Although the heartbeat is a periodical phenomenon, heart rate varies continuously with time. These variations adjust the cardiac output to the demand of the organism and are essentially controlled by regulatory reflexes mediated by the autonomic nervous system (1–3).

During the past decades, numerous efforts have been invested to extract clinically relevant information from the analysis of the variations of heart rate on a beat-to-beat basis (4,5). A widespread method for the analysis of heart rate variability (HRV) is to decompose the variations of heart rate into a spectrum of distinct frequency components by means of Fourier analysis or autoregressive modeling (6–8). When applied to human subjects, this analysis reveals two frequency ranges in which periodic variations in heart rate prevalently occur: a high frequency range (centered around 0.25 Hz) and a low frequency range (centered around 0.1 Hz). It is well established that these periodic variations are related to breathing and the regulation of arterial blood pressure, respectively, and that they are mediated by the autonomic nervous system (9,10). The quantitative analysis of the high frequency and the low frequency components is thus widely applied to evaluate the balance between the sympathetic and the parasympathetic branches of the autonomic nervous system in a variety of clinical settings (4,11).

However, the major component of HRV occurs at frequencies below the low frequency range (<0.04 Hz). It is known that in the range 10^–4–10^–2 Hz, the power spectrum of HRV does not exhibit periodic components but exhibits a so-called power law behavior (12–15). The term "power law" describes the relationship between two quantities x and y if this relationship can be expressed as y = A x x^k or y x^k, where k is called the "power law exponent". Power laws are among the most frequent laws that describe fractal objects, i.e., objects that have a fine structure at each scale and appear similar at all scales of magnification (16). In a time series (e.g., of heart rate or interbeat intervals), this property of self-similarity reflects the presence of long-term correlations across a wide range of temporal scales (17,18). In the context of HRV, the notion of power law behavior refers to the power law relating the power spectral density (PSD) of HRV and frequency (f), expressed as PSD(f) f^ß, where ß is the power law exponent. The power law, reflecting the fractal properties of HRV, can be graphically visualized in a double logarithmic plot of PSD(f) versus f, in which the spectrum follows a line with a slope corresponding to ß. In healthy subjects, the exponent ß is close to –1.

In clinical studies of HRV, it has been shown empirically that certain disease states are associated with significant modifications of the power law exponent ß. In particular, a more negative ß (a steeper slope of the power spectrum) has been found in myocardial infarction patients and in patients after cardiac transplantation (19). It has also been shown that a more negative ß occurs in patients prone to ventricular fibrillation (20) and discriminates patients with increased mortality (21). Thus, ß represents a predictive risk-stratification marker, as more negative values of ß are associated with an increased probability of sudden cardiac death (19,20,22–24). In other studies, HRV was evaluated using detrended fluctuation analysis, a fractal analytical method to quantify long-term correlations analogous to the determination of the slope ß (25). These studies demonstrated alterations of detrended fluctuation parameters in, e.g., sleep apnea patients (26) and in asymptomatic relatives of patients with familial dilated cardiomyopathy (27).

These findings indicate that the analysis of HRV using fractal methods represents a clinically relevant evaluation tool. However, to date, the physiological mechanisms underlying the power-law behavior of HRV have scarcely been studied. To obtain the full benefit from these analyses, it is, therefore, essential to gain understanding of the underlying physiology to improve the specificity of risk stratification during clinical evaluation and the development of directed patient follow-up protocols and specific therapies.

It is an accepted notion that the aperiodic variations characteristic of HRV at low frequencies are the result of complex nonlinear interactions between extracardiac regulatory mechanisms (autonomous, hemodynamic, and endocrine) (28–30). However, in a previous study, we demonstrated that beat rate variability (BRV) in cardiac cell cultures (representing an in vitro model of a spontaneously active network of coupled cells) also follows a power law behavior in the absence of any extrinsic regulatory mechanisms (31). This suggests that not only extracardiac influences but also mechanisms intrinsic to cardiac tissue contribute to the power law behavior of HRV observed in vivo. Within an isolated network of pacemaker cells, BRV might be determined by the complex nonlinear dynamics of interactions among a large number of cell oscillators. However, it might also be determined by the integration of stochastic phenomena occurring at the cellular level.

To gain deeper insight into the phenomena underlying the power law behavior of HRV, our study pursued the following aims: 1), to analyze the spatiotemporal patterns of spontaneous electrical activity on a beat-to-beat basis in cultures of neonatal rat ventricular myocytes, 2), to assess whether spontaneous activity in these cultures shares ionic mechanisms similar to those in the SAN, and 3), to investigate, using mathematical modeling of SAN cells and networks, whether stochastic phenomena occurring at the cellular level may contribute to beating variability and/or explain the characteristics of BRV. Specifically, we sought to determine the effects i), of fluctuations of transmembrane currents (as caused by the stochastic gating of sarcolemmal channels), ii), of the stochastic gating of Ca²⁺ release channels in the sarcoplasmic reticulum (SR), and iii), of the dynamic turnover of ion channels in the cell membrane.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Patterned growth cardiomyocyte cultures on microelectrode arrays
Monolayer cultures and patterned disk-shaped cultures of neonatal rat ventricular myocytes (Wistar; 1–2-day old; Zentrale Tierställe, University of Bern, Bern, Switzerland) were prepared and grown on microelectrode arrays (Sensors, Actuators and Microsystems Laboratory, University of Neuchâtel, Neuchâtel, Switzerland) following previously published photolithographic procedures (31–33). Dissociated ventricular myocytes were preplated for 2 h in culture flasks to reduce the percentage of nonmyocardial cells. Before seeding the dissociated myocytes, the microelectrode arrays were mounted in custom-fabricated chambers (diameter 12.3 mm; volume 1 ml), incorporating a connection interface for the microelectrode arrays. The cardiomyocytes were seeded at a density of 1600–2100 cells/mm². The preparations were incubated at 36°C and 0.9% CO₂ in medium 199 with Hanks' salts (M-199, Sigma-Aldrich, Buchs, Switzerland) supplemented with 5% neonatal calf serum (Bioconcept, Allschwil, Switzerland), penicillin (20,000 units/L; Fakola, Basel, Switzerland), streptomycin (34 µmol/L; Fakola), vitamin B12 (15 µmol/L; Sigma-Aldrich), bromodeoxyuridine (100 µmol/L; Sigma-Aldrich), and epinephrine (10 µmol/L; Sigma-Aldrich).

Recording of extracellular electrograms
Electrophysiology experiments were performed in 3–4-day-old cultures (31). Before an experiment, the culture medium was renewed. The chambers were then connected to a custom-designed amplifier array (gain: 1000x) and placed back into the incubator. A period of 60 min was allowed for equilibration of temperature (36.0°C ± 0.1°C). Control measurements of temperature using a thermocouple probe (Fluke 80TK, Wallisellen, Switzerland) were performed to verify that fluctuations of temperature were <0.1°C (peak-to-peak). Extracellular electrograms (V_e) were recorded from an array of 61 indium-tin oxide microelectrodes (diameter 40 µm each), arranged along a hexagonal lattice (layout shown in Fig. 1 A, interelectrode distance (IED) = 1.425 mm). The amplified signals were sampled at 10 kHz (Data Acquisition Processor 4400a, Microstar Laboratories, Bellevue, WA) and analyzed using custom-written software.

View larger version (40K): [in this window] [in a new window]   FIGURE 1  Recording of spontaneous activity of cardiac cell cultures and data analysis. (A) Microelectrode array layout. The dots indicate the positions of the microelectrodes. The circle denotes the culture area. (B) Representative extracellular electrogram (Ve) recorded during spontaneous activity (top) and its first time derivative (dVe/dt, bottom). The red line marks activation time (occurrence of dVe/dtmin). (C and D) Top: Signals recorded during two successive beats normalized to signal amplitude. Vertical red markers denote dVe/dtmin. Bottom: Corresponding isochrone maps and localization of the focus of activation. The hexagon denotes the area mapped by the electrodes. The colored map is based on triangulation and linear interpolation of activation times. The overlaid solid (C) and shaded (D) circles represent isochrones (interval: 2 ms) based on the surface fit described by Eq. 1. The central dot denotes focus position based on this fit. (E) Overlay of the fitted maps for the beats shown in C (solid) and D (shaded). (F) Fitted focus coordinates (origin at the center of the array) and computed standard deviations SDx and SDy (error bars) determined based on the standard deviation of the residuals for the spontaneous activations shown in C (solid) and D (shaded), respectively (see Materials and Methods for details concerning the analysis). (G) Top: Interbeat interval series (20 successive beats) for the corner and central electrodes (colors corresponding to C and D). Bottom: Corresponding focus coordinates. Solid and shaded dashed lines mark the beats depicted in C and D, respectively. (H) Evenly resampled time series (2 Hz; used for spectral analysis), corresponding to the data shown in G.

For every experiment, the acquired data and all fiducial markers of dV_e/dt_min (detected automatically using a custom-written computer program) were carefully inspected visually and edited to ensure that all activations had been detected and to remove erroneous detections.

The top diagrams in Fig. 1, C and D, represent electrograms, fiducial markers, and corresponding isochrone activation maps for two successive beats. Isochrone maps were constructed based on triangulation and linear interpolation of the activation times corresponding to a given propagated action potential (34,35).

The origin of electrical activation, which we hereafter term "focus", was identified by a gradient-expansion iterative least squares surface fit of a hyperboloid function in the (x,y,t) space (Eq. 1 below) to the activation times t_act(x,y) using the Levenberg-Marquardt method (36,37):

(1)

where is conduction velocity at infinite distance from the focus, t_focus is the activation time at the focus, (x_focus,y_focus) are the coordinates of the focus, and is the hyperbolic term (we observed that the least squares error is smaller and the fit more accurate if this term is included; if = 0, the function represents a cone in the (x,y,t) space). Because the cells were seeded on a uniform substrate, conduction velocity was assumed to be isotropic. This method to locate the focus is illustrated in Fig. 1, C and D, where isochrones based on this surface-fitting procedure are represented by concentric circles overlaid on the color-coded isochrone maps and where the central dots correspond to the identified foci. The least squares fits were assumed to be successful when the ² test associated with the iteration process converged in <100 iterations. Because this analysis was convergent only if enough operational electrodes were available for a reliable determination of activation times, we excluded experiments in which >50% of the electrodes were not utilizable because of deficient connections or low signal/noise ratios. The number and location of the electrodes for which accurate focus localization was essential are described in the Supplementary Material.

To assess the applicability of Eq. 1 to describe the spread of focal activation, we computed the residuals between the measured activation times and the activation times fit with Eq. 1. We repeated this procedure for every focal activation in all experiments and applied the Kolmogorov-Smirnov test to evaluate whether the residuals differ from a Gaussian distribution. For all the focal activations analyzed, the hypothesis that the distributions of the residuals follow a normal distribution was never rejected (at p < 0.05 level), indicating that Eq. 1 is adequate to describe the spread of focal activation. We also computed the correlation coefficient of the fit (r), which was on average >0.99.

The two-dimensional (2-D) fits represented by circles in Fig. 1, C and D, are overlaid in Fig. 1 E. The fact that the two fitted maps are not concentric is suggestive of a displacement of the focus. Since conduction is never absolutely uniform in a biological preparation (even in a cell culture that is microscopically homogeneous) and since nonuniform conduction can affect the precision of the determination of x_focus and y_focus, we quantified conduction nonuniformity (CNU) and evaluated the influence thereof on the accuracy of focal point localization. The nonuniformity of conduction was reflected by the fact that the standard deviation of the residuals (SD_resid) was in the range of a couple of milliseconds and thus several times larger than the precision of activation time determination (see above).

SD_resid thus provides quantitative insight into CNU at the spatial scale corresponding to the resolution of the system (IED = 1.425 mm). To quantify CNU at this scale, we computed the ratio between SD_resid and the time required by the wavefront to cover the IED (conduction time) as CNU = SD_resid/(IED/), where is the conduction velocity obtained from the fit to Eq. 1 and IED/ corresponds to the expected conduction time. CNU is a unitless parameter (expressed as percentage) and can be interpreted as a coefficient of velocity variation. A similar method has been used previously to quantify CNU in cultured cell strands (38). In this study, CNU ranged from 6% to 20%.

SD_resid was subsequently used to quantify the standard deviation of the fitted focus coordinates x_focus and y_focus (for details, see Press et al. (37)). The computed standard deviations of focus coordinates (SD_x and SD_y), therefore, characterize the error in the localization of the focus, which is essentially due to the nonuniformity of conduction intrinsic to the biological preparations. In the two activations depicted in Fig. 1, C and D, SD_resid was 0.84 and 0.74 ms, CNU was 14.3% and 12.8%, and r was 0.9926 and 0.9938, respectively. As shown in Fig. 1 F, SD_x and SD_y amounted to 43 and 47 µm for the beat illustrated in Fig. 1 C, and to 38 and 40 µm for the beat depicted in Fig. 1 D, for which the computed focus location shifted by 130 µm.

The determination of SD_x and SD_y was verified and validated by Monte Carlo simulations of synthetic data sets (37) in which a perturbation having a normal distribution with a standard deviation equal to SD_resid was added to the fitted activation times at each electrode. These simulations confirmed the computed SD_x and SD_y with a precision of 1%. The statistical probability that focus coordinates were equal for the two beats depicted in Fig. 1, C and D, was p = 0.028. This example therefore illustrates the limit of our algorithm to discriminate between two different foci.

This analysis yielded similar results when applied on different beats in different experiments and thus shows that changes in computed focus position several times the spatial SD_x and SD_y can be considered significant, although much smaller than the spatial resolution of the microelectrode array. However, accurate assessment (SD_x and SD_y < 100 µm) of the focus position was possible only when the focus was located within the area mapped by the electrodes and the algorithm permitted only an approximate localization when the focus was outside this area. To confirm that the system and the analysis algorithm are able to accurately identify the origin of focal electrical activation in monolayer cultures, we paced the preparations from defined locations. Subsequent analysis of the focus of activation pinpointed the centers of the stimulation electrodes with a precision of 50 µm (see the Supplementary Material for details).

Fig. 1 G illustrates interbeat interval series for the six corner electrodes and the center electrode as well as focus coordinates for 20 successive beats. Interbeat interval variations (in the range of milliseconds) were significantly larger than the maximal error of their determination (0.3 ms). The crossing of the different traces reflects the movement of the focus on a beat-to-beat basis.

Beat rate (BR) and the PSD of BRV were computed as described previously (31), using established methods applied for the frequency domain analysis of HRV in vivo (4). In brief, evenly sampled series (2 Hz) of instantaneous BR were constructed by linear interpolation from the series of activation times, as shown in Fig. 1, G and H. The PSD of BRV was then computed from the fast Fourier transform of the resampled BR series of one electrode. The power law exponent of the PSD was computed as the slope of the linear regression of log(PSD) versus log(frequency). The discrete Fourier PSD data points were given a weight inversely proportional to their density on the logarithmic abscissa (31). A power law was assumed to be present if r > 0.85. Similarly, evenly sampled series (2 Hz) of focus coordinates (see Fig. 1, G and H) were constructed by linear interpolation and were used to compute and characterize the power spectra of these coordinates.

Experimental protocols
The preparations were inspected using phase contrast microscopy to verify their structural intactness (i.e., uniformity of cell density, absence of localized damage). Structurally nonuniform preparations were discarded.

Experiments in which the spatiotemporal characteristics of beating variability were investigated were conducted in disk-shaped cultures (diameter: 11 or 12.3 mm) after enhancing sustained spontaneous electrical activity with 0.4 µmol/L Bay K 8644 (31). After an equilibration period of 60 min, the electrical activity was recorded continuously for 30 min. To verify that fluctuations of temperature do not contribute to BRV under our recording conditions, control experiments (n = 5 different cultures) were performed to evaluate what temperature difference (T) is needed to produce a given change in BR (BR). In these experiments, electrical activity and temperature were continuously recorded and temperature was increased from 34°C to 38°C and decreased from 38°C to 34°C by modifying the set point of the incubator. The sensitivity of BR to temperature, determined as BR/BR_average/T during the linear phases of temperature change, amounted to 5.1% ± 1.7%/°C (range 2.3–7.5%/°C). In the other experiments presented in this study, BR variations were in the range of 2.5% and occurred rapidly over timescales of a few seconds. If they had been caused by temperature changes, such BR variations would have necessitated fast temperature variations in the range of 0.5°C. Thus, the contribution of temperature fluctuations (measured to be <0.1°C in the incubator) to BRV was minimal.

Experiments in which the effect of pharmacological agents (nifedipine (Sigma-Aldrich), Bay K 8644, ryanodine, or ZD7288 (Tocris, Bristol, United Kingdom)) on BR was tested were performed on monolayer cultures covering the entire area at the bottom of the chambers. After a control recording (10 min), the compounds were added using a custom-built injection setup. Their effect on BR was assessed by comparing BR during 10-s periods immediately before and 1.5 min after injection, respectively. Control experiments were performed with vehicle medium to verify that the injection itself did not induce significant modifications of BR. The injection led to a minimal decrease in BR by 2% ± 1% (n = 4), possibly caused by a slightly lower temperature of the injected vehicle.

Mathematical modeling of pacemaker cells and networks
The Wilders-Jongsma-van Ginneken (WJvG) model of the rabbit SAN cell (39) was used to simulate the electrical activity of spontaneously active cells and cell networks. The model takes into account 10 transmembrane currents, intracellular compartmentation, and changes of [Na⁺]_i, [K⁺]_i, and [Ca²⁺] in different compartments. In cell networks, intercellular gap junctional conductance was set to 2.58 nS, following values reported previously (40).

Membrane potential, gating variables of ion currents, and ion concentrations were integrated explicitly using a constant time step (dt) of 0.005 ms. Activation times were defined for every cell as the time when membrane potential reached –30 mV during depolarization. In the networks, focus position was defined for every activation as the location of the earliest activated cell.

Stochastic fluctuations of transmembrane current
The effects of fluctuations of transmembrane current due to stochastic open-close kinetics of ion channels was investigated by introducing into the model a fluctuation current, I_fluct, with zero mean and no long-term correlations. I_fluct was constructed as a stochastic Ornstein-Uhlenbeck process (41) by iteration of the following equation:

(2)

with ß the rate constant determining the extent of temporal correlations, the standard deviation of the fluctuation current, and (t) a normally distributed random variable with zero mean and unitary standard deviation. We set ß = 0.7 ms^–1, which resulted in I_fluct with an autocorrelation function that decreased below 0.05 for time lags >4 ms. This reflects the notion that ion channels remain open for a few milliseconds and have no long-term memory of their history of conformational states.

Stochastic gating of ryanodine receptors
To test the contribution of the stochastic open-close kinetics of the SR Ca²⁺ release channels (ryanodine receptors (RyRs)) to beating variability, we modified the original formulation of the Ca²⁺ release current, I_rel, in the WJvG model given as (39)

(3)

by the following formulation:

(4)

with P_o,RyR the fraction of open RyRs, [Ca²⁺]_Rel the Ca²⁺ concentration in the release store, [Ca²⁺]_i the myoplasmic Ca²⁺ concentration, A_rel the maximal permeability (permeability of all RyRs together), and k_m,Ca the binding constant of myoplasmic Ca²⁺ to RyRs. The term ([Ca²⁺]_Rel – [Ca²⁺]_i) represents the [Ca²⁺] difference between the release store and the myoplasm. The gating of individual RyRs and thus the fraction of open RyRs (P_o,RyR) was computed as a stochastic Markov process following a two-state model with a second-order dependence on [Ca²⁺]_i. Following Sobie et al., we incorporated a constant closing rate, ß, of 0.5 ms^–1 (42). The Ca²⁺ dependence of the opening rate, , was formulated as

(5)

At steady state (constant [Ca²⁺]_i), P_o,RyR is thus /( + ß) = 1/[1 + (K_m,Ca/[Ca²⁺]_i)²] and corresponds to the third term of the original I_rel formulation. Following Bers (43), we assumed a number of 5000 RyR channels in a cell (number required to attain a typical peak [Ca²⁺]_i during Ca²⁺-induced release in a cell having the volume of the WJvG cell (5 pL)). At each time step dt of the simulations (satisfying dt and ßdt << 1), for each closed RyR (respectively, open RyR), one random number, A, was drawn from the uniform distribution [0,1]. If A < dt (respectively, A < ßdt) then the channel was opened (respectively, closed).

Ion channel turnover
Incorporation of channels into the plasma membrane and their internalization from the membrane (turnover) were formulated as a stochastic "immigration-death" model, in which the incorporation and the internalization of channels were formulated as Poisson processes (41) for the two principal ion channel types of the WJvG model, the L-type Ca²⁺ channels, and the delayed rectifying K⁺ channels.

The processes were formulated based on the average number, N_average, of channels in the plasma membrane (10,000 and 1000, respectively (44)) and the half-lives of the channels (16 and 4 h, respectively (45)). If t_1/2 is the half-life of a given type of channel, then the recycling rate () is related to t_1/2 by = (ln 2)/t_1/2. To maintain N_average channels in the membrane, channels must therefore be incorporated at a rate, , which is N_average times larger, i.e., = N_average.

At each time step, dt, of the simulation (satisfying dt and dt << 1) for each cell and for each channel type, one random number, A, was drawn from the uniform distribution [0,1]. If A < dt then one channel was incorporated. Then, for each individual channel, one random number, B, was drawn from the uniform distribution [0,1]. If B < dt then the channel was removed from the cell membrane.

Statistics
Values are given as mean ± SD. Data were compared using the Wilcoxon signed rank test or the Student's t-test, and p < 0.05 was considered significant.

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Origin of spontaneous activity in disk-shaped cultures
To characterize the patterns of spontaneous activity in the neonatal rat cardiomyocyte cultures, we first investigated whether this activity originates from foci located at preferential locations. Sustained activity was induced in disk-shaped monolayer cultures by application of Bay K 8644 (0.4 µmol/L) and recorded during 30-min periods using the electrode layout shown in Fig. 2 A. In 11 out of 12 different preparations, the activity was of focal origin, as illustrated by the isochrone maps shown in Fig. 2, B and C. In one experiment, spontaneous activity was due to reentry of the action potential.

View larger version (44K): [in this window] [in a new window]   FIGURE 2  Origin of focal activity in disk-shaped monolayer cultures of neonatal rat cardiomyocytes. (A) Layout of the recording electrodes (hexagonal lattice, marked with dots). The circle denotes the culture area. (B) Isochrone map of a focal activation originating within the hexagonal mapping area. The grayscale map was constructed by linear interpolation and the concentric circles by surface fitting (interval: 2 ms); the focus is marked by the solid dot. (C) Map of an activation originating beyond the hexagonal limit (arrow). (D) Histogram of focus occurrence density in the regions labeled a–f in the schematic (n = 11 different cultures).

In the analysis presented in Fig. 2 D, the culture area was divided into concentric regions in which the density of focus occurrences was calculated (i.e., the relative number of occurrences divided by the area of the corresponding region). This analysis, covering >96,000 focal activations in 11 different cultures, revealed a largest density of focus occurrences in the paracentral regions (b and c) and in the periphery of the preparations (f). These two maxima were, however, not significantly different.

Spatiotemporal analysis of beating variability in disk-shaped cultures
In a former study, we showed that BRV in the monolayer cultures exhibits a power law behavior with fractal features in the temporal domain, similar to HRV in vivo (31). We extended this study with the analysis of beating variability in the spatial domain by quantifying changes in focus position on a beat-to-beat basis.

Because the precise extrapolation of focus position beyond the area mapped by the electrodes was not possible, we conducted this analysis on 5 out of 11 experiments in which the focus remained in the hexagonal mapping area during the entire duration of the recording.

This analysis is illustrated in Fig. 3 for one experiment in a disk-shaped culture (Fig. 3 A) exhibiting focal activity (isochrone map and focus identification in Fig. 3 B). Instantaneous BR (Fig. 3 C) exhibited an aperiodic pattern of variability, as shown by the absence of peak in the PSD of BRV (Fig. 3 D). The linear aspect of the spectrum on the double logarithmic scale reveals the power law behavior of BRV. The power law exponent was –1.40, corresponding to the slope of the regression line. In Fig. 3 E, the beat-to-beat movement of the focus is illustrated by the trajectory interconnecting successive focus positions. The movement of the focus was irregular, exhibiting not only slow drifts and excursions but also more abrupt position changes. We computed the maximal focus excursion, defined as the largest distance between all foci identified during the 30-min recording. This maximal excursion, i.e., the largest focus displacement, was 2.78 mm and was thus orders of magnitude larger than the accuracy of focal point localization.

View larger version (48K): [in this window] [in a new window]   FIGURE 3  Spatiotemporal analysis of beating variability in a disk-shaped monolayer culture of neonatal rat ventricular myocytes. (A) Photograph of the preparation. The recording sites are marked with dots. (B) Isochrone map of an activation. The grayscale map was constructed by linear interpolation and the concentric circles by surface fitting (interval: 4 ms); the focus is marked by the solid dot. (C) BR as a function of time. Note the aperiodic variations of BR. (D) PSD of BRV. The spectrum was characterized by a power law (regression line in the double logarithmic plot) with an exponent of –1.40. (E) Left: Map illustrating the continuously varying position of the focus (dots) during the 30-min recording. The hexagon denotes the mapping area. Right: Magnification of the object formed by the ensemble of focus positions (shaded symbols) and the trajectory of the focus (interconnecting line). (F) Focus coordinates (axes shown in E) as a function of time, illustrating the irregular movement of the focus. The recording was segmented into 20-s bins. The solid line marks the median value for each bin, and the shaded band represents the range (minimum to maximum) for the corresponding bins. (G) Residence time distribution of the focus, i.e., distribution of the number of successive times a focus remained within a given square domain (grid illustrated in E) before moving to another. The residence time distribution followed a power law function with an exponent of –1.26.

We also described the dynamics of focus movement by evaluating during how many successive beats the focus remains at or close to a given location, thus giving indications about how frequently the focus moves to a remote location. For this purpose, we segmented the x-y plane into a grid with a given mesh size (0.49 x 0.49 mm in Fig. 3 E) and counted the number of successive times the focus remained within a given grid segment before moving to another. The histogram in Fig. 3 G illustrates the normalized distribution of this "residence time" (hereafter termed "residence time distribution"). The residence time distribution also followed a power law function with an exponent of –1.26. We repeated this analysis using different mesh sizes and observed that the distribution always followed a power law behavior for mesh sizes 0.3–1.4 mm (not shown).

Altogether, these analyses reveal that beating variability displays an aperiodic scale-independent fractal behavior not only in the temporal domain (i.e., concerning BRV) but in the spatial domain as well (i.e., concerning the origin of activation). Individual fractal measures and surface fit parameters characterizing the uniformity of conduction and accuracy of focus determination are presented in Table 1 for the five preparations investigated. In all five preparations, these analyses yielded similar results. In the temporal domain, BRV was characterized by a power law behavior with an average exponent of –1.38 ± 0.14 (consistent with our previous study (31)). In the spatial domain, the fractal box dimension of focus trajectory was 1.42 ± 0.22. The average power law exponents of the spectra of focus coordinates were –1.40 ± 0.16 and –1.35 ± 0.15 for x and y, respectively, without statistical difference between them (consistent with the symmetry and isotropy of the preparations). The average power law exponent of the residence time distribution (mesh size: 0.49 mm) was –1.69 ± 0.41. The exponents of the power law behavior of BRV of focus coordinates and of the residence time distribution were all in a similar range, and no correlation was found between them.

Pharmacological investigation of depolarizing mechanisms underlying the spontaneous electrical activity in cardiomyocyte cultures
In the SAN, it is well known that the L-type Ca²⁺ current I_Ca,L, the Na⁺/Ca²⁺ exchange current I_NCX consecutive to spontaneous Ca²⁺ release from the SR, and the hyperpolarization-activated cation current I_f are the major contributors to phase 4 depolarization (46–48). However, in cultures of neonatal rat ventricular myocytes, the ionic mechanisms underlying spontaneous depolarization and generation of electrical activity have not been fully characterized.

Here we used a pharmacological approach to evaluate the contribution of these depolarizing currents to spontaneous activity in the cardiomyocyte cultures. I_Ca,L was blocked or increased using nifedipine or Bay K 8644, correspondingly. Ryanodine and ZD7288 were used to block Ca²⁺ release from the SR (I_rel) and I_f, respectively. In all these experiments, the activity was always of focal origin.

Although 38% of preparations (16 out of 42) exhibited spontaneous activity (BR; 2.46 ± 1.04 Hz), most preparations were quiescent (62%). As shown in Fig. 4 A, nifedipine (5 µmol/L) silenced the beating of cultures that were active. Conversely, Bay K 8644 (0.4 µmol/L) markedly increased BR. Moreover, the application of Bay K 8644 initiated spontaneous activity in all the cultures that were initially quiescent.

View larger version (9K): [in this window] [in a new window]   FIGURE 4  Effects of pharmacological interventions targeting ion transport on the rate of spontaneous focal activity in monolayer cultures of neonatal rat ventricular cells. (A) Whereas nifedipine silenced the preparations (n = 6), Bay K 8644 strongly enhanced spontaneous activity (n = 12). Asterisks indicate p < 0.05 versus control. (B) In the presence of Bay K 8644 (0.4 µmol/L), ryanodine (10 µmol/L) depressed spontaneous activity (n = 6), whereas at a lower concentration (0.1 µmol/L), it slightly accelerated BR (n = 6). ZD7288 slightly decreased BR (n = 6). Asterisks indicate p < 0.05 versus control. (C) Corresponding interventions simulated using the WJvG model by scaling the corresponding conductances/permeabilities. Interruptions of the curves indicate conditions under which spontaneous action potentials ceased or the cell failed to repolarize normally.

These experiments were paralleled with computer simulations using the WJvG cell model by scaling the corresponding conductances or permeabilities. As shown in Fig. 4 C, the simulations resulted in qualitatively similar changes in BR. However, the dependence of BR on I_Ca,L was smaller in comparison with the changes observed in the cultures. BR exhibited a monotonic dependence on the total permeability of RyRs, and a reduction of I_rel always resulted in a decrease of BR. Finally, BR decreased to a similar extent when I_f conductance was reduced.

In the cell model, we also investigated the effects of changes in the maximal conductance of the delayed rectifier K⁺ current I_K (Fig. 4 C). The sensibility of BR on I_K density was prominent, in particular for decreased I_K conductance, which led to a marked reduction of BR.

Theoretical investigation of BRV induced by transmembrane current fluctuations
Voltage-gated currents are mediated by conformational changes of ion channels (e.g., openings and closings) in response to membrane potential changes. The ensemble current carried by channels of a given type is characterized by an average response to a given modification of membrane potential (e.g., a given voltage-clamp protocol or the action potential itself). However, the behavior of individual channels of the same type is different, because their conformational changes occur in a stochastic manner (50). This stochastic behavior will lead to beat-to-beat differences between the actual current and the average response. These deviations will vary from beat-to-beat and will also be different for distinct cells during their excitation by the same action potential.

To evaluate the effects of such fluctuations in transmembrane current on beating variability, we introduced into the model a transmembrane fluctuation current, I_fluct, mimicking the beat-to-beat differences caused by this stochastic open-close behavior of channels (see Materials and Methods).

Fig. 5 A illustrates the behavior of a single model cell in which an I_fluct with a standard deviation of 0.162 pA/pF was introduced. I_fluct resulted in BRV exhibiting a flat power spectrum (absence of power law behavior). The standard deviation of BRV was 0.037 Hz. As shown in Fig. 5 B, the standard deviation of BRV was proportional to the standard deviation of I_fluct. Excessive I_fluct (SD > 2 pA/pF) resulted in large and unrealistic oscillations of membrane potential and a disruption of this linear relationship.

View larger version (27K): [in this window] [in a new window]   FIGURE 5  Beating variability induced by stochastic fluctuations of transmembrane current in a single model cell (A and B) and in a network of pacemaker cells (C). (A) Action potentials, BRV, and power spectrum of BRV in the model cell after introduction of fluctuations (Ifluct) with a standard deviation of 0.162 pA/pF. The shaded line indicates BR in the absence of Ifluct. The spectrum was flat. (B) Standard deviation of interbeat intervals as a function of the standard deviation of Ifluct. The relationship is linear for Ifluct with standard deviation <2 pA/pF (slope of the regression line: 0.99). Larger Ifluct resulted in excessive perturbations of membrane potential (shaded area). The arrow corresponds to Ifluct used in A and C. (C) Map of focus occurrence, focus coordinates, BRV, PSD of BRV, and residence time distribution in the network. The spectrum exhibited an ascending behavior (positive power law exponent of +0.71). The residence time distribution was characterized by an exponential function (note the linear abscissa and the logarithmic ordinate).

Theoretical investigation of beating variability induced by stochastic open-close kinetics of RyRs
To evaluate the patterns of beating variability caused by variations in Na⁺/Ca²⁺ exchange current due to the stochastic nature of Ca²⁺ release events from the SR (51), the deterministic formulation of the average behavior of RyRs was replaced by a stochastic two-state gating model of individual RyR channels (see Materials and Methods).

In a single cell, as shown in Fig. 6 A, stochastic RyR gating led to fluctuations in the Na⁺/Ca²⁺ exchange current (I_NCX) and to corresponding variations in the rate of phase 4 depolarization (dV/dt). The dV/dt variations resulted in BRV (BR: 2.521 ± 0.011 Hz). However, the power spectrum of BRV was flat, indicating the absence of a power law behavior of BRV.

View larger version (28K): [in this window] [in a new window]   FIGURE 6  Beating variability induced by stochastic open-close kinetics of RyRs in a single model cell (A) and in a network of pacemaker cells (B; same network as in Fig. 5 C). (A) Left: Action potentials, Ca2+ release current (Irel), Na+/Ca2+ exchange current (INCX), and rate of rise of membrane potential during phase 4 depolarization (dV/dt). Variations in Irel caused by stochastic release events resulted in variations in INCX and modulated diastolic dV/dt (arrows). Right: BRV and its PSD in the single cell. The spectrum was flat (no power law). (B) BRV and its PSD in the cell network. The spectrum exhibited a power law with a positive exponent (+0.28, regression line, shaded).

Theoretical investigation of beating variability induced by ion channel turnover
In living cells, membrane proteins undergo an incessant turnover: synthesized membrane proteins are stored in transport vesicles and incorporated into the membrane. In parallel, membrane proteins are continuously internalized. This trafficking concerns ion channels as well (45). Because of the discrete nature of incorporation and internalization, the number of channels in a cell varies with time, and the average number of channels depends on the quantitative relation between the rates of synthesis/incorporation and recycling/degradation. Furthermore, distinct cells, being "on average identical", may have different densities of channels at a given instant. The variation of ion current densities in space and in time might lead to variations in the rate of a spontaneously active network and in focus position.

The results presented in Fig. 4 indicate that I_Ca,L density plays a significant role in determining BR in both the cultures and the WJvG model. Previous studies indicate that BR of pacemaker cells also depends on the density of the delayed rectifier K⁺ current (46), a finding corroborated by the simulations shown in Fig. 4 C. Therefore, if channel turnover leads to beating variability, this variability should be principally caused by the turnover of these two types of channels.

We first investigated the dynamics of ion channel turnover and its influence on beating variability in a single model cell. We formulated incorporation and recycling of ion channels as independent Poisson processes, based on the average number of channels of these two types and their corresponding half-lives reported in the literature (see Materials and Methods). Fig. 7 illustrates the variation of the number of channels in a cell and its consequence on BRV.

View larger version (35K): [in this window] [in a new window]   FIGURE 7  Model of ion channel turnover and effects of the turnover of delayed rectifier K+ channels and L-type Ca2+ channels on BR in a single cell. (A) Left: Variation of the number of delayed rectifier K+ channels, with a channel half-life of 4 h and an average number of 1000 channels. The turnover process is illustrated for five different 24-h periods. The number of channels followed a Poisson distribution (shown on the right). Right: Power spectrum of the variation of the number of channels. The spectrum followed a power law with an exponent of –2 and flattened for frequencies <10–5 Hz. (B) Turnover of L-type Ca2+ channels, with a channel half-life of 16 h and an average number of 10,000 channels. Same layout as in A. (C) Left: BRV in a single cell induced by the variations highlighted in A and B. Because the relative variation of the number of K+ channels SD/mean = 3.16%) was larger than that of the number of L-type Ca2+ channels (1%), BR was principally determined by the number of K+ channels. Right: The power spectrum of BRV followed a power law behavior. Over the range 10–5–1 Hz, the exponent was –1.97 (light regression line). However, the spectrum was less steep for frequencies >10–3 Hz (exponent of –1.78 in the range 10–3–1 Hz, darker regression line).

The contribution of the turnover of these two channel types to beating variability was then evaluated in a 2-D network of 156 WJvG cells (see Fig. 8, same network as in Fig. 5 C). Introduction of channel turnover led to BRV around the nominal frequency of the WJvG cell (Fig. 8 A). In the frequency domain, these variations were characterized by a power law behavior (Fig. 8 B) with an exponent of –1.53. In the spatial domain, activity was of varying focal origin, as illustrated in Fig. 8, C and D. The residence time distribution (Fig. 8 E) followed a power law function with an exponent of –1.69. The beat-to-beat movement of the focus (x and y coordinates as a function of time) is illustrated in Fig. 8 F. Fourier analysis of these coordinates revealed an absence of periodic behavior and the PSD of the coordinates followed a power law behavior with exponents of –1.21 and –1.36 for the x and y coordinates, respectively.

View larger version (27K): [in this window] [in a new window]   FIGURE 8  Spatiotemporal characteristics of beating variability in a network of WJvG model cells with turnover of L-type Ca2+ and delayed rectifier K+ channels. (A) BR of the network (layout in C and D). The dotted line indicates BR in the absence of channel turnover. (B) PSD of BRV. The exponent of the power law was –1.53. (C) Isochrone map of an activation. (D) Map illustrating the relative occurrence of the focus in different cells of the network. The line connects successive focus positions. (E) Residence time distribution. The distribution followed a power law (exponent of –1.69). (F) Focus coordinates (origin at the center of the network), illustrating the irregular focus movement.

Effects of network dimensionality and size on the power law exponent ß
The finding that the exponent ß is different in single cells versus cellular networks indicates that the power law characteristics of BRV are influenced by intercellular interactions. These interactions are likely to be determined by both the dimensionality (i.e., one, two, or three dimensions) and the size of the network. Therefore, we investigated the effects of network dimensionality and size on the exponent ß, as presented in Fig. 9.

View larger version (12K): [in this window] [in a new window]   FIGURE 9  Effects of network dimensionality and size on the power law exponent. (A) Power law exponent as a function of dimension in symmetrical networks of approximately the same size. (B) Power law exponent as a function of dimension in networks of approximately the same diameter (seven cells). Asterisks denote significant differences (n = 5).

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Under our experimental conditions, spontaneous activity in cultures of densely seeded neonatal rat ventricular myocytes originates from a mobile focus and spreads radially across the entire tissue. The generation of this rhythm is a highly dynamic process and is based on ionic mechanisms also involved in the generation of the sinoatrial rhythm. The variability of both BR and of focus position is aperiodic, and the respective power spectra exhibit power law behavior. The numerical simulations suggest that these characteristics cannot be accounted for by "white noise" in the biological system, e.g., by fluctuations of transmembrane currents or stochastic release of calcium from intracellular stores. Rather, they necessitate the presence of mechanisms exhibiting correlations occurring over a broad range of temporal scales, e.g., as investigated in this study, the number of channels in the cell membrane.

Mechanisms of beating variability
In a previous study, Wilders et al. modified their original SAN cell model (which we used in the theoretical part of this work) by incorporating a detailed formulation of the stochastic gating of all individual ion channels (44). Their model mimicked the irregular beating observed experimentally in isolated SAN cells. Stochastic channel gating led to uncorrelated interbeat intervals. A similar study was conducted by Guevara and Lewis (52), supporting the notion that the stochastic behavior of a population of channels could contribute to irregular beating. The spatiotemporal patterns of beating variability in cell networks were, however, not addressed in these studies.

To address the question of whether modulation of BR caused by stochastic ion channel gating may lead to the power law behavior of BRV, we incorporated into the model a transmembrane fluctuation current with temporal correlations that did not extend beyond the intrinsic beating period of the cell. This resulted in BRV with a flat power spectrum and thus an absence of correlations in the interbeat interval series, equivalent to the result of Wilders and Jongsma (44). We opted for this approach because long-term simulations incorporating the stochastic gating of every channel in every cell of the network would have been computationally prohibitive.

The application of I_fluct in the cell network resulted in a power law exponent of BRV of the network that became clearly positive, indicating that the intercellular interactions in the network act as a "high-pass filter", attenuating preferentially the slow variations in BR. A similar result was observed for different intensities of I_fluct and for I_fluct formulated as Gaussian white noise (not shown).

The fact that Ca²⁺ release from the SR (in the form of Ca²⁺ sparks) is strongly involved in the modulation of BR both in SAN cells (47,53) and in our cultures (Fig. 4) motivated us to investigate how the stochastic nature of Ca²⁺ release may contribute to BRV. Although we used a formulation of Ca²⁺ release that does not account for the subcellular features of the release system, our result nevertheless suggests that BRV induced by the stochastic opening of RyRs is characterized by an absence of temporal correlations. Accordingly, the spatiotemporal characteristics of beating variability were similar as those induced by fluctuations of transmembrane current.

Thus, the stochastic gating of sarcolemmal ion channels or of RyRs represent unlikely candidates to explain the characteristics of beating variability and the power law of BRV observed in the experiments. Other possible biophysical, biological, or molecular mechanisms underlying beating variability can then be envisioned from the fact that every cell is a dynamic structure in which a continuous self-rearrangement and numerous regulatory processes concerning all constitutive elements of the cell take place. Variations in cellular characteristics can concern all biochemical processes (concentrations of metabolites and signaling molecules; concentrations and activity of enzymes) and therefore the regulation of the biophysical properties of ion transporters (e.g., by phosphorylation or expression of different isoforms, circadian rhythms) (54). There is, however, no experimental approach that would permit a tight control of all these parameters at the cellular level and thus an experimental investigation of their involvement in beating variability. Furthermore, there is to date no mathematical model of a pacemaker cell that incorporates the details of biochemical reactions or molecular mechanisms in conjunction with the electrical function of the cell. Moreover, in electrophysiological models published to date (whether incorporating a deterministic or a stochastic approach to channel gating), the maximal conductances of ion currents, i.e., the number of channels, is considered to be time invariant.

To obtain insight into the electrophysiological consequences of one aspect of cellular self-rearrangement, we challenged this notion of invariance and developed an original formulation for the turnover of ion channels in which the number of channels results from the integration of isolated incorporation/removal events. The turnover of ion channels introduced longer time constants (several hours) into the WJvG model, leading to longer term fluctuations in the simulations. The temporal dynamics of the number of channels followed a time course characteristic of Brownian motion with long-term correlations and with a power law exponent of –2, which, at the single cell level, resulted in a power law behavior of BRV with an exponent in this range. In the 2-D 156-cell network, the high-pass filtering properties of the network then reduced the spectrum slope to values close to –1.5, similar to those observed experimentally. This result suggests that it is necessary to consider phenomena characterized by long time constants and/or long-term correlations to explain the power law behavior of BRV. Nevertheless, the amplitudes of the variations of BR were notably smaller than those observed in the cultures, indicating that ion channel turnover alone cannot account for the full intensity of variability in the cultures.

Intuitively, one could explain the high-pass filtering effect by the fact that the network adapts itself to a change in intrinsic frequency of a particular cell. Because this adaptation takes a certain time (possibly depending on the position of the cell within the network), the network can more easily accommodate slower than faster variations of the intrinsic beating rate of individual cells. This easier accommodation to slower changes may then contribute to a stronger attenuation of low-frequency variations and a weaker attenuation of high-frequency components.

In our experiments, the variability of focus position was characterized by a power law behavior similar to that of BRV, a characteristic that was closely reproduced by the model incorporating the turnover of ion channels. This reflects the fact that action potentials are not initiated at random positions but focus position is related to the previous history of spontaneous activity. This memory of the previous history of focus position can be inferred from the residence time distributions. If focus positions were completely uncorrelated from beat-to-beat (corresponding, e.g., to blind throws of an object onto a grid with a probability, a, to hit a given grid segment), then the probability of having N successive occurrences at the same position would follow an exponential distribution of the type a^–N. The power law behavior of the residence time distributions indicates that, compared to an exponential distribution, the probability of having N successive occurrences at the same position is relatively higher for large N and relatively lower for small N, respectively. In other words, the longer the focus remains at or near a given location, the less likely it is to move to another location. Conversely, if the focus just moved to another location, it is relatively less probable that it will remain there. The presence of correlations extending over a broad range of temporal scales, i.e., an influence of the past history on the current dynamics, is also reflected by the behavior of focus position coordinates as a function of time.

Based on these considerations, we thus propose that the power law behavior of beating variability must be, at the level intrinsic to cardiac tissue, essentially caused by variations of cellular characteristics exhibiting long-term correlations. Ion channel turnover can represent a stochastic mechanism contributing to such variations of cellular characteristics. However, the largest part of variability probably emerges from variations (deterministic or stochastic) in biochemical or molecular processes involving the concentrations of enzymes, metabolites, and second messengers. For example, it has recently been shown that levels of cyclic adenosine 3',5'-phosphate (cAMP), protein kinase A, phospholamban phosphorylation, and L-type Ca²⁺ channel phosphorylation are dynamically involved in the control of BR in SAN cells via modulation of Ca²⁺ cycling (55). The investigation of how these biochemical aspects are involved in beating variability was, however, beyond the extent of our study.

Relationship between spatial and temporal dynamics
In the SAN, the origin of the spontaneous action potential is not static and changes according to the autonomic nervous tone. It has been shown that the site of initiation shifts in response to vagal stimulation (56,57) and that shifts in the site of origin correlate with changes in heart rate (58,59). However, our cultures are devoid of autonomous influences, and the correlation between focus shifts and changes in BR is negligible. This absence of correlation was also observed in the simulations, irrespective of the investigated mechanism of rate variability. Therefore, it is possible, as proposed previously, that the correlation present in the SAN is due to inhomogeneous vagal innervation and variable cell sensitivity to acetylcholine (60,61).

The SAN is often described as a population of electrically coupled pacemakers having different intrinsic frequencies. Their coordinated rhythm depends on the process of "mutual entrainment", which leads to a stable periodic behavior of the pacemaker cell network (62,63). In this situation, the BR of the SAN is the result of a "democratic" process underlying the coordinated firing of the SAN cells, with each cell contributing to the overall mutual entrainment (64). Moreover, the BR of the network and the focus of origin are not necessarily that of the intrinsically fastest pacemaker cell (65). The small or absent correlation between changes in focus location and sudden changes in BR in both the cultures and the mathematical model are thus not in discrepancy with this theory of synchronization.

Ionic m