A mathematical model of action potentials of mouse sinoatrial node cells with molecular bases.

Genetically modified mice are popular experimental models for studying the molecular bases and mechanisms of cardiac arrhythmia. A postgenome challenge is to classify the functional roles of genes in cardiac function. To unveil the functional role of various genetic isoforms of ion channels in generating cardiac pacemaking action potentials (APs), a mathematical model for spontaneous APs of mouse sinoatrial node (SAN) cells was developed. The model takes into account the biophysical properties of membrane ionic currents and intracellular mechanisms contributing to spontaneous mouse SAN APs. The model was validated by its ability to reproduce the physiological exceptionally short APs and high pacing rates of mouse SAN cells. The functional roles of individual membrane currents were evaluated by blocking their coding channels. The roles of intracellular Ca(2+)-handling mechanisms on cardiac pacemaking were also investigated in the model. The robustness of model pacemaking behavior was evaluated by means of one- and two-parameter analyses in wide parameter value ranges. This model provides a predictive tool for cellular level outcomes of electrophysiological experiments. It forms the basis for future model development and further studies into complex pacemaking mechanisms as more quantitative experimental data become available.

the cardiac excitation sequence is initiated by the cardiac pacemaker, the sinoatrial node (SAN), and spreads throughout the heart. The function of the SAN is therefore understandably vital in cardiac physiology. The SAN produces spontaneous rhythmic action potentials (APs) via an emergent action of many interacting membrane ionic currents and intracellular processes (8). Experimental cell membrane ion channel patch-clamp data were first incorporated into a predictive mathematical model by Noble in 1962 (63), leading to the seminal “theory leading experiment.” It resulted in decades of extensive experimental discovery of a large number of cell membrane ion channels (2, 6, 15, 17, 18) and the study of pacemaking mechanisms in several mammalian species (9, 10, 20, 44, 81). These experimental data have been incorporated into increasingly complex mathematical models of SAN automaticity to simulate species-specific SAN AP waveforms (83) and have also been used to investigate the roles of ionic currents and Ca2+ handling in cardiac pacemaking mechanisms (50, 78, 88, 90).

Despite these extensive electrophysiological studies, understanding of the underlying molecular basis of cardiac arrythmogenesis remains limited. Pore-forming transmembrane proteins facilitating ionic currents are encoded by specific ion channel genes (54, 68). Tissue-specific expression of these ion channel genes and their functional impact on the tissue's physiological (e.g., SAN pacemaking) and pathological (e.g., bradycardia) states are important issues to be addressed. In recent times, genetically modified mouse models have gained popularity in the study of clinical SAN dysfunction (39, 42). The mouse heart rate [∼294 beats/min (10)], however, is much higher than that in larger mammals [rabbit: ∼189 beats/min (14) to ∼214 beats/min (74); human: ∼73 beats/min (81)]. To fully appreciate such observations, the biophysical properties of ion channels and their genetic bases that regulate mouse SAN APs need further investigation.

A biophysically detailed mathematical model of mouse SAN APs will complement experimentation in underpinning the correlation between ion channel function and their genetic bases. However, most contemporary SAN models (13, 16, 36, 67, 89) do not consider isoform-specific channel gating and kinetics. A recent membrane delimited mouse SAN cell model (56) partially reproduced some experimental data of the mouse SAN but has fundamental limitations in simulating the functional roles of ion currents on mouse cardiac pacemaking (31). The present study aimed to further develop a biophysically detailed model of mouse SAN isolated cell APs with biophysically distinct ion channel isoforms and a detailed dynamic intracellular ionic mechanism.

MODEL DEVELOPMENT

The present model was primarily based on experimental data from isolated mouse SAN cells, from which biophysical parameters of Hodgkin-Huxley formulations of ionic currents were estimated. The model consists of 38 coupled ODEs describing voltage-gated ion channel currents, pump and exchanger currents, dynamic variations of intracellular ionic concentrations of Ca2+ ([Ca2+]i and [Ca2+]sub), K+ ([K+]i), and Na+ ([Na+]i), and Ca2+ dynamics. A schematic diagram of the model is shown in Fig. 1. Model parameters for the basal cell model are shown in Table 1. The Glossary shows all abbreviations used in this report.

Fig. 1.

Schematic diagram of the mouse SAN cell models. The depolarizing (inward arrows) and repolarizing (outward arrows) cell membrane currents with isoformal detail along with the Ca2+ handling mechanism are shown. For abbreviations in this and all other figures, see the Glossary.

Table 1.

Model parameter values in the basal cell model

Model Parameters	Values
Cm, pF	25
[Ca2+]o, mM	1.8
[CM]tot, mM	0.045
[CQ]tot, mM	10
ECaL, mV	47
ECaT, mV	45
Est, mV	17
ENa,1.5, mV	41
F, C/M	96.485
gb,Ca, nS/pF	0.0006
gb,Na, nS/pF	0.00486
gb,K, nS/pF	0.001
gCaL,1.2, nS/pF	0.24
gCaL,1.3, nS/pF	0.72
gCaT, nS/pF	0.5586
gh, nS/pF	0.228
gK1, nS/pF	0.0324
gKr, nS/pF	0.09456
gKs, nS/pF	0.01196
gNa,1.5, nS/pF	2.37 × 10−4
gNa,1.1, nS/pF	2.37 × 10−4
gst, nS/pF	0.0024
gsus, nS/pF	0.0156
gto, nS/pF	0.1968
I_NaK¯, pA/pF	5.698
kbCM, ms−1	0.542
kbCQ, ms−1	0.445
kbTC, ms−1	0.446
kbTCM, ms−1	0.00751
kbTMM, ms−1	0.751
kfCM, mM−1/ms	237.7
kfCQ, mM−1/ms	0.534
kfTC, mM/ms	88.8
kfTCM, mM/ms	237.7
kfTMM, mM/ms	2.277
KmfCa, mM	0.00035
Km,K, mM	1.4
Km,Na, mM	14
KNaCa, pA/pF	220
K1ni, mM	395.3
K1no, mM	1628
K2ni, mM	2.289
K2no, mM	561.4
K3ni, mM	26.44
K3no, mM	4.663
Kci, mM	0.0207
Kco, mM	3.663
Kcni, mM	26.44
Qci	0.1369
Qco	0
Qn	0.4315
[K+]o, mM	5.4
Kup, mM	0.0006
[Mg2+]i, mM	2.5
[Na+]o, mM	140
Pup, mM/ms	0.04
nup	2
Kmf, mM	0.000246
Kmr, mM	3.29
R, J · mol−1 · K−1	8.314
T, K	310.5
[TC]tot, mM	0.031
[TMC]tot, mM	0.062
Vcell, pL	3
Vi, pL	1.346
Vrel, pL	0.0036
Vsub, pL	0.03328
Vup, pL	0.0348
αfCa	0.021
τCa,dif, ms	0.04
τtr, ms	40
EC50_SR	0.45
KoCa, Mm−2	10
Kom, ms−1	0.06
KiCa, mM−1 ms−1	0.5
Kim, ms−1	0.005
Ks, ms−1	130 × 104
HSR	2.5
MaxSR	15
MinSR	1

