Single bio-pacemaker cell model

Previous experimental studies [10,11,20–22] implemented the suppression of Kir2.1, the incorporation of HCN channels and the expression of TBX18 to induce pacemaking in VMs. In this study, we used a human VMs model built by Ten Tusscher et al. in 2006 (TP06 model) [52] as the basal model to investigate possible pacemaking mechanisms in VM-transformed pacemaking cells. In brief, the basal VM cell model can be described by the following ordinary differential equation: (1) where V is the voltage across cell membrane surfaces, t is time, I_ion is the sum of all transmembrane ionic currents, and C_m cell capacitance.

The I_ion in the original ventricular model is described by the following equation: (2) where I_Na is fast sodium channel current, I_K1 is inward rectifier potassium channel current, I_to is transient outward current, I_Kr is rapid delayed rectifier potassium channel current, I_Ks is slow delayed rectifier potassium channel current, I_CaL is L-type calcium current, I_NaCa is Na⁺/Ca²⁺ exchange current, I_NaK is Na⁺/K⁺ pump current, I_pCa and I_pK are plateau Ca²⁺ and K⁺ currents, and I_bCa and I_bNa are background Ca²⁺ and Na⁺ currents. The formulations and their parameters for the ionic channels of human VM cells were listed in [52,54].

To mimic the reduction of Kir2.1 expression [11,20,21] or the suppression of I_K1 by expressing TBX18 [10], in simulations I_K1 was decreased by modulating its macroscopic channel conductance (G_K1). To mimic the incorporation of I_f in VMs experimentally [22], we modified the basal model of Eq 2 by incorporating human SAN I_f formulation [53]. In simulations, I_f was modulated by changing its channel conductance (G_f).

As a result, the I_ion for the bio-pacemaker model can be described as: (3) where I_K1 could be expressed by (4) where G_K1 is the conductance of I_K1, x_k1∞ is a time-independent inward rectification factor, K_o is extracellular K⁺ concentration and E_K is the equilibrium potentials of K⁺. S_K1 is a scaling factor used to simulate the change of I_K1 expression level.

The I_f channel is permeable to Na⁺ and K⁺ ions [53]. As a result, I_f could be described by two components: (5) (6) (7) where G_f,Na and G_f,K are maximal I_f,Na and I_f,K channel conductance, y is a time-independent inward rectification factor that is a function of voltage, E_Na, E_K are equilibrium potentials of Na⁺ and K⁺ channels respectively, and S_f is a scaling factor used to simulate the change of I_f expression level.

Formulations of other channel currents for the VM cell model are the same as those in the original model in [52]. The bio-pacemaker cell model can be found in S1 Code.

I_CaT was incorporated into the I_K1/I_f-induced bio-pacemaker cell model to investigate its role in generating pacemaking action potentials. Its formulations can be found in [50]. As a result, the I_ion is described as: (8)

Evaluating criterion of the pacemaking stability and ability

To analyse the effect of I_K1 and I_f on pacemaking activity, we simulated the membrane potential under different current densities of I_K1 and I_f, with I_K1 being reduced systematically by from between 60–100% (i.e., I_K1 density at -80 mV changed from 0.396 to 0 pA/pF while the I_K1 density in the original basal model is 0.99 pA/pF at -80 mV in I-V curve). The representative I-V relation curve under different inhibition of I_K1 is shown in S1A Fig. I_f density was increased by from 0 to 10 folds with a basal value of -0.63 pA/pF at -80 mV in I-V curve (i.e., I_f density changed from 0 to -6.3 pA/pF at -80 mV). The representative I-V relation curve under different incorporation of I_f is shown in S1B Fig.

Two characteristics were used to quantify the state of membrane potentials generated by the ventricular pacemaker model: the continuity and validity of spontaneous APs. The continuity was used to quantify whether or not the automaticity of membrane potential could sustain with time; whilst the validity was used to characterize whether every automatic wave was biologically-valid or not. As such, we defined the following:

W: a valid wave. An action potential whose wave trough was less than -20 mV and wave crest was more than 20 mV could be considered as a valid wave.

αW, α<1: an incomplete wave.

R: a resting period lasting 1000 ms.

Wⁿ: the concatenation of n W’s.

(W R): the concatenation of W and R.

As such, the pacemaking behaviour of the I_K1/I_f pacemaker model can be categorized as FIVE states whose definitions were as followings:

None pacemaking behaviour during the entire simulation period could be described as State-1: (9)

Transient spontaneous pacemaking behaviour could be described as State-2: (10)

Bursting pacemaking behaviour could be described as State-3: (11)

Persistent pacemaking activity with periodically incomplete depolarization could be described as State-4: (12)

Stable pacemaking activity could be described as State-5: (13)

With regard to the pacemaking ability, when the pacemaking behaviour was stable, the cycle length (CL) under varied I_K1 and I_f was calculated. The CL was defined as the averaged wavelength of pacemaking activity over a period of simulation of more than 400 s, ensuring the accuracy of the computed CL. As the basal model was for VMs, a long simulation period was necessitated to achieve a completely stable pacemaking status and minimize the effect of the transition period.

Characteristics of pacemaking during diastolic interval

The length of the diastolic interval (DI) is an important measure to characterize the pacemaking ability. In this study, we defined that DI as the time interval between the time of MDP (S2A Fig, t₁) and the time when the membrane potential reaches at -55 mV (i.e., around the activation potential of I_CaL) (S2A Fig, t₂). The diastolic upstroke velocity during DI was defined as the change rate of the membrane potential, taking the following formulation: (14)

