Computational simulations of the effects of the G229D KCNQ1 mutation on human atrial fibrillation

Abstract

Atrial fibrillation (AF) is related to mutations at the genetic level. This includes mutations in genes that encode KCNQ1, a subunit of the I _Ks channel. Here, we investigate the mechanism of gain-of-function in I _Ks towards the occurrence of AF. We used the Courtemanche–Ramirez–Nattel (CRN) human atrial cell model (Am J Physiol Heart Circ Physiol 275:H301–H321, 1998) and applied the modification proposed by Hasegawa et al. (Heart Rhythm 11:67–75, 2014) to fit the behavior of I _Ks due to the G229D mutation in KCNQ1 under a heterozygous mutant form. This was incorporated into two-(2D) and three-dimensional (3D) tissue models, where the mutation sustained a reentrant wave. However, under the wild-type condition, the reentrant wave terminated before the end of our simulations (in 2D, the spiral wave terminated before 10 s, while in 3D, the spiral wave terminated before 13 s). Sustained reentry under the mutation conditions also resulted in a spiral wave breakup in the 3D model, which was sustained until the end of the simulation (20 s), indicating AF.

Introduction

Atrial fibrillation (AF) is the most commonly sustained cardiac arrhythmia, and causes morbidity and mortality. AF is characterized by rapid and irregular activation of the atrium, so that the heart cannot pump blood effectively to the ventricle. Therefore, AF can increase the probability of blood clotting in the atria, and hence the risk of stroke [1]. Several diseases possibly promote AF, including congenital heart disease, congestive heart failure, and hypertension [2]; however, AF may also occur without these diseases. It has been shown that this (also known as idiopathic or lone AF) is influenced by genetic factors [3, 4].

In 2003, Chen et al. [4] reported a correlation between gene mutation and the behavior of I _Ks. The mutation resulted in a change of the amino acid serine to glycine at position 140 (S140G) in transmembrane segment 4 of KCNQ1 resulting in the gain-of-function of I _Ks. Gain-of-function of I _Ks was also observed as the effect of a novel KCNQ1 mutation that was reported by Hasegawa et al. [5] in 2014. This mutation caused an amino acid change from glycine to aspartic acid at position 229 in transmembrane segment 4 of KCNQ1 (see G229D in Fig. 1). To investigate the role of the G229D mutation, Hasegawa et al. [5] carried out electrophysiological experiments and computer simulations using a one-dimensional (1D) model. Hasegawa et al. [5] concluded that the mutation shortened the action potential (AP) duration (APD) in the atria, inducing AF. However, they did not show how the reduction in APD affected the occurrence of AF. As the simulations carried out in that study were in 1D, they could not mimic the reentrant dynamics. However, AF occurs in three-dimensional (3D) spaces such as the atria tissue, and therefore, a 3D simulation should be carried out to investigate AF due to the mutation.

Fig. 1

Schematic representation of the ionic channels, pump, and exchanger in the CRN model. I _b,Ca background Ca²⁺ current, I _Ca,L L-type Ca²⁺ current, I _NaK Na⁺/K⁺ pump current, I _p,Ca sarcolemmal Ca²⁺ pump current, I _NaCa Na⁺/Ca²⁺ exchanger current, I _Na fast inward Na⁺ current, I _b,Na background Na⁺ current, I _Kr rapid component of delayed rectifier K⁺ current, I _Kur ultra-rapid component of delayed rectifier K⁺ current, I _Ks slow component of delayed rectifier K⁺ current, I _to transient-outward K⁺ current, I _K1 inward rectifier K⁺ current. The squares illustrate ionic channels and transporters (only for I _up and I _p,Ca), and circles illustrate ionic exchangers. The empty symbols correspond to channels and Na⁺/Ca⁺ exchanger that do not require ATP, and the filled symbols correspond to ionic transporters and Na⁺/K⁺ pump that require ATP. E1 transmembrane domain in KCNE1, S1–S6 transmembrane segments in KCNQ1, and P P-loop in KCNQ1

Several studies have been carried out to simulate AF phenomena using either two-dimensional (2D) or 3D models [6–8], but not for the G229D mutation. The purpose of this study was to investigate the effects of the G229D mutation on AF using both 2D and 3D computational models of the human atria.