For model parameter definitions, see the Glossary.

Cell Morphology

The cell volume of mouse pacemaking SAN myocytes ranges from 0.9 pl in spider-shaped cells to 5 pl in spindle-shaped cells (10, 39, 55). Our previous experimental recordings showed that cell Cm ranges from 20 to 60 pF (39). Therefore, a cell volume of 3 pl and a Cm of 25 pF were assumed in the present model for a primary SAN cell. The intracellular volumes relate to the cell volume in proportion to those in the Kurata et al. (36) model.

Descriptions of Membrane Currents and Ionic Homeostasis

The membrane currents considered in the mouse SAN cell model consist of INa (INa,1.1 and INa,1.5), ICaL (ICaL,1.2 and ICaL,1.3), ICaT, Ist, If (If,HCN4), IKr, IK1, INaCa, INaK, and other miscellaneous currents (IKs, I, Ito, and background currents). Wherever possible, biophysical parameters for these ionic currents were derived from experimental data of mouse SAN cells. Due to the limited availability of experimental data from the mouse SAN for some ionic channel currents, some parameters in the model were either based on parent models (36, 37, 89) or estimated to reproduce ion channel I-V relationships, experimentally observed mouse SAN AP features, and the functional role of individual modeling components. The formulations for cyclic oscillations of intracellular ionic concentrations and Ca2+ buffering in the cell compartments were adapted from Kurata et al. (36). To improve the effectiveness of SR function on modulating pacemaking APs, SERCA Ca2+ uptake and RyR Ca2+ release function mechanisms (i.e., Ca2+ pump-release mechanisms) were based on the formulation by Shannon et al. (73). A summary of the sources of experimental and modeling data is shown in Supplemental Material in Table S1.

INa (INa,1.1 and INa,1.5).

The Nav1.1 and Nav1.5 isoforms predominantly compose mouse SAN INa (41). Nav1.1 is uniform across the SAN. Nav1.5 is spatially heterogeneous, with its channel current density increasing from the center toward the periphery (39). Both activate during the pacemaker potentials; Nav1.1 activates at more positive potentials than Nav1.5.

The models for the steady-state gating variables of INa were based on our experimental studies on isolated mouse SAN cells (Fig. S1, A and B) (39, 41). INa time kinetics in the present model were based on parent SAN models (36, 89). Maximum ionic current conductances for Nav1.1 and Nav1.5 isoforms were computed by comparing simulated INa I-V relationships to experimental data (Fig. S1C) (39, 41). Current traces of both isoforms during the voltage-clamp simulation are shown in Fig. S1D.

ICaL (ICaL,1.2 and ICaL,1.3).

The two isoforms contributing to ICaL in the mouse SAN are the weakly expressed Cav1.2 isoform and the strongly expressed Cav1.3 isoform (52, 56).

In the model, the steady-state inactivation of Cav1.2 was based on experimental data (52), and the steady-state activation was adopted from a study by Mangoni et al. (56) (Fig. S2A). The steady states of activation and inactivation of Cav1.3 (Fig. S2B) were based on rabbit SAN experimental data as determined by Zhang et al. (89). Experimental data for mouse ICaL time kinetics are as yet unavailable. It was found that SAN AP profiles depend strongly on ICaL inactivation time kinetics at higher voltages (32). In the parent models, the inactivation time constant had a very large peak value (∼310 ms at −40 mV) in the Kurata et al. (36) model and was kept almost constant at ∼45 ms (over −40 mV) in the Zhang et al. (89) model. Since such a slow voltage-dependent inactivation of ICaL produces a prolonged APD, the inactivation time kinetics of ICaL of the Zhang et al. (89) model were adapted in this study. The reversal potentials for Cav1.2 and Cav1.3 were set to +47 mV based on experimental I-V data (10, 52).

The ionic current models were validated by simulating experimental ICaL I-V relationships (52) for Cav1.2 and Cav1.3 isoforms (Fig. S2C). Current traces are shown in Fig. S2D. The maximum ionic current conductances so determined for both ICaL isoforms (52) were found to give current densities significantly lower than in other reports [ICaL,1.2: 3.1 pA/pF at +10 mV (56) and ICaL,1.3: 10.4 pA/pF at −20 mV (56)]. The conductances were therefore augmented to give appropriate current densities during model pacemaking APs.

ICaT.

ICaT is a voltage-dependent gated Ca2+ current with fast inactivation kinetics. ICaT is expressed in the SAN of most species (22, 56) with a predominant isoform of Cav3.1 (7). In the mouse SAN, ICaT has been shown to activate at approximately −65 mV, peak at approximately −30 mV, and have a reversal potential of +45 mV (10). The steady-state activation of ICaT was adopted from previous modeling studies (36, 89). The steady-state gating properties are shown in Fig. S3A, which shows a small window current at −45 mV, close to the late phase of mouse SAN AP diastolic depolarization. The time kinetics of ICaT were derived from experimental data (22) in previous models of ICaT in the rabbit SAN. The ionic current model was validated by quantitatively comparing the simulated ICaT I-V data with experimental data (10, 52, 56), as shown in Fig. S3B. The fast kinetics of ICaT are evident from the current traces shown in Fig. S3C.

If.

In the mouse SAN, the isoforms contributing to the native If are a heteromeric combination of HCN1, HCN2, and HCN4, with all being permeable to Na+ and K+ (3). It has been shown that the biophysical properties of individual isoforms coding If channels are different from each other (2), and their tandems (e.g., HCN1–HCN4 tandem) have distinct kinetics to native If. In the mouse SAN, it has been shown that HCN1 and HCN2 activate 25 times faster than HCN4 (61), and HCN4 is the predominant isoform compared with HCN1 and HCN2 (45). Therefore, a single native If was modeled in the present study. With experimental data showing a large range in the measured V1/2 of If [−107 mV (10), −87.4 mV (55), −105 mV (26), −101 mV (1), and −106.8 mV (1)], a value of V1/2 = −106.8 mV was adopted in the present study (Fig. S4A). The slope of If activation was also based on experimental data (10, 66) and was taken to be −16.1 mV in the basal model. If time kinetics were based on experimental data of native If from mouse SAN cells (10), as shown in Fig. S4B. The model for If was validated by quantitatively comparing the simulated If I-V data with experimental data (10, 55, 66) (Fig. S4C). Simulated current traces for If are shown in Fig. S4D.

Ist.

A novel SAN Na+ current with ICaL-like kinetics was first reported by Guo et al. (21) and has been observed in the SAN of several species (10, 21, 74). The voltage dependence of its gating kinetics are not available except for in the rat SAN. The gating kinetics of Ist were therefore based on experimental data obtained from the rat SAN (74). Ist in the mouse SAN activates around −80 mV and peaks at −60 mV (10). Both values are ∼10 mV negative to those of the rat SAN. Therefore, a −10-mV shift was implemented in the steady-state gating kinetics to fit experimental mouse SAN I-V data. Mouse SAN Ist has a reversal potential of around +17 mV. The Ist model reproduced experimental I-V data (10) (Fig. S4E) and current traces (Fig. S4F).

IKr.