The unit of diastolic upstroke velocity was V/s.

We evaluated the contributions of all of the inward currents to membrane potential during DI and found that the main inward currents which helped to depolarize membrane potential during DI are I_Na, I_NaCa and I_f. Their contribution can be described by an average integral during DI: (15)

Similarly, the main outward currents which held membrane potential at diastolic potential during DI are I_K1, I_NaK, I_Kr and I_Ks, the integral of which can be described as: (16)

Dynamic analysis in I_K1/I_f parameter space

Dynamic pacemaker AP behaviour was dependent on the balance of I_K1 and I_f interactions. Simulations were conducted in the I_K1/I_f parameter space to characterize this dependence. Results are shown in Fig 1. With differing combinations of I_K1 and I_f density, five different regions for distinctive pacemaking dynamics could be discerned, including stable pacemaking activity (blue area, State-5 defined by Eq 13), intermittence of failed depolarization (yellow area, State-4 defined by Eq 12), bursting pacemaking behaviour (orange area, State-3 defined by Eq 11), transient pacemaking activity (green area, State-2 defined by Eq 10), and no automaticity (grey area, State-1 defined by Eq 9). In each category, representative membrane potentials are illustrated at the bottom panel of Fig 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Dynamic behaviours of pacemaking action potentials in I_K1/I_f parameter space.

Blue: stable pacemaking activity (State-5); Yellow: persistent pacemaking activity with periodic incomplete depolarization (State-4). Orange: bursting pacemaking behaviour (State-3). Green: transient spontaneous pacemaking behaviour (State-2). Gray: no spontaneous pacemaking behaviour (State-1). In each category, the typical pacemaking action potentials are illustrated at the bottom panel.

https://doi.org/10.1371/journal.pcbi.1008177.g001

With a fixed I_f density, alterations to I_K1 could produce different types of pacemaking activities, and this also applied when I_K1 was fixed whilst I_f was changed. When I_f density was fixed at a density between -0.63 and -2.52 pA/pF, with a 60–80% block of I_K1 (i.e., I_K1 density at -80 mV was in the range of 0.198–0.396 pA/pF), pacemaking activity was generated but with self-termination (Fig 1, green area). Then, a further reduction in I_K1 or a slight increase in I_f induced bursting pacemaking behaviour, as shown by the orange area in Fig 1, which was between the boundaries marking the persistent automaticity and transient pacemaking activity regions. A further increase in I_f or suppression in I_K1 could produce persistent automaticity (Fig 1, yellow and blue area). But when I_K1 was greater than about 0.248 pA/pF, incomplete depolarization appeared periodically (Fig 1, yellow area). Finally, a stable and spontaneous pacemaking activity could be generated when I_K1 was decreased to less than 0.248 pA/pF at -80 mV with I_f included (Fig 1, blue area).

To illustrate possible reasons responsible for the effect of each changed ion current density on modulating pacemaking states, we chose four typical cases to analyze the corresponding pacemaking mechanisms in different pacemaking states in the I_K1/I_f parameter space. Cases are listed in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Typical cases and corresponding state.

https://doi.org/10.1371/journal.pcbi.1008177.t001

Initiation of transient spontaneous depolarization

In the basal VM cell model with the suppression of I_K1 by 70% (the density of I_K1 at -80 mV was 0.297 pA/pF), incorporation of I_f (with a current density of -0.63 pA/pF) was unable to depolarize the membrane potential and lead to spontaneous pacemaking activity because the excessive outward current of I_K1 counteracted the inward depolarizing current. This state can be described by State-1 as shown in Eq 9. When the density of I_f was increased to -1.89 pA/pF, spontaneous depolarization was provoked at the beginning of the transition period, however, the automaticity self-terminated after 163 s (Fig 2A), showing a State-2 behaviour as described in Eq 10.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Transient spontaneous pacemaking behaviour.

(A-F) Membrane potential (V), intracellular Na⁺ concentration ([Na⁺]_i), intracellular Ca²⁺ concentration ([Ca²⁺]_i), fast sodium current (I_Na), Na⁺/K⁺ pumping current (I_NaK) and L-type calcium channel current (I_CaL) with the current densities of (I_K1, I_f) at (0.297 pA/pF, -1.89 pA/pF) (‘CASE 1’) during the entire simulating period of 800 s. (G) Expanded plots of (A-F) for the time course of pacemaking self-termination (160–164 s). H: Maximum diastolic potential of spontaneous pacemaking behaviour of (A).

https://doi.org/10.1371/journal.pcbi.1008177.g002

We analysed possible ion channel mechanisms responsible for unstable and self-terminating pacemaking APs with the current densities of (I_K1, I_f) at (0.297pA/pF, -1.89 pA/pF) (defined as ‘CASE 1’). Results in Fig 2 showed that during the time course of the spontaneous pacemaking, there were changes of intracellular ionic concentrations and the MDP. There was an accumulation of the intracellular Ca²⁺ concentration ([Ca²⁺]_i, Fig 2C) during the time course of spontaneous pacemaking APs. Such an accumulation of [Ca²⁺]_i was because the automaticity in VMs shortened the DI between two successive APs, leaving insufficient time for Ca²⁺ in the cytoplasm to be extruded to restore to its initial value after each cycle of excitation. This consequentially led to overload in [Ca²⁺]_i, which suppressed the extent of the activation degree of the L-type calcium current (I_CaL, Fig 2F), especially during the phase 0 of the pacemaking action potential. Furthermore, the overloaded [Ca²⁺]_i increased the I_NaCa (S3A Fig) gradually with time, resulting in an elevated MDP (Fig 2H) that inhibited the activation degree of the fast sodium channel current (I_Na, Fig 2D). Through the Na⁺ permeability of I_f, extra Na⁺ flowed into the cytoplasm [53] during each of the APs. The increase of I_NaCa also accelerated the accumulation of the intracellular Na⁺ concentration ([Na⁺]_i) to increase from 7.67 to 11.8 mM (Fig 2B). The increased [Na⁺]_i augmented the feedback mechanism of Na⁺/K⁺ pumping activity, by which the Na⁺/K⁺ pump current (I_NaK) increased gradually with time (Fig 2E). All of these factors worked together, inhibiting the membrane potential to reach the take-off potential, leading to self-terminated automaticity at 163 s (Fig 2G).

