Synchronization of Pacemaking in the Sinoatrial Node: A Mathematical Modeling Study

Computational studies using mathematical models of the sinoatrial node (SAN) cardiac action potential (AP) have provided important insight into the fundamental nature of cell excitability, cardiac arrhythmias, and potential therapies. While the impact of ion channel dynamics on SAN pacemaking has been studied, the governing dynamics responsible for regulating spatial and temporal control of SAN synchrony remain elusive. Here, we attempt to develop methods to explore cohesion in a network of coupled spontaneously active SAN cells. We present the updated version of a previously published graphical user interface LongQt: a cross-platform, threaded application for advanced cardiac electrophysiology studies that does not require advanced programming skills. We incorporated additional features to the existing LongQt platform that allows users to (1) specify heterogeneous gap junction conductivity across a multicellular grid, and (2) set heterogeneous ion channel conductance across a multicellular grid. We developed two methods of characterizing the synchrony of SAN tissue based on alignment of activation in time and similarity of voltage peaks among clusters of functionally related cells. In pairs and two-dimensional grids of coupled cells, we observed a range of conductivities (0.00014–0.00033 1/Ω-cm) in which the tissue was more susceptible to developing asynchronous spontaneous pro-arrhythmic behavior (e.g., spiral wave formation). We performed parameter sensitivity analysis to determine the relative impact of ion channel conductances on cycle length (CL), diastolic and peak voltage, and synchrony measurements in isolated and coupled cell pairs. We also defined measurements of evaluating synchrony based on peak AP voltage and the rate of wave propagation. Cell-to-cell coupling had a non-linear effect on the relationship between ion channel conductances, AP properties, and synchrony measurements. Our simulations demonstrate that conductivity plays an important role in regulating synchronous firing of heterogeneous SAN tissue, and demonstrate how to evaluate the role of coupling and ion channel conductance in pairs and grids of SAN cells. We anticipate that the approach outlined here will facilitate identification of key cell- and tissue-level factors responsible for cardiac disease.

Introduction

The sinoatrial node (SAN) generates the electrical impulse that coordinates mechanical contraction of the heart [1, 2]. Proper SAN function is an essential component for normal pacemaking at baseline and heart rate variation in response to external regulators such as exercise or stress [3, 4]. SAN dysfunction is common in a wide variety of cardiac diseases, and is characterized by sinus bradycardia, sinus pause, and/or inappropriate heart rate responses to exercise and stress [5]. Regulation of SAN activity has great therapeutic potential for a rapidly aging population where SAN disease affects 1 in 600 heart patients over the age of 65 [6]. The only effective treatment for patients with chronic symptomatic sinus node dysfunction is pacemaker implantation [7].

SAN cells demonstrate spontaneous action potential (AP) activity and exhibit a wide variety of dynamic phenomena similar to other coupled oscillators, including collective synchronization [8]. One challenge for studying synchronization of cardiac pacemaking activity is the multiscalar and heterogeneous nature of the sinus node. Pacemaking is governed by a delicate source-sink relationship between the SAN and surrounding atria defined by the need for a relatively small structure (SAN) to excite a much larger tissue mass (surrounding atria) [9, 10]. This source-sink relationship is altered in disease due to increased fibrosis and/or cell loss leading to a shift of the primary pacemaker site, emergent behavior of ectopic foci, or otherwise reduced capacity for SAN pacemaking [11–13]. There is a critical need to expand knowledge regarding regulation of membrane ion channels in the SAN, as well as to further develop quantitative tools to assess the sensitivity of the SAN to changes in coupling and ion channel regulation [14].

Mathematical modeling has been used to investigate and advance our understanding of cardiac electrophysiology, arrhythmia mechanisms, and potential therapies [15]. Models have been particularly helpful in elucidating the ionic basis of SAN activity and cardiac pacemaking [16, 17]. For example, SAN cell models have furthered our understanding of the relative importance of coupling between Ca²⁺ cycling and membrane ion channels in automaticity [18] and the genetic basis of human SAN disease [19]. At the same time, multicellular models of coupled SAN cells have demonstrated dynamic changes in the location of the primary pacemaker site in response to β-adrenergic stimulation [20]. Other studies have coupled heterogeneous cell types in multicellular preparations to examine the influence of inexcitable cells (e.g., fibroblasts) on pacemaking [21, 22].

Although mathematical modeling has undoubtedly advanced our understanding of SAN function, there remain significant barriers to more widespread use of mathematical modeling and simulation in the field. To reduce these barriers a cross-platform user interface called LongQt has been developed for advanced cardiac electrophysiology and arrhythmia simulations [23]. Here, we present an extended user interface for LongQt, which adds support for performing advanced multicellular simulations. With this added utility, the influence of perturbations in SAN cell electrophysiology (ion channel conductance and gap junction conductivity) on synchronization of coupled heterogeneous SAN cells was investigated. Two complimentary measures were defined to quantify the level of synchronization between coupled cells: synchrony factor and peak transmembrane potential (V_m,peak) similarity. Parameter sensitivity analysis was performed in single cells and in coupled cell pairs to determine the impact of ion channel and gap junction properties on cycle length (CL), peak and diastolic V_m, and measures of synchrony. Finally, we performed two-dimensional simulations in a network of 7 × 7 heterogeneous SAN cells to examine the influence of cell properties and gap junction conductivity on pacemaking. These studies generate a number of interesting findings, including: (1) conductivity caused small but potentially important differences in the relative impact of ion channel conductances on AP properties in coupled cells; (2) a specific coupling range promoted emergent asynchronous behavior; and (3) quantitative measurements were defined to evaluate synchrony based on peak AP voltage and the rate of propagation within a group of coupled firing cells. While our findings highlight the difficulty of relating events at the single cell level to an emergent behavior like pacemaking, they also point to more robust methods for understanding the ionic basis of cardiac pacemaking.

