Effects of the heterozygous E161K mutation in SCN5A were implemented in the CellML code (Cuellar et al., 2003) of the Fabbri-Severi human SA nodal cell model (Fabbri et al., 2017) by scaling the fully-activated conductance of I_Na (g_Na) and shifting the steady-state I_Na activation and/or inactivation curves, as detailed below, based on the experimental data from literature described in the Introduction and summarized in Table 1. The slight change in steepness of the steady-state activation curve was also included. These modifications were applied to half of the intrinsic I_Na channels, thus representing the heterozygous nature of the mutation, taking into account that a functional I_Na channel is built with a single Na_V1.5 protein.

Identical changes were applied to the Maleckar et al. (2009) human atrial cell model. The latter model, which is also known as the “human atrial myocyte with new repolarization” (hAMr) model, was selected because it includes well-validated equations for the acetylcholine-activated potassium current I_K,ACh (Maleckar et al., 2008), thus allowing simulations over a wide range of acetylcholine levels. Action potentials were elicited with a 1 ms, ≈50% suprathreshold stimulus current at a frequency of 1 Hz. The threshold stimulus current amplitude was determined by increasing the stimulus current amplitude in 0.1 pA/pF steps until a train of 200 action potentials could successfully be elicited. To prevent slow drifts in ion concentrations, the intracellular Na⁺ and K⁺ concentrations were fixed, as were the cleft ion concentrations.

In H558 background, the E161K mutation was implemented by shifting the steady-state activation curve of the model I_Na by +19 mV. In R558 background, this curve was shifted by +13 mV, together with a -10 mV shift in the steady-state inactivation curve. In either background, a 50% decrease in g_Na was applied to arrive at an almost 4-fold decrease in peak current density of mutant I_Na during voltage clamp simulations. Furthermore, the slope factor of the steady-state activation curve of the model I_Na was increased by 20% to account for mutation-induced decrease in steepness of this curve.

The CellML code of both models, as available from the CellML Model Repository (Lloyd et al., 2008), was edited and run in version 0.9.31.1409 of the Windows based Cellular Open Resource (COR) environment (Garny et al., 2003b). All simulations were run for a sufficiently long time, i.e., for the duration of a train of 200 action potentials, to reach steady-state behavior. Data from the final ten action potentials were used for analysis.

First, we characterized I_Na of the Fabbri-Severi model of a single human SA nodal pacemaker cell (Fabbri et al., 2017) and its effect on the spontaneous activity of this cell. Figure 1A shows the steady-state activation and inactivation curves of I_Na (left) and the associated “window” (area of overlap; right). The window is in the voltage range of diastolic depolarization, between -60 and -40 mV. Because of the relatively slow changes in membrane potential of the SA nodal cell and the relatively fast kinetics of I_Na, this window is the major determinant of the course of I_Na during the action potential, of course in combination with the fully-activated conductance of I_Na (g_Na), which amounts to 22.3 nS in the Fabbri-Severi model. The associated steady-state current-voltage relationship of I_Na is shown in Figure 1B.

With a value near 0.6 pA, the amplitude of I_Na during diastolic depolarization is small. However, this amplitude is relatively large in comparison to that of the net membrane current (I_net) (Figure 1C, solid lines). It is therefore not surprising that full block of I_Na results in a 117-ms increase in cycle length from 813 to 930 ms (Figure 1C, dotted lines), corresponding with a 13% decrease in beating rate from 74 to 64 beats/min. In the simulated presence of 10 nM acetylcholine (ACh), the amplitude of I_net during diastolic depolarization becomes considerably smaller, whereas that of I_Na does barely change. Accordingly, full block of I_Na now results in a more prominent increase in cycle length, by 245 ms from 1027 to 1272 ms (Figure 1D), corresponding with a 19% decrease in beating rate from 58 to 47 beats/min. Of note, I_Na hardly affects the cycle length through a change in action potential duration.

Now that the simulations of Figure 1 had shown that changes in I_Na may modify the cycle length of the Fabbri-Severi model cell to a considerable extent, we first assessed the changes in peak and window I_Na caused by the E161K mutation in SCN5A, either in H558 or R558 background. Next, we tested the effects of these changes on the spontaneous pacemaker activity of the SA nodal model cell.