It was also possible to generate automaticity in the model by fixing the current density of I_f at a low value, but with a further reduction in I_K1 density. S4 Fig shows the results when I_f was held at -0.63 pA/pF and the current density of I_K1 was reduced to 0.178 pA/pF. In this case, pacemaking activity appeared in the model, but the automaticity was unstable and self-terminated due to similar mechanisms as shown in Fig 2 for the increased-I_f situation.

Bursting pacemaking behaviour

Fig 3 shows the intermittent bursting behaviour, which is generated with a different combination of I_K1 and I_f current densities in the model. In the figure, the current densities of (I_K1, I_f) were at (0.297 pA/pF, -2.52 pA/pF) (defined as ‘CASE 2’). Such kind of pacemaking state can be classified as State-3 (Eq 11).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Bursting pacemaking behaviour.

(A-F) Membrane potential (V), intracellular Na⁺ concentration ([Na⁺]_i), intracellular Ca²⁺ concentration ([Ca²⁺]_i), fast sodium channel current (I_Na), Na⁺/K⁺ pumping current (I_NaK) and L-type calcium channel current (I_CaL) with the current densities of (I_K1, I_f) at (0.297 pA/pF, -2.52 pA/pF) (‘CASE 2’) during the entire simulating period of 800 s. (G) Expanded plots of (A-F) for the time course of pacemaking resumption (383–385 s). H: Maximum diastolic potential of automatic pacemaking activity when I_f is -2.52 and -1.89 pA/pF (‘CASE 2’ vs. ‘CASE 1’, solid and dotted line respectively).

https://doi.org/10.1371/journal.pcbi.1008177.g003

In ‘CASE 2’, the spontaneous oscillation was unstable, characterized by self-termination and then resumption after a quiescent period (Fig 3A). During the first pacemaking stage (simulation period from 0 to 236 s in Fig 3), the self-termination was accompanied by the overload of [Na⁺]_i, the accumulation of [Ca²⁺]_i, which caused the reduction of I_Na, the increase of I_NaK and the decrease of I_CaL (Fig 3B–3F). An increased [Na⁺]_i, accumulated [Ca²⁺]_i, reduced I_Na, increased I_NaK and decreased I_CaL can also be observed in the pacemaking stage in Fig 2. This suggested that the underlying mechanisms responsible for the self-termination of the pacemaking APs were similar to those of ‘CASE 1’.

It is of interest to analyse the mechanism(s) for the resumption of the pacemaking APs after a long pause. It was shown that, during the time course of the quiescent interval (236–384 s in Fig 3A–3F), the intracellular Na⁺ (Fig 3B) continued to be extruded out of the cell by I_NaK (Fig 3E). As such, the [Na⁺]_i (Fig 3B) gradually decreased over time. The decrease in [Na⁺]_i led to a gradually reduced I_NaK over the time course of quiescence (Fig 3E), which decreased the suppressive effect of I_NaK on depolarization. Via I_NaCa (S3B Fig), the intracellular Ca²⁺ was kept to be extruded out of the cell, leading to a decreased [Ca²⁺]_i (Fig 3C). The reduction in [Ca²⁺]_i reduced the Ca-induced inactivation gate of I_CaL, resulting in an increased I_CaL. Ca²⁺ concentration in the sarcoplasmic reticulum ([Ca²⁺]_SR) also decreased via a leakage Ca release (I_leak) from SR to cytoplasm (S8C–S8D Fig) because of the reduced [Ca²⁺]_i. Moreover, as compared with ‘CASE 1’, the increase in I_f helped to produce a more depolarized MDP (Fig 3H), allowing the membrane potential more easily to reach the take-off potential for initiation of the upstroke. During the quiescent stage, I_CaL and I_Na kept at a small magnitude (Fig 3D and 3F) because of the quiescent membrane potential at which little I_CaL and I_Na were activated. Until the Ca²⁺ oscillation produced a full course of action potential with a sufficiently depolarised MDP, I_f and I_Na was fully activated, facilitating the resumption of the spontaneous pacemaking activity at 385 s (Fig 3G). This process of self-termination and resumption repeated alternately, which constituted bursting behaviour.

Persistent pacemaking activity

A further increase in I_f ((I_K1, I_f) at (0.297 pA/pF, -3.15 pA/pF)) produced a series of persistent spontaneous APs. Results are shown in Fig 4 (grey lines) for APs (Fig 4A), together with a phase portrait of membrane potential (V) and total membrane channel current (I_total) (Fig 4B), I_K1 (Fig 4C), and I_CaL (Fig 4D). Although the spontaneous APs were sustained during the entire simulation period of 800 s, there were some incomplete depolarizations observed periodically (Fig 4A, grey line) which can be classified as State-4 according to Eq 12. This pacemaking situation was termed as ‘CASE 3’.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Persistent pacemaking activity.