Materials and Methods

Ion channel kinetics were simulated using an existing well-validated model of the rabbit SAN cell implemented in LongQt simulation software [23, 24] (Figure 1A). Briefly, the LongQt simulation software has three main user interfaces: the grid editor, the main user interface, and the grapher. The simulations performed for this study were set up using the grid editor, which allows the user to select tissue geometry and gap junction conductivity for a set of simulations. The files generated by the grid editor can be selected to run in the main user interface, which also allows the user to select the cell model and measure properties of the simulation. Simulation results generated at the end of the simulation can be visualized by the grapher interface.

FIGURE 1

Figure 1. (A) LongQT simulation software uses an object-oriented design where a “Cell” class describes the workings of an individual cell, i.e., Rabbit SAN cell vs. Human Ventricular cell. Cell classes then are used by “Protocol” classes to define what type of simulation is being run on the cell. This could be a simulation with constant stimulus or periodic stimulus. Simulations of fiber and grid tissue require additional objects that go between Protocols and Cells to define properties of the tissue like its geometry and conductivity. Finally there is a ”Measure” class which tracks values such as Peak Voltage or Ca²⁺ minimum. (B) Schematic diagram of the rabbit SAN cell model, which includes mathematical representations of ion currents important for generating the sinoatrial node cell action potential: L-type and T-type Ca²⁺ currents (I_CaL and I_CaT), Rapid and slow components of the delayed rectifier K⁺ currents (I_Kr and I_Ks), 4AP-sensitive transient outward and sustained K⁺ currents (I_to and I_Sus), acetylcholine-sensitive K⁺ current (I_Kach), Hyperpolarization-activated funny current (I_h), Na⁺/Ca²⁺ exchanger (I_NCX), Na⁺/K⁺ pump (I_NaK), and Ca²⁺ signaling uptake (j_up) into the network sarcoplasmic reticulum (NSR), transfer (j_tr) into the junctional SR (JSR), and release (j_rel) into the subspace.

Multicellular Simulations

Multicellular simulations were performed in either a cell pair or 7 × 7 grid. The two-dimensional cable equation was solved using the Peaceman-Rachford alternating direction implicit method. The level of conductivity between cells was perturbed about a default value of 0.33 1/Ω-cm. A heterogeneous population of SAN APs were created by varying eight ion channel conductances (I_Ca,L, I_Catt, I_h, I_Kr, I_Ks, I_NCX, I_NaK, I_to) lognormally, with a mean of 1 and a standard deviation of 0.2. Spontaneous APs were simulated for 50–100 s until steady state was reached. AP properties such as CL, V_m,_peak, and maximum diastolic potential (MDP) were measured using LongQt. All other analysis was performed with scripts written in Python version 3, which are available on Github.

Synchrony Measurements

In order to develop measurements of synchrony in multicellular simulations, we organized APs from individual cells in the grid into “activation clusters,” which represent a group of neighboring cells whose activation could be considered related both in time and space. To organize cells into clusters, we ordered them sequentially according to their respective activation times. The sequence was then processed in order and a cell was added to a cluster if its position was within three cells of any cell already in the cluster. For higher gap junction conductivities, we increased the spatial window to five cells to account for increased communication between cells. If an activated cell was not spatially close enough to any existing group of firing cells, then it was marked as the focus of a separate and distinct cluster. Any cluster would be considered complete when one of its constituent cells fired again. This allowed for characterization of multiple clusters simultaneously within the same grid (common in lower conductivity grids).

For example, given an ordered sequence of cell activation times (Cell_1,1, Cell_1,2, Cell_5,5, Cell_1,1, Cell_1,3), where Cell_{i, j} is located in the i^th row and j^th column of a two-dimensional grid, the clustering algorithm would select Cell_1,1 as the beginning of a new cluster, C₁, as there are no other existing clusters. The second cell in the sequence, Cell_1,2, would then be added to C₁ as it is within three cells of Cell_1,1. Cell_5,5, however, is too far away from C₁ and so would be marked as the beginning of a new cluster, C₂. Since Cell_1,1 is already assigned to C₁, its appearance a second time in the sequence causes the algorithm to initialize a new cluster C₃. Finally, the last element in the sequence, Cell_1,3, would be added to C₃.

V_m,peak similarity is calculated as the inverse of the standard deviation of V_m,peak in each cluster. The synchrony factor is calculated as the inverse of the slowest propagation time between the closest cells in a cluster. These two measurements are then weighted by the size of the cluster in order to account for the number of oscillators in the network.

Software and Hardware