The left panel of Figure 2A illustrates the changes in the steady-state I_Na activation and inactivation curves due to the E161K mutation in H558 background. These changes are limited to a +19 mV shift in the steady-state I_Na activation curve and a slight decrease in steepness of this curve, which is poorly discernible in Figure 2A. As a result, the window of overlap is significantly reduced and shifted to less negative values of membrane potential, now roughly ranging from -50 to -30 mV (Figure 2A, right) and thus limiting the contribution of I_Na to diastolic depolarization. This is reflected in the current–voltage relationship of steady-state I_Na (Figure 2B, solid line), which is strongly reduced in comparison with wild-type I_Na.

In R558 background, the shift in the steady-state I_Na activation curve is less pronounced (+13 vs. +19 mV), but accompanied by a -10 mV shift in the steady-state inactivation curve (Figure 2C, left). As a result, the window of overlap is smaller than in H558 background (Figure 2C, right). However, it better fits with the voltage range of diastolic depolarization. Yet, the steady-state I_Na in this voltage range is even more reduced in comparison with wild-type I_Na than in H558 background (Figure 2D, solid line).

As set out in the Introduction, voltage clamp experiments on SCN5A channels expressed in tsA201 and HEK-293 cells (Smits et al., 2005; Gui et al., 2010a,b) had revealed that the E161K mutation induces a 2.5- to 4-fold decrease in I_Na peak current density, both in H558 and R558 background. To assess the extent to which this reduction is due to the changes in the I_Na steady-state activation and inactivation curves per se, we reconstructed current traces in response to voltage clamp steps from a holding potential of -120 mV to test potentials ranging from -100 to +50 mV. Figure 3A shows examples of wild-type and mutant current traces at test potentials of -40, -20, and 0 mV. These already demonstrate that the peak current density of the wild-type channels is higher than that of mutant channels, despite the identical value of g_Na used in the voltage clamp simulations. Figure 3B summarizes the simulation data. The changes in the steady-state activation and inactivation curves per se reduce the peak current density by ≈35%. To arrive at an almost 4-fold decrease in I_Na peak current density, in line with the experimental observations, we reduced the mutant g_Na by a factor of 2 in our further simulations. Thus we split the original g_Na of 22.3 nS of the Fabbri-Severi model cell into 11.15 nS for the wild-type I_Na channels and 5.575 nS for the E161K mutant I_Na channels, in either background.

We applied the mutation-induced changes in I_Na, i.e., the shifts in steady-state activation and inactivation curves, the reduction in steepness of the steady-state activation curve, and the reduction in fully-activated conductance, to half of the I_Na channels in the Fabbri-Severi model to assess the effects of the heterozygous E161K mutation on the spontaneous pacemaker activity of a human SA nodal cell. Figures 4A–G, shows the effects on the action potential (Figure 4A), intracellular calcium concentration (Figure 4B) and underlying membrane currents (Figures 4C–G) under control conditions (no ACh), whereas the corresponding simulation data in the presence of 10 nM ACh are shown in Figures 4H–N. The simulated administration of ACh to the Fabbri-Severi model cell does not only activate the acetylcholine-activated potassium current I_K,ACh, which is zero under control conditions, but also reduces the hyperpolarization-activated “funny current” I_f (through a -4.95 mV shift in its voltage dependence of activation), the L-type calcium current I_CaL (through a 3.1% decrease in its maximal conductance), and the rate of Ca²⁺ uptake into the sarcoplasmic reticulum (SR) by the SERCA pump (through a 7.0% decrease in its maximum activity). Of note, the changes in I_CaL and Ca²⁺ uptake rate have only marginal effects at this level of ACh (Fabbri et al., 2017).