The primary protein regulating mouse SAN IKr is mERG1 (11). The model for IKr steady-state activation was based on experimental data (10, 58), whereas the steady state of inactivation was adapted from a previous study (36) (Fig. S5A). The activation time kinetics were based on experimental data (11). Since the experimental data were acquired at 28°C, a Q10 temperature correction of 1.4 (11, 60) was applied (Fig. S5B). IKr inactivation kinetics were adapted from a previous rabbit model (36). The IKr model was validated by simulating experimental I-V data (Fig. S5C) and current traces (Fig. S5D) (10, 55).

IK1.

Maximum ionic current conductance of the mouse SAN IK1 model was based on experimental data of Cho et al. (10) obtained from mouse SAN cells and included a dependence on [K+]o.

Miscellaneous membrane currents.

The other currents that the model consists of are IKs, 4-AP-sensitive Isus and Ito, INaCa, INaK, and Ib, each of which were either based on experimental data or inherited from parent models. For example, the steady-state activation and kinetics of IKs were based on experimental data obtained from guinea pig SAN cells (24, 58). Equations for Isus and Ito were based on parent rabbit SAN cell models (36, 89). The model for INaCa was adopted from Kurata et al. (36) and simulated the saturation characteristics of INaCa at large values of [Ca2+]sub. The scaling parameter kNaCa was set to 220 pA/pF in the basal model to simulate physiological oscillations of [Ca2+]i and [Ca2+]sub (Fig. S6) and to facilitate a strong coupling between the intracellular and membrane pacemaking mechanisms.

Ionic Homeostasis in Mouse SAN Cells

Na+ and K+ concentrations.

The simultaneous direct measurement of APs and intracellular ion concentrations is a formidable task experimentally (59), providing limited data to facilitate a quantitative description of the dynamic intracellular variations of [Na+]i and [K+]i (see for, e.g., Refs. 16 and 28) in a mathematical model of cardiac cells. The present basal model incorporated the dynamics of [Na+]i and [K+]i based on an ionic material balance. INaK was optimized to simulate stable oscillations of the intracellular ion concentrations.

Ca2+ homeostasis, buffers, and SR function.

There is growing experimental evidence of the role of intracellular Ca2+ in modulating pacemaking APs of the mouse SAN (57, 62) via INaCa, which contributes to the diastolic depolarization phase of the AP (9, 84). Intracellular Ca2+ has also been shown to play an important role in the late diastolic depolarization phase of cat pacemaking APs (29, 44) via ICaT. Depletion of Ca2+ by challenging RyR function (i.e., reduced Ca2+ release from the SR) has been shown to reduce pacemaking in isolated mouse SAN cells (62). Mathematical modeling has also shown that a strong interaction between intracellular Ca2+ and membrane ionic processes regulates pacemaking APs (6, 50). In the present model, we started by incorporating the Ca2+ buffering formulations of Kurata et al. (36). Second, the SERCA Ca2+ uptake mechanism was based on the formulation by Shannon et al. (73), which models Ca2+ uptake from the cytosol to the NSR. The values of parameters for the SERCA pump vary widely among various modeling studies (69–71, 73). In this study, the Hill coefficient was set to 2 and the affinity of the forward mode of the pump to 0.246 μM to reflect measurements in rodents (69). The other parameters were optimized to maintain an ∼0.1 μM peak free [Ca2+]i in the cytosol (36). Finally, the Ca2+ release mechanism to simulate RyR function was also based on the formulation of the Shannon et al. (73) multistate RyR “Ca2+ clock” model. The RyR Ca2+ release model was optimized by Maltsev and Lakatta (50) for rabbit SAN cell pacemaking APs with longer a CL of ∼ 315 ms. Since the CL of mouse SAN cell APs is much shorter (∼200 ms) and is regulated by several ionic channel currents with much larger densities, we therefore further modified the parameters of the Shannon et al. (73) model for modeling the Ca2+ clock in our mouse SAN cell model.

Numerical Integration and Stability of Solutions

Integration method.

The presented mouse SAN cell model consists of a 38 variable coupled stiff nonlinear ODE system. A standard explicit Euler integration method with a constant time step (dt = 10−3 ms) was found to give stable solutions. Further reduction of the time step or implementation of a fourth-order Runge-Kutta integration method gave the same simulation results within 0.1% tolerance levels while increasing the computational time and effort. Simulations were carried out on a Sun Sparc with Solaris 9 using C programming.

Free parameters values, stability of solution, and determination of initial conditions.

Biophysically detailed cardiac models consist of a large numbers of variables (∼10–100) and parameters (in excess of 100) (32). The basal model was primarily based on experimental data associating the majority of the modeling parameters to experimentally observed electrophysiological data. Biophysical parameters not available in experimental data were taken from the parent models (36, 37, 89) and adapted to give typical AP characteristics in mouse SAN cells. For details of all inherited features, see Table S1.

For the determination of stable intracellular ionic oscillations, mass balance based on the ionic conservation principle was applied (33, 77). During the time integration, the accumulation or depletion of all intracellular ions was monitored. Oscillations of all intracellular ionic concentrations were stabilized within a tolerance of 1% variation. For example, the concentrations of [Na+]i and [K+]i were stabilized around average physiological values of 8 and 140 mM, respectively, by minimally adjusting the maximum INaK. It should be emphasized that the free parameters were estimated by minimal modifications to basal values. Consideration of dynamic variations in [Na+]i and [K+]i did not give rise to increased degeneracy and drift in the presented models (35), as seen predominantly in atrial and ventricular models (28, 46). This is due to the mass balance-based parametric fine tuning and the spontaneously beating SAN APs not requiring an extracellular stimulus. Such parametric fine tuning will nonetheless allow some drift and has been minimized in this study, as described above. The degeneracy and drift could be limited by maintaining constant [Na+]i and [K+]i. This would, however, have rendered the function of INaK to be basically a background current. It should be noted that model development is, essentially, an iterative process, with novel experimental data being incorporated into the models continually.

Once the free parameters were determined, the initial conditions and stability of the solutions were estimated. To determine the initial conditions, all 38 variables were set to values expected at MDP. Thus, the value of Vm was set to −65 mV. Initial values of 0 and 1 were assigned to the activation and inactivation gates, respectively. All ionic concentrations were set at anticipated physiological diastolic levels (e.g., [Ca2+]i = 0.1 μM, [Na+]i ∼ 8 mM, and [K+]i ∼ 140 mM). The initial conditions for the RyR model four variables (O, I, R, and RI) were set to suitable nonzero values. The dynamic spontaneously oscillating model was then integrated for 1,000 s. In the final 100 s of simulated activity, the beat-to-beat variation of all the 38 variables was monitored. When the variation was below tolerance values, the solution was deemed to be stable. Values of all variables at MDP were then noted as initial conditions, as shown in Table 2.

Table 2.

Initial conditions for the basal cell model