The current densities of (I_K1, I_f) of stable pacemaking activity and periodically incomplete pacemaking activity are at (0.248 pA/pF, -3.15 pA/pF) (‘CASE 4’, black lines) and (0.297 pA/pF, -3.15 pA/pF) (‘CASE 3’, grey lines) respectively. (A-D) The membrane potential (V), phase portraits of membrane potential against the total membrane channel current (I_total), inward rectifier potassium channel current (I_K1) and L-type calcium channel current (I_CaL) during the simulating time course of 760–763 s. (Inset) Expanded plot during V from -75 to -45 mV and I_total from -0.1 to 0 pA/pF marked by the asterisk in (B).

https://doi.org/10.1371/journal.pcbi.1008177.g004

When the density of I_K1 was further reduced to 0.248 pA/pF (Fig 4C, black line), a stable pacemaking activity was established (Fig 4A, black line), with an average CL of 895 ms and MDP of -72.63 mV. This kind of pacemaking state can be described as State-5 by Eq 13. In this condition, the pacemaking activity was robust and the pacing CL was close to that of the native human SAN cells (approximately 800–1000 ms [55]). We termed stable pacemaking activity with (I_K1, I_f) at (0.248 pA/pF, -3.15 pA/pF) as ‘CASE 4’.

In order to understand potential mechanism(s) underlying the genesis of incomplete pacemaking potentials in CASE 3, phase portraits of membrane potential against I_total for incomplete (grey line) and complete (black line) depolarization APs were plotted and superimposed for comparison, as shown in Fig 4B, with a highlight of phase portraits during the diastolic pacemaking potential range from -75 mV to about -45 mV being shown in the inset. In the case of incomplete depolarization (CASE 3), there was a greater I_K1 (Fig 4C) that counteracted the inward depolarizing current, leading to a smaller I_total during the diastolic depolarization phase (see the grey line in Fig 4B and the inset). Consequentially the membrane potential failed to reach the take-off potential for the activation of I_CaL (Fig 4D, grey line), leading to an incomplete course of action potential (Fig 4A, grey line). In CASE 4, with a reduced I_K1 (Fig 4C, black line), there was a greater I_total during the diastolic depolarization phase (Fig 4B and the inset, black line), which drove the membrane potential to reach the take-off potential for the activation of I_CaL (Fig 4D, black line), leading to the upstroke of the action potential.

The frequency for the appearance of the incomplete AP was dependent on the density of I_K1. Incomplete depolarization occurred less frequently, with progressively smaller I_K1. By way of illustration, the incomplete depolarization appeared once every three cycles with the I_K1 density at 0.297 pA/pF in CASE 3 (S5A Fig), but this became once every five cycles when I_K1 density was reduced to 0. 277 pA/pF I_K1 (S5B Fig). This suggested that a large residual of I_K1 in the VM model might result in the failure of complete depolarization.

Pacemaking cycle length in I_K1/I_f parameter space

A systematic analysis of the relationship between the calculated CL in the I_K1 and I_f density parameter space is presented in Fig 5. In the figure, the measured CL was coloured from 650 ms in dark red to 1000 ms in yellow. In this study, we regarded persistent pacemaking action potential with CL 1000 ms or less as ‘valid pacemaking activity’ (corresponding to ‘State-5’ with CL < = 1000 ms in Fig 1), therefore, only the CLs of the valid pacemaking potentials are shown in Fig 5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Measured cycle length in I_K1/I_f parameter space.

The density of I_K1 is from 0 to 0.396 pA/pF and the density of I_f is from 0 to -6.3 pA/pF at -80 mV. The measured cycle length (CL) is coloured from 650 ms in dark red to 1000 ms in yellow. White means that there is no definitely measured pacemaking CL (i.e., the pacemaking action potential is not persistent or stable) or the pacemaking activity is ‘invalid’ (CL > 1000 ms).

https://doi.org/10.1371/journal.pcbi.1008177.g005

It was shown that a sufficient depression in I_K1 (up to 75%; I_K1 density < 0.248 pA/pF) was required to produce a stable pacemaking action potential with a ‘valid’ pacemaking frequency. With the increase of the I_K1 inhibition level, the CL became shortened at all I_f densities considered. Also, the more that I_K1 was inhibited, the less I_f was needed to provoke ‘valid’ spontaneous pacemaking activity. By contrast, the effect of I_f on the CL presented two phases, which was dependent on the I_K1 density. When I_K1 was less than 0.198 pA/pF, the pacemaking ability became robust with the increase in I_f density. However, when I_K1 was increased from 0.198 to 0.248 pA/pF, an increase in I_f actually slowed the pacemaking activity, leading to an increased CL.

Reciprocal role of I_f and I_K1 in generating pacemaking APs

Further analysis was conducted to investigate the reciprocal role of reduced I_K1 and increased I_f in generating pacemaking APs. This section analysed the role of I_f and I_K1 in the parameter space which produced stable pacemaking activity (the ‘State-5’ in Fig 1) with the density of I_K1 at 0.05 and 0.198 pA/pF. By sufficiently reducing I_K1 to a density of 0.05 pA/pF alone (in the absence of I_f), the model was able to generate stable spontaneous APs with a CL of 1011 ms (S6A and S6E Fig, dotted lines). In this case, the incorporation of I_f with a small density helped to boost the pacemaking activity and increase the pacemaking frequency. It was shown that with the incorporation of I_f at a density of -0.63 pA/pF, the CL was reduced by 233 ms, changing from 1011 ms to 778 ms (S6A Fig). Compared with the case of I_f absence, incorporation of I_f—even with a small density—helped to depolarize cell membrane potential during the early DI phase (S6A and S6C Fig). Moreover, a similar phenomenon to Fig 2 was that the incorporation of I_f led to an accumulation of [Na⁺]_i (S6B Fig) and [Ca²⁺]_i (S6D Fig and Inset C). The accumulation of [Ca²⁺]_i increased I_NaCa (S6F Fig) which contributed to the depolarization of membrane potential. As a result, the incorporation of I_f facilitated genesis of spontaneous APs and shortened the DI significantly thus decreased the CL.

