Reciprocal interaction between IK1 and If in biological pacemakers: A simulation study

Pacemaking dysfunction (PD) may result in heart rhythm disorders, syncope or even death. Current treatment of PD using implanted electronic pacemakers has some limitations, such as finite battery life and the risk of repeated surgery. As such, the biological pacemaker has been proposed as a potential alternative to the electronic pacemaker for PD treatment. Experimentally and computationally, it has been shown that bio-engineered pacemaker cells can be generated from non-rhythmic ventricular myocytes (VMs) by knocking out genes related to the inward rectifier potassium channel current (IK1) or by overexpressing hyperpolarization-activated cyclic nucleotide gated channel genes responsible for the “funny” current (If). However, it is unclear if a bio-engineered pacemaker based on the modification of IK1- and If-related channels simultaneously would enhance the ability and stability of bio-engineered pacemaking action potentials. In this study, the possible mechanism(s) responsible for VMs to generate spontaneous pacemaking activity by regulating IK1 and If density were investigated by a computational approach. Our results showed that there was a reciprocal interaction between IK1 and If in ventricular pacemaker model. The effect of IK1 depression on generating ventricular pacemaker was mono-phasic while that of If augmentation was bi-phasic. A moderate increase of If promoted pacemaking activity but excessive increase of If resulted in a slowdown in the pacemaking rate and even an unstable pacemaking state. The dedicated interplay between IK1 and If in generating stable pacemaking and dysrhythmias was evaluated. Finally, a theoretical analysis in the IK1/If parameter space for generating pacemaking action potentials in different states was provided. In conclusion, to the best of our knowledge, this study provides a wide theoretical insight into understandings for generating stable and robust pacemaker cells from non-pacemaking VMs by the interplay of IK1 and If, which may be helpful in designing engineered biological pacemakers for application purposes.