Model Variables	Initial Values
V, mV	−64.6508005232
qa	0.6185734202
qi	0.4571815638
dT	0.0015910627
fT	0.4307958088
pa	0.4033900505
pi	0.9254742549
xs	0.0127072059
fL,1.2	0.9969251824
dL,1.2	0.0000044417
fL,1.3	0.9812943664
dL,1.3	0.0001993197
fCa	0.7679518810
r	0.0045587632
q	0.6140737444
m1.5	0.3991763968
h1.5	0.2765325493
j1.5	0.0252682289
m1.1	0.1069853058
h1.1	0.4548050890
j1.1	0.0272302605
y	0.0281231662
[Ca2+]i, mM	0.0000326892
[Ca2+]rel, mM	0.1038669181
[Ca2+]up, mM	0.1038669181
[Ca2+]sub, mM	0.0000534741
fTC	0.0000534741
fTMC	0.1771999779
fTMM	0.7268211576
fCMs	0.0231294040
fCMi	0.0142601584
fCQ	0.1067137942
[Na+]i, mM	8.10311
[K+]i, mM	139.9221122
R	0.7717873638
O	0.0000000737
I	0.0000000207
RI	0.2173099253

For model variable definitions, see the Glossary.

Typically found in gating kinetics and in ion current formulations, which are all functions of Vm, division by small numbers close to machine precision makes the models unsuitable for numerical analysis. All above integrations were carried out by eliminating singularities in the modeling logical expressions (46) by performing a first-order Taylor expansion around all the relevant singularity values, and the approximations held valid with a range of 0.05 mV around the singularity.

RESULTS

Figure 2 shows the simulated time course of mouse SAN APs together with dynamic variation of the major ion channel currents and intracellular ion concentrations. The simulated APs have a rapid depolarizing upstroke, short APD, and fast pacemaking rate, which agree with the unique features of isolated mouse SAN cell APs. Figure S6 shows details of the dynamic variations of intracellular Ca2+ variables. It can be seen that the model produces [Ca2+]i oscillations with an amplitude of 0.1 μM, which is close to experimental data (36). [Ca2+]sub oscillations have systolic and diastolic values of 1,744 and 69.6 nM, respectively, which are comparable with experimental observations in rabbit SAN cells (51). To validate the model, simulated APs and ion channel blocking simulations were compared with experimental data as detailed below.

Fig. 2.

Model-generated APs, ionic currents, and [Ca2+]i transients in the basal model. A: model-generated APs. B–M: ionic currents during APs, normalized to Cm (in units of pA/pF). N: simulated [Ca2+]i transients. O: dynamic [K+]i variation during the AP. P: dynamic [Na+]i variation during the AP.

Simulated APs and Experimental Recordings

Simulated APs were validated against experimental recordings, as shown in Fig. 3. Simulated AP characteristics, i.e., CL, APD50, APD90, OS, MDP, dV/dtmax, TOP, and DDR, agreed well with the experimental data, which further validated the use of parameters in the equations of individual ionic currents. Quantitative results are shown in Table 3 along with quantitative experimental AP characteristics and Mangoni et al. (56) model simulation.

Fig. 3.

Mouse SAN APs and characteristics. A: experimental AP recording from mouse SAN cells (solid line) (Cm ∼ 26 pF) from our laboratory. B: simulated SAN AP in the basal model. C: simulated AP from Mangoni et al. (56) mouse SAN model. D: simulated MDP, OS, and TOP (open bars), experimental data (●), and Mangoni et al. (56) model data (×). The same nomenclature for modeling and experimental data apply to E and F. E: simulated APD50, APD90, and CL, experimental data, and Mangoni et al. (56) model data. F: simulated dV/dtmax and DDR along with experimental data and Mangoni et al. (56) model data. Data in D–F were obtained from quantitative results or from analysis of single cell AP profiles by Cho et al. (10), Lei et al. (39, 41), Mangoni et al. (52, 55), Clark et al. (11), Zhang et al. (91), Verheijck et al. (80), Alig et al. (1), Wu et al. (84), and Chen et al. (9).

Table 3.

Experimental data and simulated mouse SAN AP characteristics from Mangoni et al. (56) and our basal models

	CL, ms	APD50, ms	APD90, ms	OS, mV	MDP, mV	dV/dtmax, V/s	TOP, mV	DDR, V/s
Experimental data range	106.1 to 255	26.3 to 56.88	45.8 to 107.14	22.7 to 42.0	−70 to −52	8.2 to 32	−54.6 to −40	0.172 to 0.34
Mangoni et al. model	212.31	104.00	121.1	24.6	−66.87	5.31	−47.3	0.3701
Present cell model	212.21	35.04	62.00	22.60	−64.52	9.20	−48.45	0.1786

Experimental data were obtained from quantitative results or from analysis of single cell AP profiles by Cho et al. (10), Lei et al. (39, 41), Mangoni et al. (52, 55), Clark et al. (11), Zhang et al. (91), Verheijck et al. (80), Alig et al. (1), Wu et al. (84), and Chen et al. (9). For definitions of the AP characteristics, see the Glossary.

Roles of Individual Ionic Channel Currents in the Genesis of Mouse SAN APs

The functional roles of ion channels were investigated by isoform-specific ion channel blocking simulations and comparison with experimental data.

Role of INa.

Simulation of block of INa,1.1 and INa,1.5 in the model is shown in Fig. 4. As shown in Fig. 4, block of INa,1.1 did not affect the model pacemaking (CL increased from 212.12 to 217.01 ms), but block of INa,1.5 augmented CL by 32.5% (CL = 281.14 ms). Block of both isoforms simultaneously augmented CL by 38.2% (CL = 293.27 ms). The simulations are in agreement with our experimental results (41), where a 27% increase of CL was observed due to inhibition of INa using 30 μM TTX. A similar functional role of INa in the rabbit SAN has been observed by Kurata et al. (37). Although INa is not the main upstroke current in SAN cells and does not affect dV/dtmax (Fig. 4, bottom), it contributes to pacemaking since both isoforms are active in the late diastolic depolarization pacemaker potential range (Fig. S1, A and B). This was reflected in the reduced DDR when INa was blocked in our basal model (Fig. 4). Since INa acts predominantly in the late diastolic depolarization phase, block of INa did not affect APD and marginally reduced dV/dtmax in the model. Experimentally, however, we have observed that low (100 nM) or high (30 μM) concentrations of TTX used to block INa,1.1 or INa,1.5, respectively, caused a reduction of dV/dtmax. This may be due to TTX affecting other depolarizing Ca2+ currents in SAN cells, as has been seen in experimental studies (75, 87).

Fig. 4.

Functional role of INa in the mouse SAN cell model. A, top: AP profiles under control (solid line), Nav1.1 blocked (dashed line, almost coincidental with the control), and Nav1.5 blocked (dashed-dotted line) conditions. Bottom, dV/dt during APs. B: CL (top) and DDR (bottom) under control and Nav1.1, Nav1.5, and INa blocked conditions. Shown are the model simulation (open bars), experimental data (●) (41), and Mangoni et al. (56) model simulation data (×).

To further investigate the functional role of INa,1.5 in mouse SAN pacemaking, the effects of alterations of gNa,1.1 and gNa,1.5 were simulated in the model (Fig. S7). Reduction of either isoform-specific conductance slowed down pacemaking, whereas increased conductance caused the expected progressive acceleration of pacemaking.

Role of ICaL.