With a fixed I_K1 density of 0.05 pA/pF, the relationship between the computed CL of pacemaking APs and I_f density was shown in Fig 6A. In a range from 0 to -2.52 pA/pF, an increase in I_f density produced a marked decrease in the CL (Fig 6A), which was associated with an increase in the rate of membrane depolarization during the DI (diastolic depolarizing rate) (S2B Fig). In this range, an increase in I_f caused an accumulation of [Na⁺]_i (Fig 6C), which enhanced I_NaCa (Fig 6D). Particularly, there was a precipitous decline of I_Na when I_f was in the range from 0 to -1.26 pA/pF (Fig 6F). This was because the augmented I_f elevated the MDP of the action potential (Fig 6B), which affected the activation level of I_Na. Furthermore, when I_f density was over -2.52 pA/pF, there was a less dramatic decrease in CL with an increase of I_f (Fig 6A). This was attributable to a reduced I_Na (Fig 6F) as a consequence of the gradual elevation of the MDP (Fig 6B). Another factor was that the maximum density of I_f was not increased linearly (Fig 6E) because of elevated MDP.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Effect of I_f density on the pacemaking cycle length under different I_K1 density.

(A) Change of cycle length with the increase of I_f from 0 to -6.3 pA/pF when I_K1 density is 0.05 (solid line) and 0.198 pA/pF (dotted line). (B-F) Change of maximum diastolic potential (MDP), maximum intracellular Na⁺ concentration ([Na⁺]_i), maximum Na⁺/Ca²⁺ exchange current (I_NaCa), maximum funny current (I_f) and maximum fast sodium current (I_Na) during a pacemaking period with the increase of I_f from 0 to -6.3 pA/pF. (G-H) Change of the normalized total integral of main inward currents (Integral I_in) and normalized total integral of main outward currents (Integral I_out) during DI phase with the increase of I_f when I_K1 density is 0.05 (solid line) and 0.198 pA/pF (dotted line). The inward currents include I_Na, I_NaCa and I_f; the outward currents include inward rectifier potassium channel current (I_K1), Na⁺/K⁺ pumping current (I_NaK), rapid delayed rectifier potassium channel current (I_Kr) and slow delayed rectifier potassium channel current (I_Ks).

https://doi.org/10.1371/journal.pcbi.1008177.g006

Depending on the I_K1 density, the relationship between CL and I_f density could also be bi-phasic. With a small I_K1 density (0.05 pA/pF), the measured CL decreased monotonically with the increase in I_f density (Fig 6A, solid line). However, with a large I_K1 (0.198 pA/pF), the measured CL first decreased with an increased I_f density, but then increased with it when I_f density was greater than -3.15 pA/pF, implicating a slowdown in the pacemaking activity with the increase of I_f (Fig 6A, dotted line). The action potentials with the current densities of (I_K1, I_f) at (0.198 pA/pF, -3.78 pA/pF) and (0.198 pA/pF, -5.04 pA/pF) indicated that such a slowdown of pacemaking APs with an increased I_f was mainly due to the prolonged DI (S1 Text). The extra I_f accelerated pacemaking activity, leaving insufficient time for intracellular Ca²⁺ to be extruded, which resulted in a higher amplitude of [Ca²⁺]_i (Fig A in S1 Text). As a consequence, I_NaCa amplitude increased, leading to an increased [Na⁺]_i, and finally an increased outward I_NaK (Fig B in S1 Text). At the same time, I_f elevated the MDP (Table A in S1 Text), so the amplitude of I_Na decreased (Fig B in S1 Text). The changes in these currents caused a prolonged CL at a great I_f density. This indicated that the bi-phasic effect of I_f may be attributable to the delicate balance between inward currents and outward currents during the DI phase. To extend this hypothesis to the whole I_f parameter space, I_in and I_out were defined as the total integral of three dominant inward currents during the DI phase (I_Na, I_NaCa and I_f) (Eq 15) and the total integral of four dominant outward currents during the DI phase (I_K1, I_NaK, I_Kr and I_Ks) (Eq 16). At a low I_K1 density of 0.05 pA/pF, with the increase in I_f density, there was a monotonic decrease of CL (Fig 6A, solid line). However, at a large I_K1 density (e.g., 0.198 pA/pF), more than -3.15 pA/pF I_f resulted in an increased CL (Fig 6A, dotted line). This observation can be explained by the unsaturated I_in. The density of I_in was an integrated consequence of MDP, the activation degree of I_Na and I_f, as well as the density of I_NaCa. At 0.05 pA/pF I_K1 density, I_out was small, therefore I_Na, I_NaCa and I_f was saturated (i.e., I_in was saturated). However, when I_K1 was large, increased I_out led to a more negative MDP which resulted in an increased I_in due to more activated I_Na, I_NaCa and I_f. But I_in reached at a asymptotic value after I_f density was greater than -3.15 pA/pF (Fig 6G, dotted line) because the a similar level of the diastolic [Ca²⁺]_i (i.e., minimum [Ca²⁺]_i) at different I_K1 caused a decreased normalized integral of I_NaCa (Fig C in S1 Text), therefore I_in was not saturated at large I_K1 (Fig 6G, dotted line). Consequentially, the I_out outbalanced the I_in, leading to a prolonged DI that increased the CL at large I_K1.