Introduction Currently, electronic pacemaker implantation is the only non-pharmacological therapy for some patients with pacemaking dysfunction. But electronic pacemakers have some possible limitations [1]. Implantation of a pacemaker device may have complications for patients, especially for aged ones because of their infirm health [2]. Pediatric patients can receive electronic pacemakers; however, the device has to be replaced as they grow, and repeated surgeries are needed [3]. Electronic devices can be subject to electromagnetic interference [4], which causes inconvenience to the patients. A further issue is that classical electronic pacemakers are insensitive not only to hormone stimulation [5] but also to autonomic emotion responsiveness [4], although there are some attempts to make them respond to autonomic nervous control [6]. In addition, the long-term use of electronic pacemakers has been reported to increase the risk of heart failure [7]. Appropriately designed biological pacemakers (bio-pacemakers) have the potential to overcome some of the limitations of electrical device use [8], such as the lack of pacing flexibility. Also, engineered bio-pacemakers could potentially involve only minor surgical trauma for implantation as well as facilitating chronotropic responses [9]. In previous experimental studies, it has been shown that a bio-pacemaker can be engineered via adenoviral gene transduction [10][11][12] or lentiviral vector [13,14] techniques, by which non-pacemaking cardiac myocytes (CMs) can be transformed to the rhythmic pacemaker-like cells.
The native cardiac primary pacemaker, sinoatrial node (SAN), is a special region comprised of cells with distinct electrophysiological properties to those of cells in the working myocardium. Such intrinsic and special electrophysiological properties of SAN cells are mainly manifested by their small if not absent inward rectifier potassium channel current (I K1 ) [15], but a large "funny" current (I f ) [16] that is almost absent in atrial and ventricular cells. I K1 helps to hyperpolarize cell membrane potential and contributes to maintaining a stable negative resting potential in ventricular myocytes (VMs). But in the SAN, absence of I K1 makes a more positive maximum diastolic potential (MDP) than the resting potential in VMs, which helps its depolarization. I f plays an important role in phase 4 depolarization of SAN action potentials, thus influences the pacemaking rate in the SAN [17]. In addition to the absence of I K1 and presence of I f , T-type Ca 2+ channel current (I CaT ) [18] and sustained inward current (I st ) [19] also contribute to spontaneous pacemaking activity in SAN cells. Similar to I f , I CaT plays a role in the depolarization in the SAN but is usually undetectable in VMs. Such unique electrophysiological properties of SAN cells formed a theoretical basis to engineer non-pacemaking CMs into spontaneous pacemaker cells. These non-pacemaker cells include native CMs, such as ventricular [11,[20][21][22], atrial [23] or bundle branch myocytes [24]. They can also be stem cells, such as embryonic stem cells [25][26][27], bone marrow stem cells [13,28,29], adipose-derived stem cells [30][31][32], or induced pluripotent stem cell [14,33,34].
With gene therapy, these non-rhythmic cells have been manipulated to provoke automaticity. In previous studies, knocking down the Kir2.1 gene to reduce the expression of I K1 promoted spontaneous rhythms in newborn murine VMs [20]; by reprogramming the Kir2.1 gene in guinea-pigs, VMs also produced pacemaker activity when I K1 was suppressed [11,21]. As I f plays an important role in the native SAN cell pacemaking, a parallel gene therapy manipulation to create engineered bio-pacemaker has been carried out by expressing the HCN gene family in non-rhythmic CMs [35]. It has been shown that expressing HCN2 produced escape beats in canine CMs [23,24] and initiated spontaneous beats in neonatal rat VMs [22]. HCN expression in stem cells-induced-CMs also enhanced their pacemaking rate [13,28,29,36,37]. Overexpressing HCN4 can also induce spontaneous pacemaking activity in mouse embryonic stem cells [38]. However, acute HCN gene expression might have a side effect on the normal cardiac pacemaking activity [39][40][41]. The overexpression of the HCN gene in non-rhythmic CMs can cause ectopic pacemaker automaticity and even arrhythmicity [42].
It has been suggested that a combined manipulation of I K1 and I f may be a better alternative for creating a bio-pacemaker [43]. The expression of the transcriptional regulator TBX18, which influenced both I f and I K1 expression, generated appropriate automatic responses in non-pacemaking CMs [10,32,44]. In addition, reprogramming TBX18 in porcine VMs did not show the increase of arrhythmia risk [12], indicating the probable superiority of manipulating I K1 and I f jointly for generating a bio-pacemaker. Furthermore, an experimental study showed that in HEK293 cells with expression of the Kir2.1 and the HCN genes, not only I f but also moderate I K1 was necessary to induce spontaneous rhythmic oscillations [45].
Computational modelling offers a means to investigate different approaches to generating stable pacemaking activity. Based on the cardiac cell model, bifurcation analysis was widely used to explore the effect of changes on some individual ion channel current on the pacemaking activities, such as the role of down-regulating I K1 in the pacemaking initiation in VMs [46][47][48] and the role of I f in the pacemaking activity in SAN cells [49]. The interplay between multiple ion channel currents in the pacemaking initiation was also investigated in VMs [50] or SAN [51] models, with focus on their roles in sustaining the spontaneous oscillation. Up to now, the dynamic interplay between I K1 and I f on modulating the stability of bio-pacemaker action potentials (APs) has not yet been comprehensively investigated in biophysically detailed cardiac single-cell models.
Considering the availability of experimental data on transformation of VMs into rhythmic cells [11,[20][21][22] and the cardiac dysfunctions arising from atrioventricular heart-block, we chose the VM cell model for investigating bio-pacemaking activity. In this study, we constructed a bio-pacemaker model based on a human VMs model [52] by manipulating I K1 and incorporating I f [53] into the model. This study aimed to investigate (i) possible mechanism(s) underlying the pacemaking activity of the VMs in the I K1 /I f parameter space; and (ii) the reciprocal interaction of reduced I K1 and increased I f in generating stable pacemaking APs. In addition, possible factors responsible for impaired pacemaking activity due to the inappropriate ratio between I K1 and I f were also investigated. Moreover, theoretical analysis of the role of I CaT in I K1 /I f pacemaker model was investigated. This study provides insight into generating stable and robust engineered bio-pacemaker.

Single bio-pacemaker cell model
Previous experimental studies [10,11,[20][21][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: 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: 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 2+ exchange current, I NaK is Na + /K + pump current, I pCa and I pK are plateau Ca 2+ and K + currents, and I bCa and I bNa are background Ca 2+ 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: where I K1 could be expressed by where G K1 is the conductance of I K1 , x k11 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: 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:

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 n : 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: Transient spontaneous pacemaking behaviour could be described as State-2: Bursting pacemaking behaviour could be described as State-3: Persistent pacemaking activity with periodically incomplete depolarization could be described as State-4: Stable pacemaking activity could be described as State-5: 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 1 ) and the time when the membrane potential reaches at -55 mV (i.e., around the activation potential of I CaL ) (S2A Fig, t 2 ). The diastolic upstroke velocity during DI was defined as the change rate of the membrane potential, taking the following formulation: 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: 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:

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. 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 selftermination (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.

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.
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 2+ concentration ([Ca 2+ ] i , Fig 2C) during the time course of spontaneous pacemaking APs. Such an accumulation of [Ca 2+ ] i was because the automaticity in VMs shortened the DI between two successive APs, leaving insufficient time for Ca 2+ in the cytoplasm to be extruded to restore to its initial value after each cycle of excitation. This consequentially led to overload in [Ca 2+ ] 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   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.  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 2+ ] 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 2+ ] 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'.

Bursting pacemaking behaviour
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 2+ was kept to be extruded out of the cell, leading to a decreased [Ca 2+ ] i (Fig 3C). The reduction in [Ca 2+ ] i reduced the Ca-induced inactivation gate of I CaL , resulting in an increased I CaL . Ca 2+ concentration in the sarcoplasmic reticulum ([Ca 2+ ] SR ) also decreased via a leakage Ca release (I leak ) from SR to cytoplasm (S8C- S8D Fig) because of the reduced [Ca 2+ ] 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 2+ 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'.
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 (

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. 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) 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 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 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.
Depending on the I K1 density, the relationship between CL and I f density could also be biphasic. 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 Change of maximum diastolic potential (MDP), maximum intracellular Na + concentration ([Na + ] i ), maximum Na + / Ca 2+ 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 (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 2+ to be extruded, which resulted in a higher amplitude of [Ca 2+ ] 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) (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 . 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 monophasic, 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 guineapig'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 2+ ] i because there was not enough time to extrude Ca 2+ from the cytoplasm [60]. The accumulated [Ca 2+ ] 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 biopacemaking 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], coexpression 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.

Ca 2+ 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 2+ clock [65]. A biological experiment demonstrated that Ca 2+ -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 2+ 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 2+ clock, which provided a probability for the creation of Ca 2+ clock-based bio-pacemaker. The role of Ca 2+ 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 2+ ] i . This indicated that the Ca 2+ 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 2+ 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.