Block of ICaL,1.2 in the model (Fig. 5) marginally increased CL by 0.2%, which was accompanied by reduced dV/dtmax, APA, and APD. Our simulations were similar to experimental data from rabbit SAN cells (34), where block of ICaL by nifedipine in small balls from central SAN tissue accelerated pacemaking, reduced dV/dtmax, and reduced APA. Since the activation potential of ICaL,1.2 is higher than TOP, its functional role in pacemaking activity is primarily in the modulation of OS, rather than the pacemaking rates (Fig. S8). The limited role of ICaL,1.2 on APs is due to its relatively small window current and small current density compared with ICaL,1.3 (Fig. S2). The Mangoni et al. (56) model showed a small reduction of pacemaking rates upon inhibition of ICaL,1.2.

Fig. 5.

Functional effects of block of ICaL isoforms. A, top: AP profiles under control (solid line) and Cav1.2 blocked (dashed line) conditions. Bottom, dV/dt. B, top: AP profiles under control (solid line), partial block of Cav1.3 by 45% (dotted line), and total Cav1.3 blocked (dashed line) conditions. Bottom, dV/dt during APs. C: CL under control, Cav1.2 blocked, and CaV1.3 blocked conditions. Shown are the model simulation (open bars), experimental data (●), and Mangoni et al. (56) model simulation data (×).

The functional role of Cav1.3 on mouse SAN pacemaker activity was studied by blocking ICaL,1.3 (Fig. 5). Experimentally, Cav1.3−/− mice showed significantly slower pacemaking with prolonged CL or even an arrest in pacemaking activity (52). Such a negative chronotropic effect due to reduced ICaL,1.3 on heart rates has also been observed in vivo, where Cav1.3 knockout mice showed prolonged PR intervals with a high incidence of sinus bradycardia (91). In the present model, block of ICaL,1.3 by 45% produced a 17.33% increase in CL. Complete block abolished autorhythmic APs, leading to a stable resting potential of −40.47 mV. The Mangoni et al. (56) model also showed a similar arrest of pacemaking upon ICaL,1.3 block (stable resting potential at −29.3 mV). When both ICaL,1.2 and ICaL,1.3 were blocked in the model, pacemaking in the model was arrested with a stable resting potential of −40.8 mV. Changes in CL are shown in Fig. 5.

The effects of altered gCaL,1.3 on pacemaking were simulated in the model. ICaL,1.3 is a vital regulator of OS and APD, as shown in Fig. S8. A decrease in ICaL,1.3 slowed down pacemaking APs, resulting in increased CL along with dramatically reduced APAs. Augmenting ICaL,1.3 accelerated pacemaking APs, resulting in reduced CL along with increased APAs and dV/dtmax.

Role of ICaT.

Experiments using Cav3.1 knockout mice have shown a pacemaking rate reduction by 34% (56). Studies in cat atria by Huser et al. (29) and Zhou and Lipsius (92) have also shown that a Ni2+-induced reduction of ICaT causes a 35–230% prolongation of CL, indicating the critical importance of ICaT in SAN pacemaking. In the model, complete ICaT block dramatically slowed down the pacemaking, with CL increasing by 49.1% (Fig. 6). As observed experimentally (29), ICaT block revealed that the reduction of CL was directly related to the reduced level of diastolic Ca2+ (Fig. 6). In model simulations, ICaT block produced a 38% reduction of peak [Ca2+]i and a 42.79% reduction of [Ca2+]sub as well as reduced Ca2+ release during DDR. The simulated changes of [Ca]i and [Ca]sub in response to ICaT block were greater than the experimental data of Huser et al. (29), who observed an initial reduction of the [Ca2+]sub pedestal during the diastolic depolarization phase but no remarkable change in peak [Ca2+]sub or [Ca2+]i. This is possibly due to the fact that in the model, ICaT is a major contributor to DDR and it plays an important role (secondary to ICaL) in the early upstroke. Such a discrepancy between the model simulation and experimental data suggest a possible limitation of the model, which requires further improvements when more experimental data on the biophysical properties of mouse ICaT become available. As ICaT activates during the DDR phase, the increase in CL arising from ICaT block is primarily due to a decrease in DDR (see Fig. S9). Our simulation results are quantitatively similar to the experimentally observed role of ICaT. ICaT block produced a 6.2% increase of CL in the Mangoni et al. (56) model. Previous models of rabbit SAN cells have assumed a small maximal conductance of ICaT, limiting its role in pacemaking dynamics. However, ICaT has been experimentally observed to play a significant role in mouse SAN pacemaking and an even greater role in cat SAN cell pacemaking, as observed by Huser et al. (29). In the present model, the important role of ICaT is reflected in its effects on CL and DDR (Fig. S9, A and B). Furthermore, the role of ICaT was investigated in a two-parameter analysis, where gh and Pup (parameters regulating the diastolic depolarization, along with gCaT) were simultaneously varied in the presence and absence of ICaT. As shown in Fig. S9C, the absence of ICaT augmented the CL of stable pacemaking APs dramatically and reduced model robustness. Furthermore, the absence of ICaT caused pacemaking arrest or aperiodic firing in a larger parametric range. The parameter analysis revealed the critical role of ICaT in assuring the robustness of mouse SAN pacemaking activities.

Fig. 6.

Functional effects of ICaT block. A: AP profiles under control (solid lines) and ICaT blocked (dashed lines) conditions. B: [Ca2+]i (top) and [Ca2+]sub (bottom) profiles during APs under control (solid lines) and ICaT blocked (dashed lines) conditions. C: CL (top) and DDR (bottom) under control and ICaT blocked conditions. Shown are the model simulation (open bars), experimental data (●) (56), and Mangoni et al. (56) model simulation data (×).

Role of Ist.

Ist has been observed in pacemaking SAN cells of several species (10, 20, 21). Ist has been observed in the mouse (10) and was incorporated as a small maximum conductance current in the model. A larger conductance of Ist caused the model to become unstable. Block of Ist in the model (Fig. 7) marginally slowed the pacemaking activity and increased CL by 2.1%. The alterations in CL due to Ist block are shown in Fig. 7.

Fig. 7.

Functional effects of Ist block. A: AP profiles under control (solid lines) and Ist blocked (dashed lines) conditions in the model. B: CL under control and Ist blocked conditions. Shown are the basal cell model simulation (open bars) and Mangoni et al. (56) model data (×).

Role of If.

The role of If in cardiac pacemaking remains controversial. A recent study (50) has argued the important role of If in SAN pacemaking. However, HCN4 gene knockout experimental studies (4, 26) in the mouse SAN have shown that the role of If is significant, albeit not isolated from other mechanisms governing pacemaking. The functional role of If on mouse SAN pacemaking activity was determined by block of If (Fig. 8). Complete block of If slowed down pacemaking, with CL being increased by 13.72% in the model. It produced a more hyperpolarized MDP to −69.4 mV, a hyperpolarization of 5 mV. Since If modulates pacemaking by primarily acting during the DDR phase of the AP, it did not affect OS and APD. The measured CL and DDR under complete block conditions are shown in Fig. 8, B and C, respectively. Effects of alterations of gh on model CL and DDR are shown in Fig. S10.

Fig. 8.

