Abstract
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 Ca2+ 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 (Vm,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 Vm, 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
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 (ICa,L, ICatt, Ih, IKr, IKs, INCX, INaK, Ito) 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, Vm,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 (Cell1,1, Cell1,2, Cell5,5, Cell1,1, Cell1,3), where Celli, j is located in the ith row and jth column of a two-dimensional grid, the clustering algorithm would select Cell1,1 as the beginning of a new cluster, C1, as there are no other existing clusters. The second cell in the sequence, Cell1,2, would then be added to C1 as it is within three cells of Cell1,1. Cell5,5, however, is too far away from C1 and so would be marked as the beginning of a new cluster, C2. Since Cell1,1 is already assigned to C1, its appearance a second time in the sequence causes the algorithm to initialize a new cluster C3. Finally, the last element in the sequence, Cell1,3, would be added to C3.
Vm,peak similarity is calculated as the inverse of the standard deviation of Vm,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 website1, and are accessible through Github2. 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: ICa,L = 3.19472, ICatt = 2.59237, Ih = 2.64054, IKr = 2.3018, IKs = 2.93799, INCX = 2.79239, INKA = 2.45686, Ito = 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 3
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 Ca2+ current (ICa,L) and the transport rate of the Na+/K+ ATPase (INaK) have a large positive influence on Vm,peak in both single and coupled simulations, while perturbations in conductances of the T-type Ca2+ current (ICatt), rapid delayed rectifier K+ current (IKr), and transient outward K+ current (Ito) are inversely related, i.e., an increase in ICatt, IKr, or Ito decrease peak Vm. Despite the overall agreement, there are small but interesting differences between sensitivity of the single vs. coupled cell. First, while ICa,L has a positive effect on CL and Vm,peak in both single and coupled cells, its influence is diminished in the coupled cell. Likewise, coupling reduces the influence of ICa,L on DDR and MDP. In contrast, our simulations predict that coupling increases the effect of IKr, at least with respect to DDR and CL.
Figure 4
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 Vm,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., ICa,L, ICa,tt, Ih for CL and Ih and IKr 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
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, Vm,peak, and CL for both the coupled WT and variant cells at each coupling value. As expected, values for MDP, Vm,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
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
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.
Vm,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) Vm,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 Vm,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. Vm,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 Ih (Figure 8A). However, ICa,L and INaK both had a large negative contribution and Ito a large positive contribution to Vm,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−4 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
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 ICa,L and IKr; (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 Ca2+ 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 ICa,L and IKr 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 Ca2+ 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 INa and If, 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.
Statements
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.
Conflict of interest
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–6Movies 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).
References
1.
DobrzynskiHBoyettMRAndersonRH. New insights into pacemaker activity: promoting understanding of sick sinus syndrome. Circulation (2007) 115:1921–32. 10.1161/CIRCULATIONAHA.106.616011
2.
MangoniMENargeotJ. Genesis and regulation of the heart automaticity. Physiol Rev. (2008) 88:919–82. 10.1152/physrev.00018.2007
3.
QuinnTAKohlP. 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. 10.1016/j.pbiomolbio.2012.08.008
4.
DobrzynskiHAndersonRHAtkinsonABorbasZD'SouzaAFraserJFet 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. 10.1016/j.pharmthera.2013.04.010
5.
SemelkaMGeraJUsmanS. Sick sinus syndrome: a review. Am Fam Physician (2013) 87:691–6.
6.
GosneyMHarperAConroyS. Geriatric Medicine. Oxford: Oxford University Press (2012).
7.
VardasPESimantirakisENKanoupakisEM. New developments in cardiac pacemakers. Circulation (2013) 127:2343–50. 10.1161/CIRCULATIONAHA.112.000086
8.
AschoffJ. (ed.). (1981). Annual rhythms in man, in Biological Rhythms (Boston, MA: Springer).
9.
UnudurthiSDWolfRMHundTJ. Role of sinoatrial node architecture in maintaining a balanced source-sink relationship and synchronous cardiac pacemaking. Front Physiol. (2014) 5:446. 10.3389/fphys.2014.00446
10.
WolfRMGlynnPHashemiSZareiKMitchellCCAndersonMEet 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. 10.1152/ajpheart.00734.2012
11.
SandersPKistlerPMMortonJBSpenceSJKalmanJM. Remodeling of sinus node function in patients with congestive heart failure: reduction in sinus node reserve. Circulation (2004) 110:897–903. 10.1161/01.CIR.0000139336.69955.AB
12.
LouQHansenBJFedorenkoOCsepeTAKalyanasundaramALiNet 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) Circulation130:315–24. 10.1161/CIRCULATIONAHA.113.007086
13.
FedorovVVGlukhovAVChangR. 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. 10.1152/ajpheart.00892.2011
14.
GlynnPOnalBHundTJ. Cycle length restitution in sinoatrial node cells: a theory for understanding spontaneous action potential dynamics. PloS ONE (2014) 9:e89049. 10.1371/journal.pone.0089049
15.
GlynnPUnudurthiSDHundTJ. Mathematical modeling of physiological systems: an essential tool for discovery. Life Sci. (2014) 111:1–5. 10.1016/j.lfs.2014.07.005
16.
WildersR. 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. 10.1109/IEMBS.2007.4352245
17.
MaltsevVAYanivYMaltsevAVSternMDLakattaEG. Modern perspectives on numerical modeling of cardiac pacemaker cell. J Pharmacol Sci. (2014) 125:6–38. 10.1254/jphs.13R04CR
18.
MaltsevAVMaltsevVASternMD. 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. 10.1371/journal.pcbi.1005675
19.
FabbriAFantiniMWildersRSeveriS. Computational analysis of the human sinus node action potential: model development and effects of mutations. J Physiol. (2017), 595:2365–96. 10.1113/JP273259
20.
MichaelsDCMatyasEPJalifeJ. Mechanisms of sinoartrial pacemaker synchronization: a new hypothesis. Circ Res. (1987) 61:704–14. 10.1161/01.RES.61.5.704
21.
GreisasAZlochiverS. Modulation of cardiac pacemaker inter beat intervals by sinoatrial fibroblasts - a numerical study. Conf. Proc. IEEE Eng. Med. Biol. Soc. (2016) 2016:165–8. 10.1109/EMBC.2016.7590666
22.
KharcheSRVigmondEEfimovIRDobrzynskiH. Computational assessment of the functional role of sinoatrial node exit pathways in the human heart. PloS ONE (2017) 12:e0183727. 10.1371/journal.pone.0183727
23.
OnalBGratzDHundT. LongQt: a cardiac electrophysiology simulation platform. MethodsX (2016) 3:589–99. 10.1016/j.mex.2016.11.002
24.
KurataYHisatomeIImanishiSShibamotoT. Dynamical description of sinoatrial node pacemaking: improved mathematical model for primary pacemaker cell. Am J Physiol Heart Circ Physiol. (2002) 283:H2074–101. 10.1152/ajpheart.00900.2001
25.
Ohio Supercomputer Center. Ohio Supercomputer Center. Columbus OH (1987).
26.
SobieEA. parameter sensitivity analysis in electrophysiological models using multivariable regression. Biophys J. (2009) 96:1264–74. 10.1016/j.bpj.2008.10.056
27.
OnalBGratzDHundTJ. Ca2+/calmodulin kinase II-dependent regulation of atrial myocyte late Na+ current, Ca2+ cycling and excitability: a mathematical modeling study. Am J Physiol Heart Circ Physiol. (2017) 313:H1227–39. 10.1152/ajpheart.00185.2017
28.
UnudurthiSDWuXQianLAmariFOnalBLiNet al. Two-pore K+ channel TREK-1 regulates sinoatrial node membrane excitability. J Am Heart Assoc. (2016) 5:e002865. 10.1161/JAHA.115.002865
29.
SobieEASarkarAX. 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. 10.1109/IEMBS.2011.6091153
30.
SarkarAXChristiniDJSobieEA. Exploiting mathematical models to illuminate electrophysiological variability between individuals. J Physiol. (2012), 590:2555–67. 10.1113/jphysiol.2011.223313
31.
BoyettMRHonjoHKodamaI. The sinoatrial node, a heterogeneous pacemaker structure. Cardiovasc Res. (2000) 47:658–87. 10.1016/S0008-6363(00)00135-8
32.
Abramovich-SivanSAkselrodS. 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. 10.1007/BF02513334
33.
MichaelsDCMatyasEPJalifeJ. Dynamic interactions and mutual synchronization of sinoatrial node pacemaker cells: a mathematical model. Circ Res. (1986) 58:706–20. 10.1161/01.RES.58.5.706
34.
FedorovVVGlukhovAVChangRKosteckiGAferolHHuckerWJet al. Optical mapping of the isolated coronary-perfused human sinus node. J Am Coll Cardiol. (2010) 56:1386–94. 10.1016/j.jacc.2010.03.098
35.
MonfrediODobrzynskiHMondalTBoyettMRMorrisGM. The anatomy and physiology of the sinoatrial node-a contemporary review. Pacing Clin Electrophysiol. (2010) 33:1392–406. 10.1111/j.1540-8159.2010.02838.x
36.
JamesTN. The sinus node. Am J Cardiol. (1977) 40:965–86. 10.1016/0002-9149(77)90048-0
37.
ArchivPJoumalE. Membrane currents in the rabbit sinoatrial node cell as studied by the double microelectrode method. Pflfigers Arch. (1976) 364:45–52. 10.1007/BF01062910
38.
GlukhovAVFedorovVVAndersonMEMohlerPJEfimovIR. 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. 10.1152/ajpheart.00756.2009
39.
YadaHMurataMShimodaKYuasaSKawaguchiHIedaMet al. Dominant negative suppression of rad leads to QT prolongation and causes ventricular arrhythmias via modulation of L-type Ca2+ channels in the heart. Circ Res. (2007) 101:69–77. 10.1161/CIRCRESAHA.106.146399
40.
ZhangHSunAYKimJJGrahamVFinchEANepliouevIet al. STIM1–Ca 2+ signaling modulates automaticity of the mouse sinoatrial node. Proc Natl Acad Sci USA. (2015) 112:E5618–27. 10.1073/pnas.1503847112
41.
SwanHViitasaloMPiippoKLaitinenPKontulaKToivonenL. 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.
42.
NawathePAKryukovaYOrenRVMilanesiRClancyCELuJTet al. An LQTS6 MiRP1 mutation suppresses pacemaker current and is associated with sinus bradycardia. J Cardiovasc Electrophysiol. (2013) 24:1021–7. 10.1111/jce.12163
43.
WinfreeAT. The Geometry of Biological Time. New York, NY: Springer New York (2001).
44.
OsipovGVKurthsJZhouC. Synchronization in Oscillatory Networks. Berlin, Heidelberg: Springer Berlin Heidelberg (2007).
45.
BerridgeMJRappPE. A comparative survey of the function, mechanism and control of cellular oscillators. J Exp Biol. (1979) 81:217–79.
46.
KawamuraY. Phase synchronization between collective rhythms of fully locked oscillator groups. Sci Rep. (2015) 4:4832. 10.1038/srep04832
47.
ZhangHHoldenAVBoyettMR. Gradient model versus mosaic model of the sinoatrial node. Circulation (2001) 103:584–8. 10.1161/01.CIR.103.4.584
48.
ClohertySLDokosSLovellNH. 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.
Summary
Keywords
sinoatrial node, synchrony, computational modeling, coupling, ion channel
Citation
Gratz D, Onal B, Dalic A and Hund TJ (2018) Synchronization of Pacemaking in the Sinoatrial Node: A Mathematical Modeling Study. Front. Phys. 6:63. doi: 10.3389/fphy.2018.00063
Received
15 December 2017
Accepted
04 June 2018
Published
22 June 2018
Volume
6 - 2018
Edited by
Flavio H. Fenton, Cornell University, United States
Reviewed by
Fei Geng, McMaster University, Canada; Yanyan Ji, Georgia Institute of Technology, United States
Updates
Copyright
© 2018 Gratz, Onal, Dalic and Hund.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Thomas J. Hund thomas.hund@osumc.edu
This article was submitted to Biomedical Physics, a section of the journal Frontiers in Physics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.