Incorporation of I_CaT in I_K1/I_f pacemaker model

The incorporation of I_CaT in the I_K1/I_f pacemaker model with the current densities of (I_K1, I_f) at (0.297 pA/pF, -0.63 pA/pF) induced automaticity in non-pacemaking cell model (S7A Fig). It was shown that I_CaT initiated the oscillation of Ca²⁺ in the cytoplasm and the SR (S7B and S7C Fig) in the non-pacemaking cell model, activating the I_NaCa and I_CaL (S7D and S7E Fig).

However, when I_CaT was incorporated into a stable I_K1/I_f pacemaker model with the current densities of (I_K1, I_f) at (0.099 pA/pF, -0.63 pA/pF), an increased pacemaking CL was observed (853 ms vs. 950 ms, S8A Fig). The reason was that the addition of I_CaT in the model affected the Ca²⁺ dynamics, resulting in an accumulation of [Ca²⁺]_i as well as [Ca²⁺]_SR, leading to an increased I_NaCa and greater inactivation of I_CaL (S8C–S8F Fig). An increased I_NaCa consequently caused an accumulation of [Na⁺]_i, leading to an increased outward I_NaK (S8G Fig). More importantly, an increased I_NaCa elevated the MDP, leading to a less activated I_Na and I_f (S8H and S8I Fig). These together with an increased I_NaK produced a prolonged CL.

Summary of major findings

Biological experiments provided possibilities to transform non-pacemaking CMs or stem cells into bio-pacemaker cells by regulating the expression of I_f [13,28,29,36,37] or I_K1 [11,20,21]. Computationally, bifurcation theory was used to analyse the effect of the density of an individual ion channel current on the membrane potential of a cardiac cell [47,50], by which decisive currents in pacemaking initiation can be screened out. The interaction between multiple ion channels in cardiac pacemaking has also been considered in some prior studies. Lakatta et al. [51] investigated the effect of multiple component ion channels on cardiac pacemaking by identifying numerically minimal ensembles of ion channels in the SAN model. Their study provided a simplification of the model which may be suitable for further bifurcation analysis. However, due to the simplification, the model does not reflect the whole ionic dynamics of cardiac cells. In another study, Kurata et al. [50] simulated the combined action of I_f and I_CaT, I_st or I_CaL with three specific values of I_CaT, I_st or I_CaL considered. That study also suggested that there was a threshold of I_K1 below which automaticity can be induced in I_f-incorporated ventricular model, but potential underlying ionic mechanisms on how balanced I_K1/I_f modulates the stability of the pacemaking activity was not shown. Considering the biological possibility that the superiority of combined action of I_K1 and I_f [10,32,44] and unclear interaction between I_K1 and I_f in the initiation of stable pacemaking activity in a whole cardiac cell model, the present study developed a virtual bio-engineered pacemaking cell model based on a human ventricular model by downregulating of I_K1 and incorporating I_f. In comparison with previous computational modelling studies on bio-pacemakers [50,51], the present study provides new insights into understandings of 1) the dynamic regulation of ionic concentrations and ionic channel currents in I_K1/I_f-induced ventricular pacemaker; 2) the interaction between I_K1 and I_f in the genesis of pacemaking action potentials; and 3) the stability of pacemaking potentials in wide I_K1/I_f parameter space.

The present I_K1/I_f pacemaker model can initiate robust and stable pacemaking activity by balancing the actions of reduced I_K1 and increased I_f, though the effect of each manipulation on pacemaking activity is different. While the action of a reduced I_K1 on the pacemaking activity is mono-phasic, that of an increased I_f is biphasic. It was shown that inhibiting I_K1 promotes pacemaking ability and stability. The incorporation of I_f at an appropriate level promotes pacemaking activity, but an excessive I_f might result in abnormal pacemaking activity accompanied by abnormal intracellular ionic concentrations (such as pacemaking activity with periodically incomplete depolarization) or weak pacemaking ability, which could be proarrhythmic. As a result, the reciprocal interaction between I_K1 and I_f is crucial for producing stable spontaneous pacemaking activity in VMs. In the present model, specific I_f density at different I_K1 densities or vice versa was suggested.

This study demonstrated the feasibility of creating VM-based pacemaker cells by downregulating I_K1 and upregulating I_f and provided evidence of the superiority of I_K1/I_f-based pacemaker model theoretically. The results of this study may be useful for optimizing the future design of engineered bio-pacemakers.

Role of I_K1 suppression on pacemaking activity

Simulation results that suppressing I_K1 by 75% - 100% can initiate stable pacemaking behaviour are in consistence with those of experimental findings, where it has been found that more than 80% inhibition of I_K1 was required to produce a pacemaking phenomenon in guinea-pig’s ventricular cavity [11,21]. It is also in agreement with previous bifurcation analyses in showing that it required I_K1 to be reduced to at least 15% of the control value to transform a VM cell model to be auto-rhythmic [47]. And a complete block of I_K1 produced a spontaneous pacemaking activity with a CL of 795 ms [47], close to 833 ms when I_K1 was completely suppressed in the present study. The effect of I_K1 suppression on pacemaking activity was monotonic in our pacemaker model. The more the I_K1 is blocked, the faster the pacemaking activity is with all I_f densities considered. Similar results have also been observed in another human ventricular cell model [56] based on modifications of the O’Hara and Rudy model (ORD model) [57] (Fig A in S2 Text, solid and dotted lines).