Methods

Human atrial cell model

Many researchers have developed mathematical models of ion transport through cardiac cell membranes. These include human ventricular models [9, 10], and human atrial models [11, 12]. Here, we used the Courtemanche–Ramirez–Nattel (CRN) human atrial cell model [11]. Membrane potential behavior of a single cell was described using following relation:

$$\frac{{{\text{d}}V_{\text{m}} }}{{{\text{d}}t}} = - \frac{{(I_{\text{ion}} + I_{\text{stim}} )}}{{C_{\text{m}} }}$$

(1)

where V _m (mV) is the membrane potential, t the time (ms), I _ion (pA/pF) the total ionic current, and I _stim (pA/pF) is the total stimulus current flowing across the membrane. The total ionic current in the CRN model is given by

$$I_{\text{ion}} = I_{\text{Na}} + I_{{\text{K}}1} + I_{\text{to}} + I_{\text{Kur}} + I_{\text{Kr}} + I_{\text{Ks}} + I_{{\text{Ca}},{\text{L}}} + I_{{\text{p}},{\text{Ca}}} + I_{\text{NaK}} + I_{\text{NaCa}} + I_{{\text{b}},{\text{Na}}} + I_{{\text{b}},{\text{Ca}}} ,$$

(2)

where I _Na is Na⁺ current. The K⁺ currents consist of I _K1, I _to, I _Kur, I _Kr, and I _Ks, which represent contributions of the inward rectifier current, transient outward current, ultra-rapid component of delayed rectifier current, rapid component of delayed rectifier current, and slow component of delayed rectifier current, respectively. The term I _Ca,L represents the L-type Ca²⁺ current, and I _p,Ca is the sarcolemmal Ca²⁺ pump current. The terms I _NaK and I _NaCa are the Na⁺/K⁺ pump current and Na⁺/Ca²⁺ exchanger current, respectively, and I _b,Na and I _b,Ca are the background Na⁺ and Ca²⁺ currents, respectively. This model keeps track of the concentration of various ions inside the cell, and the concentration of ions outside the cell is assumed to be fixed (see Table 1). Figure 1 shows an illustration of various ionic channels in the CRN model.

Table 1 Ionic concentrations outside the cell [11]

G229D KCNQ1 mutation in the I _Ks channel

The CRN model is suitable for the study of reentrant arrhythmia in the human atrium [8, 13]. However, to investigate the effects of the G229D mutation in the human atrial model, a modification to the I _Ks the current model is necessary to fit the experimental data provided by Hasegawa et al. [5] as given by Eqs. (3–5).

$$I_{\text{Ks}} = \overline g_{\text{Ks}} \cdot \left({1 + \frac{0.6}{{1 + \left({\frac{{3.8 \times 10^{ - 5} }}{{[{\text{Ca}}^{2 + } ]_{\text{i}} }}}\right)^{1.4} }}}\right) \cdot x_{{\text{s}}1} \cdot x_{{\text{s}}2} \cdot (V - E_{\text{k}} ),$$

(3)

$$\frac{{{\text{d}}x_{{\text{s}}1} }}{{{\text{d}}t}} = \frac{{(x_{{\text{s}},\infty } - x_{{\text{s}}1} )}}{{\tau_{x_{{\text{s}}1} } }}$$

(4)

and

$$\frac{{{\text{d}}x_{{\text{s}}2} }}{{{\text{d}}t}} = \frac{{(x_{{\text{s}},\infty } - x_{{\text{s}}2} )}}{{\tau_{x_{{\text{s}}2} } }},$$

(5)

where [Ca²⁺]_i (mM) is the intracellular Ca²⁺ concentration, x _s1 and x _s2 are the activation and deactivation gates of I _Ks, respectively, x _s,∞ is the steady-state value of the x _s1 and x _s2 gates, $\tau_{x_{{\text{s}}1} }$ is a time constant for x _s1, and $\tau_{x_{{\text{s}}2} }$ is the time constant of x _s2. We used $\overline g_{\text{Ks}}$ = 0.0136 nS pF⁻¹ for this study. Both the wild-type (WT) and mutation conditions of I _Ks were modeled using Eqs. (3–5). The difference between the WT and mutation conditions was described using x _s1 and x _s2 kinetic formulas. For WT conditions, time-dependent behaviors of x _s1 and x _s2 can be calculated using the following relations:

$$x_{{\text{s}},\infty } = \frac{1}{{1 + \exp ({ - \frac{{V_{\text{m}} + 28.8}}{15.45}})}},$$

(6)

$$\tau_{x_{{\text{s}}1} } = 326.9 + \frac{0.4}{{2.326 \times 10^{ - 4} \cdot \exp ({\frac{{V_{\text{m}} + 65.5}}{17.8}}) + 1.292 \times 10^{ - 3} \cdot \exp ({ - \frac{{V_{\text{m}} + 227.2}}{230}})}}$$

(7)

and

$$\tau_{x_{{\text{s}}2} } = \frac{5}{{0.01 \cdot \exp ({\frac{{V_{\text{m}} - 50}}{100}}) + 0.0193 \cdot \exp({ - \frac{{V_{\text{m}} + 66.54}}{155}})}};$$

(8)

with the mutation, time-dependent behaviors of x _s1 and x _s2 can be calculated using:

$$x_{{\text{s}},\infty } = \frac{0.85}{{1 + \exp ({ - \frac{{V_{\text{m}} + 82.8}}{41.72}})}},$$

(9)

$$\tau_{x_{{\text{s}}1} } = 326.9 + \frac{0.4}{{2.326 \times 10^{ - 4} \cdot \exp ({\frac{{V_{\text{m}} + 119.5}}{17.8}}) + 1.292 \times 10^{ - 3} \cdot \exp ({ - \frac{{V_{\text{m}} + 281.2}}{230}})}}$$

(10)

and

$$\tau_{x_{{\text{s}}2} } = \frac{5}{{0.01 \cdot \exp ({\frac{{V_{\text{m}} - 50}}{100}}) + 0.0193 \cdot \exp ({ - \frac{{V_{\text{m}} + 66.54}}{155}})}}$$

(11)

The I _Ks modification induced changes in the resting state values of some variables. Table 2 shows the initial values used for this work.

Table 2 Initial conditions (state variables of the model at rest)

2D and 3D human atrial tissue models

Cardiac electrophysiological models mimic the electrical conduction phenomena of the AP in cardiac tissue using an electrical conduction equation that is based on continuum mechanics. A partial differential equation for the passive electrical conduction and an ordinary differential equation for the active ionic conduction were coupled to simulate active wave propagation through atrial tissue:

$$\frac{{\partial V_{\text{m}} }}{\partial t} = \frac{1}{{\beta C_{\text{m}} }}(\nabla \cdot \tilde \sigma_{\text{l}} \nabla V_{\text{m}} + \nabla \cdot \tilde \sigma_{\text{e}} \nabla \theta_{\text{e}} ) - \frac{1}{{C_{\text{m}} }}(I_{\text{ion}} (V_{\text{m}} ,\tilde n) - I_{\text{stim}} ),$$

(12)

where V _m is the membrane potential, t the time, β the surface to volume ration, C _m the total membrane capacitance, $\tilde \sigma_{\text{l}}$ the intracellular conductivity, $\tilde \sigma_{\text{e}}$ the extracellular conductivity, θ _e the extracellular potential, $\tilde n$ the gating variables of ionic channels, and I _ion and I _stim are the total ionic current and the stimulus current flowing across the membrane, respectively. Table 3 provides more detail about the units and constants used in this equation.

Table 3 Units and constants for Eq. (12)

Two I _Ks models were used in this study, including the WT and G229D mutation I _Ks models. The CRN cell model was incorporated into the 2D and 3D human atrial models. For the 2D simulations, a time step of Δt = 0.01 ms was used, together with a spatial discretization of Δx = Δy = 0.02 cm, and the simulated tissue was 15 × 15 cm in extent. We previously validated this simulation method using 3D simulation studies of human cardiac tissue [14]. A 3D human atrial model was acquired from the Yonsei Severance Hospital (Sinchon-dong, Seodaemun District, South Korea). The model was discretized using prismatic meshes and represented a human left atrium. The model consisted of 452,140 nodes and 904,276 elements, as shown in Fig. 2b.