The main effect of the application of ACh is a prolongation of the basic cycle length from 813 to 1027 ms (Figures 4A,H, wild-type traces), which is associated with a decrease in the intracellular calcium concentration [Ca²⁺]_i (Figures 4B,I, wild-type traces), and sodium-calcium exchange current I_NCX during diastolic depolarization (Figures 4E,L, wild-type traces). I_f is also reduced (Figures 4F,M), as a result of the direct effect of ACh on its voltage dependence of activation. The remaining inward currents, i.e., the L-type calcium current I_CaL (Figures 4D,K), T-type calcium current I_CaT (Figures 4E,L) and I_Na (Figures 4G,N), are not largely affected. The total repolarizing current I_repol is also hardly affected, with an almost constant value of ≈9 pA during diastolic depolarization (Figures 4D,K). This is because it includes the ACh-activated I_K,ACh, which is also shown separately (Figures 4F,M), in addition to the rapid delayed rectifier potassium current I_Kr, the slow delayed rectifier potassium current I_Ks, the ultrarapid delayed rectifier potassium current I_Kur, the transient outward potassium current I_to, and the sodium-potassium pump current I_NaK, which are not shown separately in Figure 4. The reduction in the “pacemaker currents” I_f and I_NCX (Lakatta and DiFrancesco, 2009), in combination with the minor changes in other inward currents as well as I_repol, result in the reduction in the net membrane current I_net (Figures 4C,J) that underlies the observed increase in cycle length.

In absence of ACh, the mutation-induced reduction in I_Na (Figure 4G) causes a reduction in I_net, which becomes most prominent during the second half of diastolic depolarization (Figure 4C) and in turn results in a prolongation of the cycle length by 53 and 54 ms in H558 and R558 background, respectively (Figure 4A, horizontal arrow), from the basic cycle length of 813 ms, corresponding with a 6% decrease in beating rate from 74 to 69 beats/min in either background. Simulations with the wild-type model in which g_Na is reduced by 50% (labeled “50% g_Na”), thus simply blocking 50% of the channels, result in a cycle length prolongation of 55 ms and traces that are barely distinguishable from the mutant ones, thus demonstrating that the contribution of the mutant component of I_Na to total I_Na is almost zero.

In the simulated presence of 10 nM ACh, the effects of the mutation are more pronounced, as reflected by the prolongation of the cycle length by 104 ms in H558 background and 105 ms in R558 background (Figure 4H, horizontal arrow), from the basic cycle length of 1027 ms, corresponding with a 9% decrease in beating rate from 58 to 53 beats/min in either background. Again, the contribution of the mutant component of I_Na to total I_Na is almost zero, as demonstrated by the highly similar prolongation of cycle length, by 107 ms, in case of 50% g_Na. The more pronounced effect of the highly similar reduction in I_Na in the simulated presence of 10 nM ACh (Figures 4G,N) can be explained by its occurrence in the setting of a smaller I_net during diastolic depolarization (Figures 4C,J).

As illustrated in Figure 5A, the mutation-induced increase in cycle length increases with increasing levels of ACh. At 25 nM, this increase amounts to 347 and 355 ms in H558 and R558 background, respectively, from a basic cycle length of 1448 ms, slightly less than the increase of 361 ms in case of 50% g_Na. In terms of beating rate (Figure 5B), an ACh level of 25 nM results in a rate of 41 beats/min in case of wild-type sodium channels and 33 beats/min in case of the heterozygous E161K mutation, in either background. The percent decrease in beating rate relative to wild type is shown in Figure 5C. In absence of ACh (Figure 5C, leftmost bars), the E161K mutation results in a 6.1% decrease in beating rate in H558 background and 6.2% decrease in R558 background. This percentage increases with increasing levels of ACh. In the simulated presence of 10 nM ACh, the decrease in beating rate amounts to 9.2 and 9.3%, respectively. At an ACh level of 25 nM (Figure 5C, rightmost bars), the mutation-induced decrease in beating rate is almost 20%, with values of 19.3% in H558 background and 19.7% in R558 background.

Sinus node dysfunction may also be related to changes in atrial excitability, potentially leading to sinus node exit block or atrial conduction block. Therefore, we assessed the effects of the E161K mutation on atrial excitability using the Maleckar et al. human atrial cell model (Maleckar et al., 2009), which we used in combination with its well-validated equations for I_K,ACh (Maleckar et al., 2008). Changes in I_Na, simulating a heterozygous mutation in SCN5A, were implemented as in the Fabbri-Severi model. Of note, the effect of ACh on the Maleckar et al. model cell is limited to the activation of I_K,ACh.