Though our simulation results suggest an important role of sufficient suppression of I_K1 in generating persistent and stable pacemaking APs, it is noteworthy that the deficiency of I_K1 has been reported to be lethal for adult rodents [58]; and loss function of Kir2 gene may prolong QT intervals as well as cause Andersen’s syndrome [59]. Consequently, suppression of I_K1 from VM for generating a biological pacemaker may only be suitable when applied to highly localized, designated ‘pacemaker’ regions.

Role of I_f on pacemaking activity

I_f has been shown to play an important role in generating pacemaking APs in both native [13,22,28,29,36,37] and engineered pacemakers [43]. Experimentally it has been shown that high expression of HCN2 can initiate spontaneous beats in neonatal rat VMs [22,36] and improve spontaneous beats in rabbit CMs [13]. HCN4 incorporation by the expression of TBX18 can also initiate spontaneous pacemaking activity in both rodent VMs [10] and porcine VMs [12].

In the present study, I_f helped to promote pacemaking activity, via its action of depolarization during the diastolic depolarization phase as well as its action on the intracellular ion concentrations. The inclusion of I_f in the VM cell model caused the accumulation of [Na⁺]_i by Na⁺ channel of I_f [53]. The enhanced pacemaking activity caused by extra I_f also induced the accumulation [Ca²⁺]_i because there was not enough time to extrude Ca²⁺ from the cytoplasm [60]. The accumulated [Ca²⁺]_i increased I_NaCa, which promoted membrane potential depolarization especially during the early stage of DI (S6 Fig). Such a promoting action of I_f in bio-pacemaking can also be seen in another independent model as shown in S2 Text.

The increase in I_f density can enhance the automaticity in most cases. However, the effect of I_f on the pacemaking activity was observed to be bi-phasic. When it was increased to be over a threshold, excessive I_f resulted in an elevated MDP (Fig 6), which caused a reduced activation of I_f and I_Na, leading to a slowdown of the ability of pacemaking activity. This phenomenon occurred when the I_K1 was not suppressed sufficiently. The impairing effect of excessive I_f on pacemaking APs was also observed in another ventricular pacemaker model based on modification of the ORD model [56] (S2 Text). It was shown that a greater increase in I_f density even terminated pacemaking activity (Fig B in S2 Text). The possible impairing effect of I_f on I_Na was verified by the fact that in the bio-pacemaker induced by HCN2 expression [61], co-expression of the skeletal muscle sodium channel 1 (SkM1), in order to hyperpolarize the action potential threshold, helped to counterbalance the negative effect of I_f overexpression, producing an accelerated depolarization phase. In fact, in the original TP06 model, the peak amplitude of I_Na was about -300 pA/pF at the resting potential of -86.2 mV [54]. In our pacemaker model, the peak amplitude of I_Na was significantly reduced because of the elevated MDP. This suggested that to counterbalance the elevated MPD and the reduced I_Na, an increase in the channel expression of I_Na may help to produce an enhanced pacemaker. This might be simulated by increasing I_Na conductance in the model study. Furthermore, when I_K1 was large, an increase in I_f even lengthened pacemaking period or caused unstable pacemaking behaviour (Fig 1). This simulation result is in agreement with a previous biological experimental study that observed arrhythmicity when acute HCN gene was expressed [42]. Another experimental study showed that HCN2-expression caused an excessive increase in the basal beating rate [62]. In our model, it has been also observed that excessive I_f may cause an overly fast pacemaking rate.

Reciprocal interaction between I_K1 and I_f

Our study demonstrates that the reciprocal interaction between I_K1 and I_f plays a crucial role in creating stable and persistent pacemaking. Only an optimal combination of I_K1 and I_f can initiate stable pacemaking activity. In the present pacemaker model, the greater the degree of I_K1 suppression, the smaller was the I_f density required for the generation of spontaneous oscillation (Fig 1). And modulation of the two currents simultaneously helps to create a physiologically-like pacemaker that is better than that produced by manipulating I_K1 or I_f alone (Fig 5). Such observation of reciprocal interaction between I_K1 and I_f in pacemaking is consistent with previous experimental observations. Previous studies have shown that although suppressing I_K1 [11,21], or incorporating sufficient I_f [22] alone was able to initiate pacemaking activity in VM cells, a pacemaker constructed by TBX18 showed greater stability, due to its combined actions of I_K1 reduction and I_f increase [10]. Another experiment in porcine VMs [12] also indicated that TBX18 expression did not increase the risk of arrhythmia, which means that a mixed-current approach is probably a superior means of producing a bio-pacemaker. Experiments in a Kir2.1/HCN2 HEK293 cell [45] and Kir2.1/HCN4 [43] showed that I_K1 may actually recruit more I_f by activating current at more negative membrane voltages because I_K1 was the only hyperpolarizing current in these experiments. Our simulations, however, did not yield such a result because the interaction of other outward currents (such as I_NaK, I_Kr and I_Ks) contributed to the hyperpolarization of membrane potential and helped the activation of I_f. We thought an integrated action between all of ionic currents in cardiac cells should be considered, rather than evaluated specific ionic currents in a partial model.