Fig. 2

Meshes and the S1–S2 protocol. a Mesh of 2D tissue model. b Mesh of 3D human atrial tissue model acquired from Yonsei Severance Hospital, South Korea. c S1–S2 protocol in 2D. Red square represents the border of the S2 area. This depolarized area will propagate to only one direction since there is conduction block on the other direction (represented by black filled arrow), making a reentry wave. d S1–S2 protocol in 3D. White dashed lines represent the border, which S2 is given (this area will be forced back to its resting potential). The wave will propagate to this area, generating a reentry (color figure online)

Simulation protocol

Single-cell simulations were carried out using a dynamic restitution protocol [15]. The stimuli used were currents with amplitude of 30 pA/pF and a 1-ms duration. Stimuli were applied at an assigned basic cycle length (BCL); we used a BCL of 700 ms. After 20 stimuli were applied, the APDs were obtained from the final two APs. The BCL was then reduced, and the stimuli were again applied 20 times, where the APDs were obtained from the final two APs. The BCL was decremented in steps of 10 ms until it reached 200 ms, and then decremented by 2 ms for BCLs of <200 ms. APD restitution (APDR) curves were obtained using two methods: the first was by plotting the measured APD on the y-axis and BCL on the x-axis, and the second was by plotting the APD on the y-axis and diastolic interval (DI) on the x-axis [15]. Here, we used the APD₉₀ (90 % repolarized from the peak potential).

2D and 3D simulations were carried out according to the S1–S2 protocol [10]. With the 2D simulations, a stimulus was applied to the tissue model three times with the same BCL (S1)—here we use a BCL of 700 ms—in an area located in the most left of the tissue (it was similar to a 15 cm long vertical line), producing a planar wave. A second stimulus (S2) was then delivered with a shorter BCL than the first three stimuli. The S2 stimulus was applied to an area of 7.5 × 7.5 cm with rectangle-like shape located in the left bottom corner of the tissue. If the duration of the S2 BCL is sufficiently long, the resulting wave will propagate in all directions and disappear. However, for some S1–S2 intervals, the electrical wave will only propagate in one direction, since the area in other direction will still be in the refractory period. Once the previously depolarized area is repolarized, the wave will then propagate to that area, and reentry will occur. If the duration of S2 BCL is too short, S2 will fail to create a propagating wave since the tissue in the S2 area will still be in the refractory period. Figure 2c shows the protocol applied in the 2D tissue model.

In the 3D tissue model, the S1–S2 protocol was applied as shown in Fig. 2d. Stimuli were applied in the right upper part of the left atrium around the Bachmann’s bundle is located three times with the same BCL. Here we used BCL of 700 ms as for the S1 BCL. These stimuli will produce waves that propagate in all directions. At the third wave, some areas will be forced back to the resting potential. The depolarization will then propagate to the repolarized tissue area and initialize the occurrence of the reentrant wave.

Results

The measured APD₉₀ was 296 and 225 ms for the WT and mutation, respectively, as shown in Fig. 3a. The APD shortening induced by this mutation was more significant than the APD shortening induced by the V141M mutation [16], and was less significant than the APD shortening induced by V241F [17] and S140G [6] mutations. Under G229D mutation conditions, I _Ks was higher than in WT I _Ks condition, as shown in Fig. 3c. This effect was consistent with the results from the study of Hasegawa et al. [5]. The mutation also induced significant changes in some other currents in addition to I _Ks, including I _Kr (see Fig. 3d), I _K1 (see Fig. 3g), I _Kur (see Fig. 3h), I _b,Ca (see Fig. 3i), and I _b,Na (see Fig. 3j).

Fig. 3

Single-cell simulation results. a AP profiles, b fast inward Na⁺ current (I _Na), c slow component of delayed rectifier K⁺ current (I _Ks), d rapid component of delayed rectifier K⁺ current (I _Kr), e transient outward K⁺ current (I _to), f L-type inward Ca²⁺ current (I _Ca,L), g inward rectifier K⁺ current (I _K1), h ultra-rapid component of delayed rectifier K⁺ current (I _Kur), i background Ca²⁺ current (I _b,Ca), and j background Na⁺ current (I _b,Na)