Functional effects of block of If. A: AP profiles under control (solid line) and total If blocked (dashed line) conditions in the model. B: CL under control and If blocked conditions. C: DDR under control and If blocked conditions. Shown are the cell model simulation (open bars) and Mangoni et al. (56) model data (×).

Role of IKr.

The role of IKr in the present model is to modulate APD and APA. Partial block of IKr by 48% prolonged APD90 by 27.96% in the model. Accompanied with APD prolongation, block of IKr also elevated MDP by +4.2 mV in the model. Note that partial block of IKr only marginally reduced pacemaking (CL increased marginally by ∼0.7%; Fig. 9,i). Such a response to IKr block has been observed experimentally in the rabbit SAN (8) and has also been reproduced in parent rabbit SAN models (37, 89). This is due to the primary function of IKr being modulation of APD and APA. As IKr was progressively blocked, APA reduced progressively (Fig. 9, A,ii, and Fig. S11). The reduction in APA is attributable to the elevated MDP, which results in incomplete activation of INa and ICaL in the upstroke phase, giving reduced APA. Complete block of IKr abolished the pacemaking activity with Vm resting at −0.35 mV. Upon block of IKr by 48%, pacemaking in the Mangoni et al. (56) model was abolished.

Fig. 9.

Role of IKr, simulated effects of E-4031, and contribution of Ito in mouse SAN pacemaking. Solid lines denote control AP profiles, and open bars show present model simulation data. ×, Mangoni et al. (56) model data; ●, experimental data. A,i: AP profiles under control, partial 48% IKr blocked (dotted line), and total IKr blocked (dashed line) conditions in the cell model. A,ii: CL (top) and APA (bottom) under control, 48% IKr blocked, and total IKr blocked conditions. B,i: AP profiles under control and E-4031 (dotted line) conditions in the cell model. B,ii: CL (top) and APA (bottom) under control and E-4031 conditions. C,i: AP profiles under control and Ito blocked (dotted line) conditions in the cell model. C,ii: CL (top) and APA (bottom) under control and Ito blocked conditions.

The present model showed a marginal reduction of CL due to IKr block, whereas our experimental data (62) showed that application of 1 μM E-4031 reduced CL by 68%. In isolated mouse hearts, Cho et al. (10) found that 1 μM E-4031 reduced CL by 36.5%. A low (0.2 μM) concentration of E-4031 was found to inhibit cellular IKr in the mouse SAN by 60% and increase CL by 51.5% in the experimental study by Clark et al. (11). In intact mouse hearts, the same study found an application of 1 μM E-4031 to increase CL in the range of 31–212%. E-4031 is presumed to be an IKr-selective drug. However, Verheijck et al. (79) showed that E-4031 not only inhibits IKr strongly but also reduces ICaL by 30% in the rabbit SAN. Admittedly, the 30% block of ICaL was observed at a high (10 μM) concentration of E-4031. In the present model, a substantial reduction of ICaL (45%) was found to be necessary in the model to give a prolonged CL. Other major K+ currents (i.e., IK1, Ito, IKs, and Isus) may also be sensitive to E-4031 but were found to play insignificant roles in the model's pacemaking behavior (Fig. S11). Therefore, the effects of E-4031 were simulated by reducing IKr density by 60% and reducing ICaL density by 50%. As shown in Fig. 9, B,i and ii, such alterations reduced pacemaking by 30.5%. APD90 was prolonged by 87.1%, and dV/dtmax was reduced by 79%. The present model as well as previous models (36, 56, 89) cannot reproduce CL prolongation by means of IKr block alone, which is a possible model limitation. To reproduce E-4031-induced prolongation of CL, reductions of IKr as well as ICaL are required in the present model. The full effect of E-4031 on mouse SAN ionic currents requires future experimental quantification.

Role of Ito.

It has been shown experimentally that Ito block causes an increase of APA and prolongs APD (40). Selective block of Ito in the model increased CL by 3.6% and prolonged APD50 by 10.2% and APD90 by 5.5% (Fig. 9, C,i and ii). Importantly, block of Ito caused OS to increase by 7 mV. The effect of Ito on OS reflects its role in regulating the early phase of repolarization. The small role of Ito in the present SAN cell model could be a model limitation (Fig. S11C). However, it should be noted that the role of Ito is greater in modulating repolarization in the ventricles compared with the atria (85, 86) and SAN (64, 65). To date, the biophysical properties of Ito in the mouse SAN have not been experimentally quantified. Furthermore, previous SAN modeling studies (36, 50, 89) have indicated that the role of Ito, compared with IKr, in regulating SAN AP characteristics is that of regulating OS and early repolarization. Further experimental quantification of Ito properties in mouse SAN cells will improve our understanding of Kv4.2/Kv4.3 Ito channels.

Role of INaCa.

INaCa strongly couples the AP to intracellular Ca2+ and contributes to the DDR phase of the SAN AP. Progressive reduction of INaCa reduced APA and DDR, eventually giving rise to aperiodic or chaotic oscillations. When INaCa was blocked by >70%, regular pacemaking activity was abolished. Complete block of INaCa set Vm to a resting value at −39.5 mV (Fig. 10). The effects of various levels of INaCa inhibition or augmentation on APA and DDR are shown in Fig. S12. With a decrease of INaCa, APA monotonically decreased. With severe block (>70%), APs became aperiodic and eventually pacemaking was arrested. The INaCa formulation in the Mangoni et al. (56) model was adopted from the Zhang et al. (89) model. Their formulation of INaCa does not provide the strong coupling between membrane and intracellular processes. Therefore, although block of INaCa in the Mangoni et al. (56) model produced a prolonged CL by 18.9%, it is not an essential component of the pacemaking mechanism in their model.

Fig. 10.

Effects of INaCa inhibition on pacemaking activity. A: AP profiles under control (solid line), partial 70% INaCa blocked (dotted line), and total INaCa blocked (dashed line) conditions in the model. APA (B) and DDR (C) under control and partial 70% INaCa blocked conditions. Shown are the model simulation (open bars) and Mangoni et al. (56) model data (×).

Role of intracellular Ca2+.

The role of intracellular Ca2+ release and uptake mechanisms was investigated in the model. In the basal model, SR Ca2+ SERCA uptake and RyR release mechanisms were based on the formulations by Shannon et al. (73). The isolated Ca2+ clock model has been previously investigated in detail (50) and produces a wide range of pacing rates suitable for the rabbit SAN cell model. We further explored the two-parameter space of the uptake rate parameter Pup and the release rate parameter ks to adapt the Ca2+ clock to our membrane clock of the mouse SAN cell model. The integrated model was then able to reproduce a wide range of stable CL APs regulated by the intracellular Ca2+-handling mechanism (Fig. 11). By identifying appropriate values of Pup and ks (shown by “+” in Fig. 11), the model was able to generate stable pacemaking APs and intracellular Ca2+ transients (Fig. 11). Upon block of SR uptake (Pup = 0), pacemaking was dramatically slowed down (Fig. 11). Ryanodine is known to slow down pacemaking in the mouse SAN (62) by means of challenging SR Ca2+ release (9). In the model, block of SR Ca2+ release (ks = 0) also had a severe effect on the AP, making it aperiodic (Fig. 11). Furthermore, when intracellular Ca2+ ([Ca2+]sub and [Ca2+]i) was buffered to diastolic concentrations ([Ca2+]i and [Ca2+]sub set to 50 nM, simulating the effects of BAPTA), pacemaking was arrested (Fig. 11).