Figure 6A shows atrial action potentials elicited at 1 Hz with a stimulus current of 1 ms duration under control conditions (no ACh). Both variants of the E161K mutation result in a reduction in the action potential overshoot and a less rapid activation, indicative of a decrease in capacitive current of the atrial cell (Figure 6A, inset). The associated I_Na, which is shown in Figure 6B, is approximately halved or – because the slower activation allows a larger fraction of the channels to inactivate during the activation process – even more than halved. In the simulated presence of 10 nM ACh, the resting membrane potential becomes more negative with a value of -79 mV vs. -74 mV (Figure 6C), which results in a larger I_Na because less sodium channels are inactivated at this more negative resting potential (Figure 6D). This does, however, not imply that action potential generation is facilitated in presence of ACh since the distance to action potential threshold is increased and an additional outward current, i.e., I_K,ACh, is activated. Again, activation is slowed and the underlying I_Na is approximately halved (Figures 6C,D, insets).

Atrial excitability was characterized, over a wide range of concentrations of ACh, by determining the threshold stimulus current of the atrial cell as well as its maximum upstroke velocity, as a direct measure of the capacitive current. Threshold stimulus current was determined with 1 ms stimuli of increasing amplitude that were applied at a frequency of 1 Hz. Results are shown in Figure 7A. The threshold increases with increasing levels of ACh and is 5–8% higher in case of the E161K mutation. Although barely visible in Figure 7A, the threshold in R558 background, with the smaller shift of the steady-state activation curve, is smaller than in H558 background.

Maximum upstroke velocity was determined by eliciting action potentials with a 1 ms, ≈50% suprathreshold stimulus current at a frequency of 1 Hz. As shown in Figure 7B, it is not strongly dependent on the level of ACh. Upstroke velocity is approximately halved, or even more than halved, in case of the E161K mutation. As for the threshold stimulus current, the mutation in R558 background is associated with (slightly) less strong effects than in H558 background. In contrast with the SA nodal cell, the mutant component of I_Na is not effectively zero, as revealed by the larger I_Na (Figures 6B,D, insets) and larger maximum upstroke velocity (Figure 7B) in comparison with the 50% g_Na simulations.

Our simulations demonstrate that the E161K mutation in SCN5A hampers both impulse generation and impulse propagation through its effects on the electrophysiological properties of human SA nodal and atrial cells. Generation of the SA nodal action potential is hampered by the strong decrease in I_Na during diastolic depolarization and the associated decrease in I_net, in particular in presence of ACh (Figures 4, 5). Impulse propagation, as comprehensively reviewed by Kléber and Rudy (2004), is hampered by the current-to-load mismatch that results from, on the one hand, the increase in threshold stimulus current and, on the other hand, the decrease in capacitive current of the atrial cell (Figures 6, 7). The highly common H558R polymorphism in SCN5A does not have a major effect on the outcome of the simulations, despite its effects on the biophysical properties of the E161K mutant channels (Figures 2, 3).

The E161K mutation has been subject of simulations since its discovery in 2005 (Smits et al., 2005). Smits et al. (2005) showed that the E161K mutation impairs conduction in linear strands of atrial and ventricular cells, using the human atrial and ventricular cell models by Courtemanche et al. (1998) and Priebe and Beuckelmann (1998), respectively. Simulations were based on the observations by Smits et al. (2005) on E161K channels expressed in tsA201 cells. Mutation effects were relatively mild in comparison to later observations by Gui et al. (2010a,b) and obtained in H558 background. Smits et al. (2005) also carried out simulations on SA nodal cells, using the model of a peripheral rabbit SA nodal cell by Zhang et al. (2000), including equations for I_K,ACh (Zhang et al., 2002). Despite the largely different time course of the human SA nodal action potential, with a much smaller I_Na, a much longer diastolic phase, and a significantly less negative maximum diastolic potential, the main findings of the present study are similar to those of Smits et al. (2005) regarding the SA nodal action potential: the E161K mutation causes a decrease in diastolic depolarization rate that results in an increase in cycle length, which is much more pronounced in presence of ACh.