The APDR curves shown in Fig. 4a and b reveal shorter APDs with the G229D mutation for all BCL and DI ranges. Alternans and steep restitution were not observed under both WT and mutation conditions.

Fig. 4

APDR curves. a APD as a function of DI, b APD as a function of BCL

Figure 5 shows the transmembrane potential map for 2D simulations using an S1–S2 interval of 400 ms. We showed the result using this interval because this interval was in the “R” range of the vulnerability window. A reentrant wave was generated under WT conditions, but terminated after 3 s since the first S1 stimulus was applied (see Fig. 5a). However, the reentrant wave was sustained until the end of the simulation (i.e., 20 s) with the G229D mutation (see Fig. 5b). The wavelengths produced under the mutation condition were much shorter than those produced under the WT condition.

Fig. 5

Spiral wave activity in the 2D human atrial tissue modeled using a WT I _Ks and b the G229D mutation I _Ks. S1 BCL is 700 ms and the S1–S2 interval is 400 ms. AP traces, which are presented on the right side, was recorded from the white point (indicated by an arrow at the first figure in each panel). Time zero was considered as the time when the first S1 was applied

Figure 6 shows grids describing the vulnerability to the re-entry window, which were constructed by investigating the consequences of premature stimulation (i.e., S1–S2 protocol). The S1–S2 interval window that produced reentry shifted toward shorter intervals with the G229D mutation than the WT condition.

Fig. 6

Vulnerability to spiral waves constructed by summarizing the outcomes at various S1–S2 intervals. F indicates “failure to propagate.” This means that the S2 failed to generate the wave. R indicates “reentry generation.” Therefore, S2 triggered the occurrence of a spiral wave. P indicates “normal propagation” which means S2 generates wave that propagates normally along the tissue and disappear

Figure 7 shows transmembrane potential maps for the 3D human atrial tissue model using an S1–S2 interval of 700 ms. A spiral wave was observed and developed into a spiral wave breakup after several seconds even though we applied an S1–S2 interval, which is similar to the normal BCL.

Fig. 7

Spiral wave activity in the 3D human atrial tissue model. a The WT I _Ks and b G229D mutation I _Ks. S1 BCL is 700 ms and S2 BCL is 700 ms. AP traces, which are presented on the right side, were recorded from the white point (indicated by arrow at the first figure in each panel). Time zero was considered as the time when the first S1 was applied

Discussion

To the best of our knowledge, the present study represents the first investigation of the mechanism of AF induced because of the KCNQ1 G229D mutation using a 3D human cardiac electrophysiology model. We observed the effects of the G229D mutation at the cellular level in the KCNQ1 subdomain on the occurrence of AF using a 3D simulation. We used the CRN model of human atrial cells [11] and applied a modification to the I _Ks proposed by Hasegawa et al. [5] to describe the behavior of I _Ks with the G229D mutation. The human atrial cell model was then incorporated into 2D and 3D human atrial models. We compared the electrical patterns and propagation wave under WT and mutation conditions.

The main findings of this work were as follows: (1) in a single cell, the APD was reduced due to the KCNQ1 G229D mutation. The peak I _Kr was also reduced. (2) With the mutation, the reentrant wave was more easily induced. (3) In the 2D and 3D simulations, the pro-arrhythmic effect of the G229D mutation resulted in a shorter wavelength, therefore, once a reentrant wave occurs, it is sustained. (4) The results of the 2D and 3D tissue models differed.

The G229D mutation reduced the APD. In addition, our data showed that the mutation significantly increased the peak value of I _Ks and reduced the peak value of I _Kr. I _Ks, together with I _Kr, are significant in the repolarization and termination of the cardiac AP (in phase 3) [18]. Due to the early repolarization of the cell under mutation condition, I _Kr could not reach its maximum open state to repolarize the cell. This induced a lower current peak than with the other condition.