Fig. 11.

Regulation of pacemaking by intracellular Ca2+ mechanisms. A: parametric analysis of basal cell firing as a function of major intracellular Ca2+ parameters for Ca2+ uptake (Pup) and Ca2+ release (ks). The “+” shows the parameter values used in the basal model. At low Pup or low ks, firing was found to be aperiodic or arrested. In B–D, the dashed gray line shows the 0-mV reference voltage. B: control APs. C: APs under Pup = 0 conditions. D: APs under ks = 0 conditions. E: APs under Ca2+ buffered at diastolic concentrations of 50 nM.

Effects of Iso on mouse SAN pacemaking activity.

The adrenergic response of the pacemaking mechanism in the mouse SAN has been studied using Iso experimentally. It has been shown that Iso causes a positive shift of the V1/2 of If. Alig et al. (1) observed a shift from −106 to −92 mV (a positive shift of +14 mV). Baruscotti et al. (4) saw a similar +7-mV shift, whereas Liao et al. (43) saw a shift from −128 to −110 mV (shift of +18 mV). The slope factors in the steady-state of activation were unaffected. Upon shift of the V1/2 from −106 to −92 mV in our model, CL was marginally reduced (CL = 208.1 ms). The effects of varied V1/2 and Pup on model pacemaking are shown in Fig. S13. As the simulated effect of shift of V1/2 was substantially smaller than the experimentally observed effects of the 20–40% increase in pacemaking rates (1, 4, 84), it might implicate a possible role of changes in the intracellular Ca2+ mechanism and other ionic currents affected by Iso. Experimentally, it has been found that ICaL is dramatically affected by Iso, resulting in augmented ICaL density (1) as well as a shift of the I-V curve to more negative potentials (72). Similarly, Iso also augments IKr density (by 12%) and shifts its steady state of activation to more negative potentials (82). Experimentally, it has also been observed that Iso increased INaK activity by increasing its sensitivity to [Na+]i (5, 15, 23), i.e., a reduction of the half-maximal [Na+]i in INaK (Km,Na). Iso also augmented SR uptake and release by stimulation of Ca2+/calmodulin-dependent protein kinase II (50, 73). In simulations of Iso, we adopt the approach of Shannon et al. (72) and reduced the [Ca2+]i affinity of the forward model of the SERCA pump (Kmf) and increased the baseline non-SR-dependent transition rate constant of RyR (KoCa) in the Ca2+ release mechanism. Therefore, in accordance with experimental data, the effects of Iso were simulated by the following: 1) shift of If V1/2 to −92 mV (1); 2) increasing current densities of both ICaL isoforms by 45% (1); 3) shift of the ICaL I-V curves toward more negative potentials by 5 mV (72); 4) reducing Km,Na from 14 to 11 mM (5, 15); 5) augmenting gKr by 12% (82); 6) shift of the IKr steady state of activation by 5 mV to more negative potentials (82); 7) reducing Kmf from 0.246 to 0.123 μM (72); and 8) increasing KoCa from 10 to 20 mM−2/ms (72). With these changes, simulated APs and Ca2+ transients under control and Iso conditions are shown in Fig. 12, A and B. Simulated actions of Iso produced a 23.21% reduction of CL (CL = 163.1 ms; Fig. 12). Along with the reduced CL, there was an increase in the [Ca2+]sub amplitude but a decrease in the amplitude of the cytosolic [Ca2+]i oscillations. In previous ventricular cell modeling studies (17, 19, 72), [Ca2+]i has been shown to be augmented due to Iso. [Ca2+]i is important in the ventricles for tension-force generation (17, 73). In the SAN, however, [Ca2+]sub, rather than [Ca2+]i, has been identified to play a more critical role in regulating pacemaking activity (9, 50, 51). With the above-detailed Iso-induced electrophysiological alterations, the Maltsev and Lakatta (50) model also showed a reduction of cystolic free [Ca2+]i and an increase of [Ca2+]sub. The differences in intracellular compartments and major currents between ventricular and pacemaking cells may explain such a difference between the two cells' responses of [Ca]i to Iso.

Fig. 12.

Effects of Iso. A: solid lines denote control, whereas dashed lines denote Iso simulations. Shown are simulated control AP (top), [Ca2+]sub (middle), and [Ca2+]i (bottom) profiles. B: CL under control (CL = 212.21 ms) and Iso (CL = 163.59 ms) simulations. C: DDR under control (DDR = 0.178 mV/ms) and Iso (DDR = 0.23 mV/ms) conditions.

DISCUSSION

The presented model successfully reproduced physiological mouse SAN pacemaking APs with the atypically short APD and high pacing rate, similar to experimental recordings (39, 41, 55, 56). The model was validated by its ability to reproduce the effects of several ion channel blocks or reduced current density due to gene knockouts as well as intracellular Ca2+ alterations on pacemaking APs (10, 25, 52, 55, 56). Its robustness was verified by means of two-parameter analyses. The model thus provides a tool for quantitatively evaluating the functional role of isoform-specific ion channel current on the generation and control of cardiac rhythm. Almost all contemporary SAN models simulate rabbit or guinea pig cell electrical activities, which have substantially different AP morphologies. In a previous study, Mangoni et al. (56) developed a model for the mouse SAN cell AP. However, the simulated APs have a long plateau and APD, a low OS, and a much larger DDR compared with experimental data (Fig. 3). The ionic currents in the Mangoni et al. model (56) were not validated against experimental data either. For example, their INa,1.1 and INa,1.5 models do not reproduce the experimental I-V data of the mouse SAN (41); the V1/2 of the steady-state If activation curve is positive compared with experimental data (10, 55, 66). In their model, ICaL does not have Ca2+-dependent inactivation, thus limiting the coupling between intracellular ionic concentrations and the pacemaking APs. Assuming constant values for intracellular ion concentrations, INaCa and INaK in the Mangoni et al. (56) model are analogs to background currents without functional relevance. Indeed, INaCa in the model has a much smaller effect on the pacemaking (increase of CL by 18.9% upon block of INaCa). Since intracellular Ca2+ dynamics were not incorporated, the model offers limited possibilities of studying the regulatory effects advanced mechanisms, such as that of the Ca2+ clock (50), on pacemaking. We therefore believe the present biophysically detailed model to be one of the first to reproduce the mouse SAN AP along with the functional roles of ionic currents and intracellular Ca2+ handling.

Comparison With Previous Models of Other Species

In recent decades, several mathematical models of the SAN AP have been developed for the rabbit (13, 16, 36, 89) and guinea pig (67). These models are based on a single channel voltage-clamp experimental data representing generic-type ion channels, without considering the species-specific and isoform-specific properties of the ion channels. Although these models provide useful tools to study the functional roles of a certain generic-type ion channel current, they are not suitable to evaluate the functional roles of isoform-specific channels on cardiac pacemaking. The present models have a similar structure to previous SAN models (13, 16, 36, 89) but incorporate the following major difference and important advances.