Butters et al. (2010) studied the effects of several mutations in SCN5A, including E161K, on cardiac pacemaking in a two-dimensional model of sino-atrial tissue based on a reconstruction of a single slice of the rabbit right atrium, using modified versions of the single cell models of rabbit SA nodal and atrial cells by Zhang et al. (2000, 2002) and Lindblad et al. (1996), respectively. Their simulations of the E161K mutation were again based on the observations by Smits et al. (2005). The mutation slowed down pacemaking and this effect was more pronounced in the presence of ACh. Furthermore, a conduction block in the direction toward the atrial septum occurred in the absence of ACh, while conduction in the direction toward the crista terminalis was sustained. With 15 nM ACh, exit block occurred in both directions and the SA node was unable to drive the surrounding atrium.

In the present study, we carried out simulations with the use of comprehensive well-validated human single cell models and implemented effects of the E161K mutation based on the experimental data by Gui et al. (2010a,b). Furthermore, we assessed the role of the highly common H558R polymorphism, also based on the experimental data by Gui et al. (2010a,b). We focused on cellular mechanisms and refrained from using two- or even three-dimensional models of sino-atrial interaction, with their coupling interactions between heterogeneous cells, because results of such simulations are critically dependent on the exact geometry and cell distribution (Garny et al., 2003a), which have not been fully elucidated.

Our simulations show that the effects of the E161K mutation are so large, both in H558 and R558 background, that the contribution of mutant I_Na to total I_Na of the SA nodal cell is effectively zero. Accordingly, it is not surprising that E161K mutation carriers show sinus node dysfunction as do carriers of several mutations in SCN5A that are associated with a complete or almost complete loss of function of the mutant channels, as reviewed by Lei et al. (2007, 2008) and Remme et al. (2008). In many cases, the sinus node dysfunction shows a reduced penetrance, as for the E161K mutation, where 8 out of 14 mutation carriers show the disease. Reduced penetrance is a common finding in primary electrical diseases, including those associated with mutations in SCN5A (Remme, 2013; Robyns et al., 2014). Several factors, including co-inherited genetic variants and alternative splicing sites, may play a role. However, our simulations suggest that it is highly unlikely that the common H558R polymorphism can explain the reduced penetrance of the E161K mutation.

Unfortunately, experimental data on the effects of the E161K mutation are limited to voltage clamp data from mutant SCN5A channels in mammalian expression systems, viz. tsA201 and HEK-293 cells, obtained at room temperature and in the presence (Smits et al., 2005) or absence (Gui et al., 2010a,b) of a wild-type β-subunit. There are some quantitative differences in the experimental data, potentially due to differences in cell type or presence of β-subunit, with a more severe pattern in the data of Gui et al. (2010a,b), with a +19 mV shift of the steady-state activation curve in H558 background vs. the +11.9 mV shift reported by Smits et al. (2005) and a 4-fold vs. 2.5-fold reduction in peak current density. However, the latter difference may partly be explained by the larger shift in the steady-state activation curve. In our simulations, we have assumed that the experimentally observed effects of the mutation in mammalian expression systems at room temperature also hold for cardiac cells at the physiological temperature of 37°C. We cannot rule out that recordings at a more close-to-physiological temperature would have revealed additional effects of the E161K mutation, as has been the case for several other BrS-related mutations in SCN5A expressed in mammalian cell lines, like T1620M and Y1795H (Dumaine et al., 1999; Rivolta et al., 2001). Furthermore, we have to keep in mind that the cardiac sodium channel is part of a macromolecular complex that does not only comprise the α- and β-subunits, but also several other proteins that may regulate channel activity, such as ankyrin, caveolin and syntrophin, and may act differently in case of mutant channels (Ruan et al., 2009). Thus, it is uncertain to which extent the experimental data on the effects of the E161K mutation acquired in mammalian expression systems represent the mutant-induced changes that occur in human cardiac cells.