Alternans is a condition where the amplitudes of the heartbeat vary from beat to beat [19]. In the present study, alternans was not observed. The slope of the APDR curve was also <1. Simulation studies carried out by ten Tusscher and Panfilov [10] confirmed that a steep restitution curve leads to instability of the wave [10]. However, they also determined that the slope of the restitution curve is not the main parameter to determine wave instability [10]. This was consistent with our results. In the present study, steep ADPR curve was not observed. However, the electrical wave instabilities were observed in the 2D and 3D simulations under mutation condition.

The G229D mutation increased the atrial tissue susceptibility to re-entry; therefore, it had a pro-arrhythmogenic effect. Qu et al. [20] argued that the APD is the key determinant of re-entrant arrhythmia. A shorter APD will enable the tissue to be excited at a high rate. This is as shown in the vulnerability window in Fig. 6. Under WT condition, the S2 stimulus failed to propagate in the short S1–S2 interval; however, the S2 stimulus was propagated under the mutation condition and induced the re-entry. In addition, tissue susceptibility to re-entry can be indexed via temporal and spatial vulnerability [6]. The temporal vulnerability of cardiac tissue can be observed through the vulnerability window [6]. The vulnerability window itself is a summarization of the outcomes of the simulation from various S1–S2 intervals. The range for re-entry generation in the vulnerability window under the G229D mutation was slightly increased, indicating increasing temporal vulnerability (Fig. 6).

Meanwhile, the spatial vulnerability of re-entry is related to the size of the tissue required to facilitate the re-entrant wave. This index is closely related to the wavelength produced. The wavelength itself is affected by APD and conduction velocity (CV) since wavelength is equal to the product of the APD and the CV. The G229D mutation significantly reduced the wavelength. The tail of the waveform left the tissue faster with a shorter than a longer wavelength. This enabled the wave from the S2 stimulus to propagate to the other areas of tissue, and it had a conduction block on one side of the propagating area that enabled the formation of the re-entrant wave. Therefore, the G229D KCNQ1 mutation increased the spatial vulnerability of the tissue to re-entrant waves.

The G229D mutation resulted in an increase in I _Ks, leading to persistent atrial tachycardia in the 2D model and AF in the 3D model. There were differences between the 2D and 3D results, especially under the mutation condition. These differences could be due to the different dimensions of the tissue model used.

Several gains of function I _Ks caused by the mutation in the KCNQ1 subunit were already identified. Two of them were S140G and V241F mutations [4, 17]. Both of these mutations had the same effect as the G229D mutation in making the immediately active I _Ks and as a result, abbreviated the APD. Computer simulations have been carried out to investigate the effect of these mutations on cardiac tissue models [6, 8]. Imaniastuti et al. [8] in their study regarding the V241F mutation using computational modeling, found that abbreviated AP caused by the mutation results in a short wavelength in the tissue model. This short wavelength makes the re-entry easier to be accommodated by the tissue model. S1–S2 protocol applied in 2D and 3D tissue models under the mutation generated sustained spiral waves until the end of the simulation. However, no spiral wave breakup was observed in that study, which indicated that the mutation in the computational model will possibly induce atrial flutter, not AF.

A computational study that investigated the effect of the S140G mutation carried out by Kharche et al. [6] also highlighted the effect of gain-of-function I _Ks caused by the mutation that makes the re-entry produced in the tissue model become more persistent compared with the WT condition. However, the S140G mutation only slightly increased the vulnerability window, which means that it only slightly increased the temporal susceptibility of the tissue to the re-entry wave. It is clear that in case of the S140G, the spatial vulnerability, which is closely related to tissue size required to accommodate the re-entry, was predominate here to determine susceptibility toward re-entry generation; this condition was similar to the effect of the G229D mutation.

In summary, our 2D and 3D results demonstrated that the G229D mutation within KCNQ1 stabilized the reentrant wave. In simulations under the WT condition, spiral waves terminated early. Simulations under mutation conditions exhibited chaotic electrical propagation, which is indicative of AF. We conclude that the G229D mutation within KCNQ1 increased the likelihood of AF occurrence.

One potential limitation of this study is that we assumed homogenous cellular electrical properties. Another is that we did not consider the effects of cardiac mechanics on tissue geometry. However, these potential limitations are not expected to influence our conclusions significantly.