A major difference between the previous rabbit SAN models and the present mouse model is that the mouse AP has much faster pacemaking rates, a higher OS, and a larger dV/dtmax, which are attributable to the emergent interactions of ion channel properties specific to the mouse SAN, e.g., a substantially larger ICaL density in the mouse model and a greater ICaT density than that of ICaL, as found by Mangoni et al. (53).

A major advance in the present model is the introduction of isoform-specific models for the major ionic channel currents present in the SAN, which include INa (consisting of INa,1,1 and INa,1.5), ICaL (consisting of ICaL,1.2 and ICaL,1.3), ICaT (consisting of ICaT,3.1). Model equations and related parameters for these isoform-specified ion channel currents were based on experimental data measured from mouse SAN cells (2, 10, 11, 26, 41, 52, 55, 56, 66). Another advance is its ability to reproduce the experimentally observed effects of ion channel blockers or gene knockouts (INa, ICaL, ICaT, IKr, and If) on pacemaking APs of the mouse SAN (11, 39, 41, 52, 56).

Key Achievements of the New Model

A key achievement of the present model is that the simulated APs (Fig. 3) were comparable with those recorded from mouse SAN cells in terms of AP characteristics, including CL, APD50, APD90, OS, MDP, dV/dtmax, TOP, and DDR (see Table 3). These models could reproduce experimental data of voltage clamp for some major currents (Figs. S1–S5) as well as the functional roles of isoform-specific channels of INa and ICaL, which were more accurate than those by the previous Mangoni et al. (56) model. The functional roles of IKr and Ito have been explored in the model and indicate future directions of experimental quantification of mouse SAN electrophysiology. The simulated functional roles of all other channels also matched experimental data (see results and Figs. 3–11). The model also incorporated an intracellular Ca2+ handling mechanism including the Ca2+ clock, which, together with the membrane clock, allows the simulation of various experimentally observed findings regarding intracellular Ca2+. Therefore, the present model provides a basis for further simulating the initiation and conduction of cardiac pacemaking activity in the intact tissue of the mouse SAN and atrium in the future.

Limitations of the Model

Limitations of model components.

In the model, the functional effects of INa,1.1 (Fig. 4) and ICaL,1.2 (Fig. 5) were limited. This is partly due to the window currents of both these isoforms being small compared with other dominant isoforms (Figs. S1 and S2). Furthermore, these window currents line in the diastolic depolarization phase of the AP, where INaCa, If, and intracellular Ca2+ mechanisms play a predominant role. The biophysical properties of these isoforms will be improved as experimental data become available. A second limitation of the model is that IKr block did not affect pacemaking rates, in contrast to our experimental observations. In most contemporary models of the SAN, IKr block causes increased pacemaking along with a reduction of APA. This limitation may be due to the inadequate modeling of IKr using a Hodgkin-Huxley formulation (36), as stated below. The modeling of IKr and ICaT can also be further improved when more experimental data on their biophysical properties become available. Finally, the intracellular Ca2+ mechanism did not modulate pacemaking as strongly in the present model compared with the study of Maltsev and Lakatta (50). Although block of SR uptake or release showed the experimentally observed alterations in pacemaking, the augmentation of Pup in the present model had a smaller effect compared with the rabbit model (50). This is a consequence of the large current densities of ICaL, ICaT, and If. Furthermore, the present model consisted of ionic currents INa and IK1, which are absent in the model of Maltsev and Lakatta (50), which also contributes the differences in the functional roles of Ca2+ clocks between the two models. The role of the Ca2+ clock in SAN pacemaking is still under debate (27), and further experimental quantification of SAN cellular intracellular mechanisms is definitely required.

Lack of complete experimental data.

There is a lack of a complete data set on ion channel gating and kinetics, especially the intracellular Ca2+ handling mechanisms from the mouse SAN. In model development, we considered mouse SAN-specific biophysical data for some major currents but implemented generic current models for other channel currents (Table S1), which either have a small effect (e.g., IKs) or are primarily fine tuning tools for the model development (e.g., INaK and background currents). Although the mouse is a widely used animal model in the study of SAN and atrial arrythmias, the quantification of the biophysical properties of mouse membrane electrophysiology is still limited. In some instances, we need to use data obtained either at low temperature or from different species. Uncertainty about the accuracy of the experimental data further limited model functionality. The models also do not include some currents known to play important roles in the mouse SAN, for example, ACh-activated K+ current (18, 47), ATP-sensitive K+ current (18), a Ca2+-activated nonselective cationic current (12), and store-activated Ca2+ channel currents (30).

Model limitations due to Hodgkin-Huxley ion currents and lack of the Ca2+ sparks.

Hodgkin-Huxley formulation as basis of ion channel currents and passive diffusion of intracellular Ca2+ have already been shown to be limiting factors for mathematical models of cardiac APs (50). For the ion channel currents, two-state Hodgkin-Huxley formulation cannot faithfully simulate the activation, inactivation, and deactivation process as a multistate Markovian chain formulation does. However, the Markovian chain formulation requires extensive experimental data obtained at various protocols for validating model parameters. Otherwise, the parameter set may not be unique. To our best knowledge, such experimental data necessary for developing a well-validated Markovian chain model for mouse SAN ion channels are not yet unavailable. Therefore, two-state Hodgkin-Huxley formulation was implemented in the present models.

Ca2+ diffusion from the subspace to the myoplasm, Ca2+ release and uptake by the SR, and Ca2+ buffering mechanisms were adopted from the parent Kurata et al. (37) model and from Shannon et al. (73). It allowed simulation of the experimentally observed effects of Ca2+/calmodulin-dependent protein kinase II, cyclopiazonic acid, and RyRs (38, 48–50). Investigating the role of Ca2+ sparks in driving SAN cell pacemaking due to the lack of spatially extended considerations of Ca2+ dynamics, as done in our recent study (76), would be the next modeling challenge.

Despite these limitations, the models presented here represent a significant progress in developing biophysically detailed mathematical models for mouse SAN cells and form a basis for further development with the advent of new experimental data.

Conclusions

In this study, a mathematical model for pacemaking APs of the mouse SAN cell has been developed. This model forms an important step leading toward a “virtual whole heart” of the mouse. Using the model, we analyzed the functional role of individual ion channel currents and intracellular Ca2+ handling on generating mouse pacemaking APs. Our simulation data have shown that the genesis of pacemaking APs is a coordinated action of all ion channel currents (involving both depolarizing and repolarizing currents) and the intracellular Ca2+ clock, rather than being regulated by a single factor. For example, block of inward depolarizing currents, such as INa, If, ICaT, Ist, or INaCa, slowed down pacemaking APs, highlighting their functional contributions to mouse SAN pacemaking APs. Block of outward repolarizing current, such as IKr, abolished the pacemaking AP, indicating its important role as well. In the present model, intracellular Ca2+ handling has been shown to contribute to pacemaking APs (50) but does not underestimate the importance of HCN-related ionic pacemaking mechanisms (4).

GRANTS

This work was supported by Welcome Trust (UK) Grant WT/081809/Z/06/Z.

DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).