In our simulations, it is assumed that the Fabbri-Severi model of a human SA nodal cell (Fabbri et al., 2017) and the Maleckar et al. model of a human atrial cell (Maleckar et al., 2009) are fully representative of the biophysical properties of the corresponding real cells. However, experimental data on the biophysical properties of cardiac cells, and human cells in particular, are often incomplete or obtained under non-physiological conditions. Functional data on I_Na in human SA nodal cells, for example, are limited to accidental observations by Verkerk et al. (2009) when performing voltage clamp experiments on single human SA nodal pacemaker cells to record I_f. Accordingly, I_Na was included in the Fabbri-Severi model, but simply adopted from its parent model, i.e., the model of a rabbit SA nodal cell by Severi et al. (2012), who, in turn, used the equations of the SA nodal cell model by Noble et al. (1989). The latter equations were adopted from the model of SA nodal cell electrical activity by Noble and Noble (1984), who based their equations on the model of cardiac electrical activity by DiFrancesco and Noble (1985). The description of I_Na in the latter model is essentially that of Hodgkin and Huxley (1952), but with modifications based on data from experiments on rabbit cardiac Purkinje fibers at 10–26°C (Colatsky, 1980) and from isolated rat ventricular myocytes at 20–22°C (Brown et al., 1981). Activation and inactivation curves were shifted along the voltage axis to account for their observed temperature dependence, whereas rate constants were scaled to arrive at time constant values near the experimentally observed ones. The I_Na equations in the Maleckar et al. (2009) model are adopted from the model of an adult human atrial cell by Nygren et al. (1998), who used the voltage clamp data of Sakakibara et al. (1992) on I_Na in isolated human atrial myocytes at 17 ± 1°C to construct their model equations. Activation and inactivation curves were shifted in the positive direction along the voltage axis to reach a realistic activation threshold and stable resting potential. The time constants of activation and inactivation are very similar to those of the rabbit atrial cell model of Lindblad et al. (1996), who based their mathematical description of I_Na on the voltage clamp data of Wendt et al. (1992), acquired from cultured rabbit atrial myocytes at 17°C, and I_Na of the DiFrancesco-Noble model (DiFrancesco and Noble, 1985). Activation kinetics of I_Na were modified from those of the DiFrancesco-Noble model with the use of a Q₁₀ of 1.7.

In our simulations, it is also assumed that mutation effects are limited to half of the intrinsic I_Na channels, thus representing the heterozygous nature of the mutation, taking into account that a functional I_Na channel is built with a single Na_V1.5 protein. This assumption is challenged by recent experimental findings of Clatot et al. (2018), who demonstrated that wild-type and mutant Na_V1.5 proteins can form dimers, enabling coupled gating of wild-type and mutant Na_V1.5 sodium channels that can be responsible for a dominant negative effect of the mutation. Thus, the cellular effects of the E161K mutation could be more severe than anticipated.

The action of ACh in the Fabbri-Severi model is, apart from the activation of I_K,ACh, limited to suppression of I_f, I_CaL, and Ca²⁺ uptake into the SR, thus simplifying the intracellular SR-based calcium clock signaling cascade through adenylyl cyclase, cyclic AMP and protein kinase A (Maltsev et al., 2014; Behar et al., 2016). This cascade is important for a detailed understanding of the action of ACh. Yet, in the present study, where the focus is on the effects of the E161K mutation in SCN5A and possible effects of the H558R polymorphism, it is sufficient that the suppression by ACh of I_f, I_CaL, and Ca²⁺ uptake into the SR, and that of I_NCX through the reduced intracellular calcium level, are included, in addition to the activation of I_K,ACh.

One should be careful in the interpretation of experimental data, not only because of the above considerations, but also because changes in biophysical parameters are not necessarily independent. This is nicely demonstrated by the decrease in peak current density shown in Figure 3. Such decrease is often interpreted as a reduction in the number of channels or a reduction in their single conductance. In this case, however, the E161K mutant-induced changes in steady-state activation and inactivation curves per se reduce the peak current density by ≈35%. This effect may, at least partly, explain the apparent discrepancy between the ≈30% decrease in cell surface expression of E161K/H558 channels in HEK-293 cells and the > 70% decrease in peak current density reported by Gui et al. (2010a). Similarly, an apparently slower inactivation of E161K/H558 channels at membrane potentials ranging from -40 to 0 mV can, at least partly, be explained by the +19 mV shift in their steady-state activation curve (Gui et al., 2010a). Because experimentally observed changes in time constants of (fast) inactivation or recovery from inactivation, if any, are already small, no attempts were made to include these changes in the simulations of the present study.