In addition, simulation results indicated that I_K1 expression level may influence the I_f’s effect on the pacemaking activity (Fig 6A). Excessive I_K1 hindered I_f’s ability to modulate pacemaking activity. This further showed that the balanced expression of I_K1 and I_f affected the balance between the inward and outward currents during the diastolic depolarization phase, thus affected the membrane potential state and the pacemaking CL of the pacemaker. An experiment showed a coincident result that the expression of HCN2 in adult rat VMs could not cause spontaneous beats due to the high expression of I_K1 [22], but in neonatal rats, the I_K1 was less so that expressing HCN2 could provoke automaticity. Similarly, such a dynamic balance between the inward and outward currents during the repolarization and the diastolic depolarization phase was also affected by other repolarization currents, such as I_Kr, I_Ks and I_NaK etc. Possible effects of modulating these repolarization currents on the bio-pacemaking warrants further studies in the future.

The present I_K1/I_f-induced pacemaker model exhibited greater robustness than I_K1-based or I_f-based pacemaker models. In the I_K1-based model, the range of I_K1 density that could initiate spontaneous beatings was from 0 to 0.0246 pA/pF, while in I_K1/I_f-based pacemaker model, this value extended to 0.248 pA/pF (Fig 1). The superiority of I_K1/I_f-based pacemaker model than I_f-based pacemaker model seemed to be more distinct. Incorporating I_f alone at a high density of -6.3 pA/pF could not provoke any spontaneous beating, but combining with the suppression of I_K1, small incorporation of I_f helped to ignite automaticity (Fig 1). The flexibility of this system also reflected in the easy modulation of CL via manipulating I_K1 and I_f density (Fig 5).

Compared with the human SAN cell model developed by Fabbri et al. [53] and human SAN cell [64], the action potential generated by the I_K1/I_f-induced pacemaker model had a longer action potential duration at 90% (APD90) and a more negative MDP when the CL was similar (see Table 2). Such differences may be attributable to the fact that there are regional differences in the functional expression of ionic currents between the SAN and ventricular myocytes [63]. In the presented I_K1/I_f-induced pacemaker model, though we have reduced I_K1 and incorporated I_f to a similar level of ion channel current densities as the SAN, other ionic currents in the present I_K1/I_f-induced pacemaker model inherited the same channel properties of VMs cell model, causing different pacemaker behaviour in the I_K1/I_f biological pacemaker model as compared to the SAN model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Comparison of pacemaker behaviours between the SAN cell and the bio-pacemaker cell.

https://doi.org/10.1371/journal.pcbi.1008177.t002

Ca²⁺ dynamics in I_K1/I_f pacemaker model

There is still debate about the relative role of two pacemaking mechanisms of membrane clock (I_f) and Ca²⁺ clock [65]. A biological experiment demonstrated that Ca²⁺-stimulated adenylyl cyclase AC1 can promote pacemaking ability in HCN2-expressed left bundle branches [62]. A model study [66] that evaluated the synergism between Ca²⁺ clock and membrane clock in SAN central cell, also showed that the synergistic system was more robust and flexible. Another study [67] showed that VMs may also have Ca²⁺ clock, which provided a probability for the creation of Ca²⁺ clock-based bio-pacemaker. The role of Ca²⁺ dynamics in bio-pacemaker was also shown in our I_K1/I_f pacemaker model. As shown in Fig 3G, the resumption of pacemaking activity in bursting behaviour was provoked by the oscillation of [Ca²⁺]_i. This indicated that the Ca²⁺ dynamics played an important role in the creation of bio-pacemaker, which warrants further study.

Furthermore, considering the role of I_CaT in the genesis of pacemaking APs in native SAN cells, a theoretical investigation of potential role(s) of I_CaT in the bio-pacemaker was conducted using the I_K1/I_f-modulated pacemaker model. It was shown the effect of I_CaT had dual aspects. On one hand, the incorporation of I_CaT might promote the pacemaking ability of ventricular pacemaker [50]. by initiating Ca²⁺ oscillation thus producing spontaneous beatings in quiet pacemaker model. On the other hand, the incorporation of I_CaT might affect the MDP, leading to secondary actions on the homeostasis of ion concentrations, as well as ion channel currents including I_Na, I_f, I_NaCa and I_NaK, which slowed down the pacemaking activity.

Limitations

Limitations of the human VMs model we used in this study has been described elsewhere [54]. In this study, the I_f formulation of human SAN [53] was incorporated into the original VMs model. The properties of I_f, including the conductance of I_f, the half-maximal activation voltage (V_1/2) and time constants of the activation, may present species-dependence. In the present version, we only consider the conductance of I_f but have not discussed other properties of I_f. Moreover, in this study, we only investigated the pacemaking action potential at the single-cell level, without considering the intercellular electrical coupling between pacemaker cells as presented in the SAN tissue. These limitations are now being addressed for future versions of the model. In addition, bio-pacemaker models developed from other cardiac cell types, such as atrial myocytes, warrant future studies. Additionally, the other pacemaking-related currents in native SAN cells, such as I_Na and I_st, could also be adjusted for creating stable bio-pacemaker.

One of possible advantages of bio-pacemaker over the traditional electronic pacemaker is at its possible sensitivity to autonomic regulation. It is of interest to study how the pacemaking action potentials are modulated by autonomic regulation by β-Adrenergic receptor stimulation or cholinergic receptor stimulation [10], which warrants further future investigation.

It is necessary to highlight these limitations, they nevertheless do not affect our conclusions on the underlying pacemaking mechanisms of engineered bio-pacemaker cells, especially regarding the reciprocal interaction of I_K1 and I_f for a robust bio-pacemaker in modified VMs.