LongQt simulation software utilizes the Qt application framework (version 5.6 or later found at https://www.qt.io) and may be compiled to run on Mac (OS X 10.10 or later), Windows (version 7 or later) or Linux systems. Python bindings for LongQt are available for more extensive simulation use. Compiled versions of LongQt are available as downloadable executable files under the “Research” section of the Hund lab website¹, and are accessible through Github². LongQt incorporates C++ code for the Kurata SAN cell model (Figure 1B). Differential equations for the simulated model are solved in LongQt using the forward Euler approach, with a maximum timestep of 0.05 ms and a minimum timestep of 0.005 ms. A subset of simulations were run using Ohio Supercomputer Center resources [25]. Data supporting conclusions of this manuscript are available upon request to the corresponding author.

Results

Effect of Coupling on Parameter Sensitivity in Coupled Pairs of Sinoatrial Node Cells

Parameter sensitivity analysis has been performed mostly on models of the single cell to elucidate mechanisms underlying cardiac AP generation [26–30]. Using the extended LongQt platform, we sought to address the extent to which coupling influences the dynamics of synchronous firing in a network of SAN oscillators with heterogeneous ion channel activity. We first coupled a simulated wild-type (WT) to a variant SAN cell (scaling factors defined as follows: I_Ca,L = 3.19472, I_Catt = 2.59237, I_h = 2.64054, I_Kr = 2.3018, I_Ks = 2.93799, I_NCX = 2.79239, I_NKA = 2.45686, I_to = 2.06311) and observed spontaneous AP properties for coupled and uncoupled pairs (Figures 2, 3). As expected, despite very different AP properties of the individual cells, a normal degree of coupling eliminated differences between the two cells. Interestingly, the steady-state CL, DDR, and MDP values of the coupled cells were closer to that of the single cell with the shortest CL (in this case, the variant cell). This phenomenon is consistent with observations of shifts of the SAN pacemaker site in cardiac disease from the original area of excitation to areas with the earliest depolarization [31].

FIGURE 2

Figure 2. Simulated spontaneous SAN action potentials in one WT cell and one cell with lognormally perturbed ion channel factors. Simulations were run to steady state (50 s). (A) At a conductivity of 0.33 (1/Ω-cm), which is the normal conductivity between two WT cells, the cells synchronize both the voltages and the times at which they fire. (B) At a conductivity of 0 (1/Ω-cm) the cells are uncoupled and act the same as if they were run independently.

FIGURE 3

Figure 3. Measured steady state values from the last action potential in the simulation. The differences between the two cells are resolved at full conductivity to a value somewhere between the two cells. (A) Cycle length, (B) diastolic depolarization rate, (C) max diastolic potential, (D) max voltage.

To provide insight into the influence of coupling on spontaneous AP dynamics, we performed parameter sensitivity analysis on the single cell by generating 616 AP variants and performing linear regression on the dataset. We compared these results to a separate regression analysis on a dataset where the variant was coupled to a WT cell (Figure 4). For the most part, regression coefficients relating ion channel conductances to AP properties were similar for single and coupled cells. For example, perturbations in maximal conductances of the L-type Ca²⁺ current (I_Ca,L) and the transport rate of the Na⁺/K⁺ ATPase (I_NaK) have a large positive influence on V_m,peak in both single and coupled simulations, while perturbations in conductances of the T-type Ca²⁺ current (I_Catt), rapid delayed rectifier K⁺ current (I_Kr), and transient outward K⁺ current (I_to) are inversely related, i.e., an increase in I_Catt, I_Kr, or I_to decrease peak V_m. Despite the overall agreement, there are small but interesting differences between sensitivity of the single vs. coupled cell. First, while I_Ca,L has a positive effect on CL and V_m,peak in both single and coupled cells, its influence is diminished in the coupled cell. Likewise, coupling reduces the influence of I_Ca,L on DDR and MDP. In contrast, our simulations predict that coupling increases the effect of I_Kr, at least with respect to DDR and CL.

FIGURE 4

Figure 4. Partial least-squares regression analysis of ionic gating variables in the rabbit sinus node model. Regression coefficients showing how changes in model parameters affect membrane dynamics in two coupled simulated SAN cells. 616 simulations were performed at a normal coupling strength (0.33) between one WT and one variant cell with random ion channel factors (L-type Ca²⁺ current I_Ca,L, T-Type Ca²⁺ current I_Ca,T, hyperpolarization-activated current I_H, rapidly activating delayed rectifying K⁺ current I_Kr, slowly activating delayed rectifying K⁺ current I_Ks, Na⁺/Ca²⁺ exchanger I_NCX, Na⁺/K⁺ ATP-ase I_NaK, transient outward K⁺ current I_to) perturbed over a lognormal distribution, with a mean of 1.0 and a standard deviation of 0.2. Parameter sensitivities of ion channel conductance parameters affect (A–B) peak membrane voltage, (C–D) cycle length, (E–F) max diastolic potential (MDP), and (G–H) diastolic depolarization rate (DDR). The impact of I_Ca,L, I_Catt, I_to, I_Kr, and I_NaK contributed highly to membrane voltage dynamics including peak (I_Ca,L, I_NaK, I_to), cycle length (I_Catt), diastolic membrane voltage (I_Kr, I_NaK), and diastolic depolarization rate (I_NaK).

To provide additional insight into the influence of coupling strength on sensitivity analysis, we performed 6160 simulations (10 conductivities, 616 simulations of a WT cell coupled to a variant cell) over a range of conductivities (Figure 5). In many instances, the regression coefficients mapping ion channel conductances to AP properties were found to be independent of coupling strength, especially for V_m,peak and MDP. However, interesting exceptions to this behavior were observed for CL and DDR, where regression coefficients for specific ion channels were highly dependent on coupling strength (e.g., I_Ca,L, I_Ca,tt, I_h for CL and I_h and I_Kr for DDR) This series of simulations suggests that the relative importance of specific ion channels for cardiac pacemaking changes in subtle but important ways across a range of coupling values.

FIGURE 5

Figure 5. Partial least-squares regression analysis of ionic gating variables vs. cell conductivity. Two cells were paired where one was lognormally perturbed with a mean of 1 and a standard deviation of 0.2, while the other cell was WT. 616 simulations were run at each conductivity for the least-squares regression at each of 10 different conductivities. Regression coefficients are shown over a range of gap junction conductivities for: (A) Max voltage, (B) cycle length, (C) max diastolic potential, (D) diastolic depolarization rate.

We sought to explore coupling effects over a range of 25 different coupling strengths with a WT cell coupled to 20 different variants (totaling 500 different simulations of two coupled cells). We plotted the average steady-state values of the MDP, V_m,peak, and CL for both the coupled WT and variant cells at each coupling value. As expected, values for MDP, V_m,peak, and CL converge in the WT and variant cell as gap junction conductivity increases (Figure 6). Interestingly, CL appears to synchronize at lower conductivity values compared to other AP properties. Furthermore, a small window of interesting dynamics characterized by increased standard deviation values for AP properties was observed around 10^−3.8 1/Ω-cm.

FIGURE 6

Figure 6. Conductivity vs. cell properties for two coupled SAN cells where one cell was perturbed lognormally, while the other cell is a WT. Data points are the average of 20 simulations with the error bars corresponding to one standard deviation. (A) Max diastolic potential vs. Conductivity. Max diastolic potential synchronizes at higher conductivities and the means become equal around 10^−1.6 1/Ω-cm. Larger standard deviations are observed around 10⁻⁴ 1/Ω-cm. (B) Peak voltage vs. conductivity. Peak voltage equalizes at higher conductivities and the means become equal around 10⁻¹ 1/Ω-cm. Larger standard deviations are observed in the 10⁻⁴ and 10⁻³ 1/Ω-cm range. (C) Cycle Length vs. Conductivity. Cycle length synchronization happens at much lower conductivities than the others, with the means equalizing between 10⁻⁴ and 10⁻³ 1/Ω-cm. It is in this range that the cells are highly interactive but not fully able to synchronize, which is indicated by large standard deviations.

Effect of Coupling in Two-Dimensional Simulations of Heterogeneous SAN Cells

We hypothesized that the range of conductivities identified in Figure 6 with large standard deviations would promote asynchronous activity in two-dimensional simulations of heterogeneous SAN cells. We simulated 7 × 7 grids of variant SAN cells with homogeneous cell-to-cell coupling, for 25 different conductivities. At low coupling values most SAN cells oscillate without interacting with each other (Supplementary Videos 1, 2). The range of values indicated by arrows in Figure 6 was also the range in which spiral wave activity was sustained in the two-dimensional grid simulation (Supplementary Videos 3, 4). For higher degrees of cell-to-cell coupling, cells across the grid were fully synchronized (Supplementary Videos 5, 6). This set of simulations suggests that coupling can promote an arrhythmogenic substrate, even in a small group of pacemaker cells.

Quantifying Synchrony in Two-Dimensional Simulations of Heterogeneous SAN Cells

We sought to quantify the level of synchrony in a two-dimensional grid in order to quantitatively distinguish between spiral wave formation (Supplementary Videos 3, 4), completely asynchronous activity (Supplementary Videos 1, 2) and synchronous firing (Supplementary Videos 5, 6). Plotting beat-to-beat CL for all 49 cells in a 7 × 7 grid for three different gap junction conductivities demonstrates a wide range of steady-state CLs when the propagating wave is random (Figure 7A), large deviations when spiral waves are formed (Figure 7B, <20 s), and a uniform CL across the grid when the tissue is oscillating synchronously (Figure 7B, >20 s and Figure 7C). Plotting beat-to-beat CL at normal conductivities shows uniform CL across all cells in the simulation (Figure 7C).

FIGURE 7

Figure 7. Simulations were run on a 7 × 7 grid for 50 s. Each cell was randomly perturbed using a lognormal distribution with a mean of 1 and a standard deviation of 0.2. (A–C) Cycle length for 49 cells is shown over time. At low conductivities (A) cells do not synchronize, and as conductivity increases the cells begin to interact and find a common cycle length (B) and eventually fire at a common cycle length immediately (C). (D–F) Synchrony factor vs. Time. The synchrony factor measurement is the longest amount of time it takes the peak of the action potential to propagate from any cell to its neighbor. Low values indicate random and low synchrony factor waves (D). The larger peaks and valleys that appear in higher conductivities (E) correlate with spiral wave formation and multiple wave fronts in the grid. When cells interact and fire in cohesive synchronization (F), synchrony factor becomes higher and stabilizes. (G–I) Peak voltage similarity vs. Time. The peak voltage similarity measurement is the inverse of the standard deviation of the voltage values across the grid. Low conductivities indicate random noisy peak voltage similarity waves (G). As conductivity increases, peak voltage similarity increases and stabilizes initial noisy values (H). At higher conductivities, peak voltage similarity display no errant behavior (I).

The synchrony factor, which represents the inverse of the longest conduction time between two cells in the same cluster, approaches zero and demonstrates noise (0.5 amplitude trace) for chaotic asynchronous simulations where the cells are not interacting with each other (Figure 7D). When coupling is in a range that sustains spiral wave activation, the spiral waves can be visualized in the peaks and valleys of synchrony factor over time (Figure 7E, <20 s). As coupling increases to normal consistent propagation across the entire grid, the synchrony factor increases and maintains a steady value (Figure 7E, >20 s and Figure 7F). Synchrony factor values above 1 consistently represented fully synchronized grids, and below 0.5 consistently represented asynchronous random poorly coupled oscillations.

V_m,peak similarity, which represents the inverse of the standard deviation of peak voltages in one beat, approaches zero with high-amplitude fluctuations for chaotic asynchronous simulations (Figure 7G). As a simulation transitions from asynchronous (Figure 7H, <5 s) to an organized, complex activation (spiral wave, Figure 7H, >5 and <20 s) V_m,peak similarity rapidly reaches a single steady-state value after a brief period of low amplitude fluctuation. Synchronized activation produces a large steady-state value for V_m,peak similarity with a brief latency (Figure 7I). These results demonstrate the utility of quantifying synchrony measures to distinguish between random, spiral, and synchronous propagating waves sustained by coupling differences in heterogeneous SAN tissue. V_m,peak similarity values above 15 consistently represented fully synchronized grids, and below 10 consistently represented asynchronous random poorly coupled oscillations.

We performed parameter sensitivity analysis on coupled cells to determine the relative ion channel contributions to the synchrony factor and peak voltage similarity measurement (same 616 simulations as Figures 4A–D). At normal gap junction conductivity, synchrony factor is not dominated by any single ion channel conductance (relatively small regression coefficients for all conductances) with surprisingly little contribution from I_h (Figure 8A). However, I_Ca,L and I_NaK both had a large negative contribution and I_to a large positive contribution to V_m,peak similarity (Figure 8B). The relationships observed in Figure 8B also seemed to be an inverse of the contributions to peak voltage in previous simulations (Figure 4A). When we performed parameter sensitivity analysis over a range of conductivities, we observed that the relationship between each individual ion channel's contribution and synchrony factor was non-linear with respect to conductivity (Figure 4C). Interestingly, any shifts from a positive to negative contribution occurred in the range of 10⁻⁴ 1/Ω-cm which is the same range that we observed sustained spiral wave activity (Supplementary Videos 3, 4) and large standard deviations in CL (Figure 6C).

FIGURE 8

Figure 8. Partial least-squares regression analysis of synchrony factor and peak voltage similarity metrics. 616 simulations were run with two paired cells where one was lognormally perturbed with a mean of 1 and a standard deviation of 0.2, while the other cell was a WT. Ion channel coefficients did not have a significant impact on (A) synchrony factor. In contrast, parameter sensitivity analysis indicates that the impact of I_Ca,L, I_NaK, and I_to greatly impact the (B) peak voltage similarity metric. (C,D) Partial least-squares regression analysis of ionic gating variables vs. cell conductivity. At each conductivity simulations were run as in (A) and (B) and a total of 10 different conductivities were examined. (C) Synchrony factor is very non-linear with a range between 10⁻⁴ and 10⁻³ 1/Ω-cm in which all the ion channel coefficients are very close to 0. (D) Peak voltage similarity is mostly linear with an exception between 10⁻⁴ and 10⁻³ 1/Ω-cm where most of the ion channel coefficients pinch toward 0.

Discussion

In this study, we use mathematical modeling to explore the role of coupling on spontaneous AP dynamics and synchronization of pacemaking. Our simulations led to a number of important findings, including: (1) While parameter sensitivity analysis reveals a similar relationship between ion channel conductances and AP properties in single and coupled cells, our simulations predict small but potentially important differences, including complicated effects of coupling on the influence of I_Ca,L and I_Kr; (2) a specific coupling range in simulations promoted complex emergent behavior (including spiral wave activation) and at values higher than this coupling range cells fired together synchronously; (3) We define an approach for first defining groups of related cells (activation clusters) and then characterizing their synchrony (synchrony factor and peak voltage similarity), which facilitates quantification and visualization of synchronous behavior in a two-dimensional heterogeneous grid of SAN cells. Our studies are distinct from previous studies investigating coupling between spontaneously activating oscillators in that we employ an AP model that describes detailed ion channel kinetics. Another novel aspect of this set of studies is the introduction of updated LongQt simulation software to explore the impact of heterogeneous ion channel expression and gap junction conductivity in multicellular simulations. LongQt is cross-platform and available for download at hundlab.org, and may be useful in future exploration of conductivity in two-dimensional simulated tissue.

Previous studies have explored coupling inhomogeneity between simulated SAN cells and observed trends of synchronous firing and heterogeneous tissue becoming homogeneous through a democratic entrainment process at sufficient coupling values [20, 32, 33]. Our studies also support the theory of a democratic entrainment process both at the level of two coupled cells (Figure 6) where both cells adjusted their transmembrane dynamics to adjust to a new value that was distinct from firing alone. All grid simulations with uniform wave propagation began activation from a cluster of cells firing together rather than a single cellular driver (Supplementary Videos 1–6). This cluster size was different for different coupling values, indicating that the multicellular simulations demonstrated mutual entrainment of SAN cells.

The SAN is a small structure that is insulated from the rest of the right atrium, and employs a limited number of conduction pathways in order to activate the surrounding tissue [34, 35]. In the SAN, cells form groups with high degrees of coupling between cells in a group and much lower amounts of coupling at the border of groups [36, 37]. Conduction barriers due to fibrosis or structural remodeling may inhibit healthy SAN activation, and initiate SAN microreentrant waves [38]. Previous simulation studies have observed that microreentrant conduction was not sustained by AP changes, and required a large center of fibrotic tissue to produce microreentry [22]. We identified a specific range of low coupling that sustained emergent spiral wave behavior in a heterogeneous grid of SAN APs, indicating that increased coupling is a crucial component to synchronization of pacemaker cells and sufficient coupling may override differences in cell-to-cell transmembrane dynamics. Notably, our simulations did not require implementing a “track” of fibrotic tissue around which the AP wave could propagate, but still resulted in emergent behavior. The synchrony factor measurement, which best represents how closely together cells are firing within a beat, demonstrated a non-linear relationship with respect to coupling (Figure 8A). This further supports the idea that coupling non-linearly alters the ability of SAN cells to fire synchronously.

In previous work, Michaels et al. examined the effects of cell-to-cell coupling strength on entrainment [20, 33]. They tested both paired cells as well as small grids. In paired cells they observed that the cells tended to synchronize to a CL closer to the faster cell. In a grid they found that the apparent wave front slowed as coupling strength decreased, however they did not see spiral waves or other conduction issues. They also tested a grid with a partial wall of inexcitable tissue and found that the cells on the other side were still trained, although slightly delayed.

Shifts in the location and size of the SAN pacemaker may occur as a compensatory mechanism in response to sinus node dysfunction, vagal nerve stimulation, or pharmacological block of the Na⁺ current or L-type Ca²⁺ current [31]. Our simulations show that a heterogeneous SAN with low coupling will sustain pro-arrhythmic behavior, but increasing coupling may help synchronize the entire grid. A shift in size and location of the pacemaker may be beneficial due to coupling changes; this shift transforms the pacemaker into a larger group of highly coupled cells, which our simulations show can synchronize through a democratic entrainment process regardless of ion channel heterogeneity. These studies also suggest that altering gap junction coupling in the SAN may promote healthy pacemaking activity.

The studies presented here perform a variety of parameter sensitivity analyses in order to deconstruct the relationship between ion channel conductance, conductivity, SAN transmembrane properties (DDR, MDP, peak voltage, and CL), and the proposed measurements of synchrony (synchrony factor and peak voltage similarity). Both I_Ca,L and I_Kr demonstrated a high contribution to transmembrane dynamics and non-linear behavior with respect to transmembrane properties and synchrony metrics. This is further supported by experimental evidence of sinus node impairment or dysfunction related to modulation of L-type Ca²⁺ current [39, 40] or hERG channel function [41, 42]. The parameter sensitivity analysis also demonstrates that the relationship between specific ion channel conductances may vary depending on cell-to-cell coupling values (Figures 5, 8). This suggests that it is not necessarily sufficient to extrapolate effects of single cell perturbation to emergent behavior at the tissue level.

The emergent behavior of coupled oscillators has been widely explored in both computational and experimental studies of multiple areas of biology such as mitochondrial, circadian rhythms, synaptic firing, and broader ecological studies. Synchronization and its quantification has been widely discussed in networks of coupled oscillators [43–46]. Our hope is that these set of studies contributes to an already diverse set of work and adds to understanding of the impact of ion channel behaviors as well as coupling in the SAN pacemaker. We also believe that the synchrony metrics presented here would be useful for quantifying dynamics in larger tissue experiments such as optical mapping experiments. Future studies quantifying synchronization of coupled SAN oscillators in tissue or determining the impact of ion channel changes on generation of microreentrant arrhythmias may help support the findings in the simulations shown here.

Limitations

While these mathematical modeling studies are based on a well-validated single cell model of the rabbit SAN AP, the two-dimensional simulations have important limitations based on experimental data. For the sake of simplicity, SAN cells were coupled in a uniform rectangular grid with homogeneous coupling strengths, but this does not match the detailed physiology of the three-dimensional atrium. Similarly, heterogeneity of the SAN is modeled as either a gradient with AP differences between the central and periphery of the node, or a mosaic with a variable mix of SAN and atrial cells from periphery to the center. The gradient model is supported by a wide range of experimental data and simulations showing a change in the transmembrane properties of the SAN between the periphery and the center, a change in the density of ion channels responsible for I_Na and I_f, and a lack of atrial cells in the center of the SAN [47, 48]. The studies presented here more closely represent the mosaic model, but are distinct in that only SAN cells are implemented (no randomly placed atrial cells are simulated in the grid simulations). It is also important to note that the grid used in our studies contains a relatively small number of cells compared to the actual SAN. However, based on previous work [36, 37], it is possible to consider each cell in the grid as representative of a group of cells so that the behavior observed in our grid should scale to larger dimensions. Finally, these simulations did not implement parasympathetic stimulation of tissue, or patch of atrial tissue surrounding the SAN to further explore activation of atrial tissue by the SAN complex. This is especially important to note in two-dimensional simulations using a rabbit SAN model, since the architecture of the rabbit sinus node and its subsequent conduction pathways is distinct from human [31]. Conductivities between cells in these simulations were fixed, so the tubular shape of cells was ignored. In addition, we observed emergent behavior at the edges of large grid simulations, which may be an artifact of the simulation setup. While outside the scope of the current study, going forward it would be interesting to design experiments to test model predictions in in ex vivo preparations.

Author Contributions

DG, BO, AD, and TH participated in design of study; DG, AD, and BO performed simulations, analyzed data, and wrote the manuscript; DG, BO, and TH revised the manuscript.

Funding

This work was supported by NIH (grant numbers HL114893, HL134824, HL135096 to TH; HL129766 to BO); James S. McDonnell Foundation (to TH); Saving tiny Hearts Society (to TH); and a TriFit Challenge grant from Ross Heart Hospital and Davis Heart and Lung Research Institute.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2018.00063/full#supplementary-material

Supplementary Videos 1–6. Movies of cell voltages in 7 × 7 grids depicting asynchronous activity (Video 1), spiral waves (Video 3), and synchronous activity (Video 5) from Figure 7. As propagation across the grid happens very quickly, the movies are slowed by a factor of 2X. The end of the spiral wave from Figure 7B thus happens at t = 40 s in the video as opposed to t = 20 s in the simulation. Additional movies depicting clusters in the same 7 × 7 grids from Figure 7 are labeled Videos 2, 4, 6 for asynchronous, spiral waves, and synchronous activity, respectively. These movies show clusters of cell action potentials forming and then being removed. Colors for the clusters repeat regularly and are not for any purpose besides distinguishing the clusters. As the time for a cluster to propagate is very fast, these movies are slowed by a factor of 4X. All movies were created at a constant 40 fps using python plotting library (matplotlib).

Footnotes

1. ^https://hundlab.engineering.osu.edu

2. ^https://github.com/hundlab

References

1. Dobrzynski H, Boyett MR, Anderson RH. New insights into pacemaker activity: promoting understanding of sick sinus syndrome. Circulation (2007) 115:1921–32. doi: 10.1161/CIRCULATIONAHA.106.616011

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Mangoni ME, Nargeot J. Genesis and regulation of the heart automaticity. Physiol Rev. (2008) 88:919–82. doi: 10.1152/physrev.00018.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Quinn TA, Kohl P. Mechano-sensitivity of cardiac pacemaker function: pathophysiological relevance, experimental implications, and conceptual integration with other mechanisms of rhythmicity. Prog Biophys Mol Biol. (2012) 110:257–68. doi: 10.1016/j.pbiomolbio.2012.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Dobrzynski H, Anderson RH, Atkinson A, Borbas Z, D'Souza A, Fraser JF, et al. Structure, function and clinical relevance of the cardiac conduction system, including the atrioventricular ring and outflow tract tissues. Pharmacol Ther. (2013) 139:260–88. doi: 10.1016/j.pharmthera.2013.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Semelka M, Gera J, Usman S. Sick sinus syndrome: a review. Am Fam Physician (2013) 87:691–6.

PubMed Abstract | Google Scholar

6. Gosney M, Harper A, Conroy S. Geriatric Medicine. Oxford: Oxford University Press (2012).

PubMed Abstract

7. Vardas PE, Simantirakis EN, Kanoupakis EM. New developments in cardiac pacemakers. Circulation (2013) 127:2343–50. doi: 10.1161/CIRCULATIONAHA.112.000086

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Aschoff J. (ed.). (1981). “Annual rhythms in man,” in Biological Rhythms (Boston, MA: Springer).

9. Unudurthi SD, Wolf RM, Hund TJ. Role of sinoatrial node architecture in maintaining a balanced source-sink relationship and synchronous cardiac pacemaking. Front Physiol. (2014) 5:446. doi: 10.3389/fphys.2014.00446

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Wolf RM, Glynn P, Hashemi S, Zarei K, Mitchell CC, Anderson ME, et al. Atrial fibrillation and sinus node dysfunction in human ankyrin-B syndrome: a computational analysis. Am J Physiol Heart Circ Physiol. (2013) 304:H1253–66. doi: 10.1152/ajpheart.00734.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Sanders P, Kistler PM, Morton JB, Spence SJ, Kalman JM. Remodeling of sinus node function in patients with congestive heart failure: reduction in sinus node reserve. Circulation (2004) 110:897–903. doi: 10.1161/01.CIR.0000139336.69955.AB

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Lou Q, Hansen BJ, Fedorenko O, Csepe TA, Kalyanasundaram A, Li N, et al. Upregulation of adenosine A1 receptors facilitates sinoatrial node dysfunction in chronic canine heart failure by exacerbating nodal conduction abnormalities revealed by novel dual-sided intramural optical mapping. (2014) Circulation 130:315–24. doi: 10.1161/CIRCULATIONAHA.113.007086

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Fedorov VV, Glukhov AV, Chang R. Conduction barriers and pathways of the sinoatrial pacemaker complex: their role in normal rhythm and atrial arrhythmias. Am J Physiol Heart Circ Physiol. (2012) 302:H1773–83. doi: 10.1152/ajpheart.00892.2011

PubMed Abstract | CrossRef Full Text

14. Glynn P, Onal B, Hund TJ. Cycle length restitution in sinoatrial node cells: a theory for understanding spontaneous action potential dynamics. PloS ONE (2014) 9:e89049. doi: 10.1371/journal.pone.0089049

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Glynn P, Unudurthi SD, Hund TJ. Mathematical modeling of physiological systems: an essential tool for discovery. Life Sci. (2014) 111:1–5. doi: 10.1016/j.lfs.2014.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Wilders R. 25 Years of SA nodal cell modeling. In: Annual International Conference of the IEEE Engineering in Medicine and Biology – Proceedings. Lyon (2007), p. 152–5. doi: 10.1109/IEMBS.2007.4352245

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Maltsev VA, Yaniv Y, Maltsev AV, Stern MD, Lakatta EG. Modern perspectives on numerical modeling of cardiac pacemaker cell. J Pharmacol Sci. (2014) 125:6–38. doi: 10.1254/jphs.13R04CR

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Maltsev AV, Maltsev VA, Stern MD. Stabilization of diastolic calcium signal via calcium pump regulation of complex local calcium releases and transient decay in a computational model of cardiac pacemaker cell with individual release channels. PLoS Comput Biol. (2017) 13:e1005675. doi: 10.1371/journal.pcbi.1005675

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Fabbri A, Fantini M, Wilders R, Severi S. Computational analysis of the human sinus node action potential: model development and effects of mutations. J Physiol. (2017), 595:2365–96. doi: 10.1113/JP273259

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Michaels DC, Matyas EP, Jalife J. Mechanisms of sinoartrial pacemaker synchronization: a new hypothesis. Circ Res. (1987) 61:704–14. doi: 10.1161/01.RES.61.5.704

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Greisas A, Zlochiver S. Modulation of cardiac pacemaker inter beat intervals by sinoatrial fibroblasts - a numerical study. Conf. Proc. IEEE Eng. Med. Biol. Soc. (2016) 2016:165–8. doi: 10.1109/EMBC.2016.7590666

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Kharche SR, Vigmond E, Efimov IR, Dobrzynski H. Computational assessment of the functional role of sinoatrial node exit pathways in the human heart. PloS ONE (2017) 12:e0183727. doi: 10.1371/journal.pone.0183727

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Onal B, Gratz D, Hund T. LongQt: a cardiac electrophysiology simulation platform. MethodsX (2016) 3:589–99. doi: 10.1016/j.mex.2016.11.002

PubMed Abstract | CrossRef Full Text

24. Kurata Y, Hisatome I, Imanishi S, Shibamoto T. Dynamical description of sinoatrial node pacemaking: improved mathematical model for primary pacemaker cell. Am J Physiol Heart Circ Physiol. (2002) 283:H2074–101. doi: 10.1152/ajpheart.00900.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Ohio Supercomputer Center. Ohio Supercomputer Center. Columbus OH (1987).

PubMed Abstract

26. Sobie EA. parameter sensitivity analysis in electrophysiological models using multivariable regression. Biophys J. (2009) 96:1264–74. doi: 10.1016/j.bpj.2008.10.056

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Onal B, Gratz D, Hund TJ. Ca²⁺/calmodulin kinase II-dependent regulation of atrial myocyte late Na⁺ current, Ca²⁺ cycling and excitability: a mathematical modeling study. Am J Physiol Heart Circ Physiol. (2017) 313:H1227–39. doi: 10.1152/ajpheart.00185.2017

CrossRef Full Text | Google Scholar

28. Unudurthi SD, Wu X, Qian L, Amari F, Onal B, Li N, et al. Two-pore K+ channel TREK-1 regulates sinoatrial node membrane excitability. J Am Heart Assoc. (2016) 5:e002865. doi: 10.1161/JAHA.115.002865

PubMed Abstract | CrossRef Full Text

29. Sobie EA, Sarkar AX. Regression methods for parameter sensitivity analysis: applications to cardiac arrhythmia mechanisms. 2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, IEEE (2011), p. 4657–60. doi: 10.1109/IEMBS.2011.6091153

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Sarkar AX, Christini DJ, Sobie EA. Exploiting mathematical models to illuminate electrophysiological variability between individuals. J Physiol. (2012), 590:2555–67. doi: 10.1113/jphysiol.2011.223313

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Boyett MR, Honjo H, Kodama I. The sinoatrial node, a heterogeneous pacemaker structure. Cardiovasc Res. (2000) 47:658–87. doi: 10.1016/S0008-6363(00)00135-8

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Abramovich-Sivan S, Akselrod S. Phase response curve based model of the SA node: simulation by two-dimensional array of pacemaker cells with randomly distributed cycle lengths. Med Biol Eng Comput. (1999) 37:482–91. doi: 10.1007/BF02513334

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Michaels DC, Matyas EP, Jalife J. Dynamic interactions and mutual synchronization of sinoatrial node pacemaker cells: a mathematical model. Circ Res. (1986) 58:706–20. doi: 10.1161/01.RES.58.5.706

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Fedorov VV, Glukhov AV, Chang R, Kostecki G, Aferol H, Hucker WJ, et al. Optical mapping of the isolated coronary-perfused human sinus node. J Am Coll Cardiol. (2010) 56:1386–94. doi: 10.1016/j.jacc.2010.03.098

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Monfredi O, Dobrzynski H, Mondal T, Boyett MR, Morris GM. The anatomy and physiology of the sinoatrial node-a contemporary review. Pacing Clin Electrophysiol. (2010) 33:1392–406. doi: 10.1111/j.1540-8159.2010.02838.x

PubMed Abstract | CrossRef Full Text | Google Scholar

36. James TN. The sinus node. Am J Cardiol. (1977) 40:965–86. doi: 10.1016/0002-9149(77)90048-0

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Archiv P, Joumal E. Membrane currents in the rabbit sinoatrial node cell as studied by the double microelectrode method. Pflfigers Arch. (1976) 364:45–52. doi: 10.1007/BF01062910

CrossRef Full Text | Google Scholar

38. Glukhov AV, Fedorov VV, Anderson ME, Mohler PJ, Efimov IR. Functional anatomy of the murine sinus node: high-resolution optical mapping of ankyrin-B heterozygous mice. Am J Physiol Heart Circ Physiol. (2010) 299:H482-91. doi: 10.1152/ajpheart.00756.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Yada H, Murata M, Shimoda K, Yuasa S, Kawaguchi H, Ieda M, et al. Dominant negative suppression of rad leads to QT prolongation and causes ventricular arrhythmias via modulation of L-type Ca²⁺ channels in the heart. Circ Res. (2007) 101:69–77. doi: 10.1161/CIRCRESAHA.106.146399

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zhang H, Sun AY, Kim JJ, Graham V, Finch EA, Nepliouev I, et al. STIM1–Ca ²⁺ signaling modulates automaticity of the mouse sinoatrial node. Proc Natl Acad Sci USA. (2015) 112:E5618–27. doi: 10.1073/pnas.1503847112

CrossRef Full Text

41. Swan H, Viitasalo M, Piippo K, Laitinen P, Kontula K, Toivonen L. Sinus node function and ventricular repolarization during exercise stress test in long QT syndrome patients with KvLQT1 and HERG potassium channel defects. J Am Coll Cardiol. (1999) 34:823–9.

PubMed Abstract | Google Scholar

42. Nawathe PA, Kryukova Y, Oren RV, Milanesi R, Clancy CE, Lu JT, et al. An LQTS6 MiRP1 mutation suppresses pacemaker current and is associated with sinus bradycardia. J Cardiovasc Electrophysiol. (2013) 24:1021–7. doi: 10.1111/jce.12163

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Winfree AT. The Geometry of Biological Time. New York, NY: Springer New York (2001).

44. Osipov GV, Kurths J, Zhou C. Synchronization in Oscillatory Networks. Berlin, Heidelberg: Springer Berlin Heidelberg (2007).

45. Berridge MJ, Rapp PE. A comparative survey of the function, mechanism and control of cellular oscillators. J Exp Biol. (1979) 81:217–79.

PubMed Abstract | Google Scholar

46. Kawamura Y. Phase synchronization between collective rhythms of fully locked oscillator groups. Sci Rep. (2015) 4:4832. doi: 10.1038/srep04832

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Zhang H, Holden AV, Boyett MR. Gradient model versus mosaic model of the sinoatrial node. Circulation (2001) 103:584–8. doi: 10.1161/01.CIR.103.4.584

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Cloherty SL, Dokos S, Lovell NH. Qualitative support for the gradient model of cardiac pacemaker heterogeneity. 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference, Shanghai: IEEE (2005) p. 133–6.

PubMed Abstract | Google Scholar