Simulation of Propagation of Action Potentials in Cardiac Tissue with an Inhomogeneous Distribution of Extracellular Potassium

Our mind works and our heart beats as long as there is propagation of action potentials in our body. Action potentials are generated by the inward and outward movements of numerous ions, such as Na+, K+, Ca2+, Clacross the cell membrane. These action potentials travel in the form of waves through the various chambers of the heart. To fully understand the heart, a detailed knowledge of action potential propagation is essential.


Introduction
Our mind works and our heart beats as long as there is propagation of action potentials in our body. Action potentials are generated by the inward and outward movements of numerous ions, such as Na + , K + , Ca 2+ , Clacross the cell membrane. These action potentials travel in the form of waves through the various chambers of the heart. To fully understand the heart, a detailed knowledge of action potential propagation is essential.
In unhealthy (ischemic) tissue, the extracellular potassium concentration, [K] 0 , is elevated. High [K] 0 is believed to predispose cardiac tissue to break excitation (Roth and Patel) [1], which in turn plays an important role in re-entry induction and defibrillation. Because fibrillation raises extracellular [K] 0 levels, break excitation may play a more important role in defibrillation than is suggested by simulations and experiments using normal [K] 0 values. In their 2003 paper, Roth and Patel theoretically simulated action potential data in a 1 cm 2 tissue verifying that the action potential may fail to propagate in an ischemic region. The wave propagates for [K] 0 =2 to 12.5mM and at higher potassium concentrations the wave dies out.
Sidorov et al. studied propagation in ischemic tissue experimentally. Their optical mapping experiment in a rabbit heart using [K] 0 values of 4, 6, 8, 10 and 12 mM indicated that spatial heterogeneity of the action potential duration (APD) restitution, created with regional elevation of [K + ] 0 , can lead to action potential (AP) instability, 2:1 block, and reentry induction [2].
The goal of this research is to simulate Sidorov et al. experiment using a mathematical model by incorporating a Luo-Rudy representation of the cell membrane current with the bidomain model of the cardiac tissue, and to examine propagation of potassium from a normal to an ischemic region.

Luo rudy model
Hodgkin and Huxley published a landmark paper that paved the way for subsequent electrophysiological models. They successfully explained the propagation of action potentials in the giant squid axon by representing the membrane resistance as resulting from ion flow through three ion channels: sodium, potassium and leakage [3]. Each pathway was accompanied by a battery that had a potential equal to the Nernst potential for that ion. Variable resistances signifying the various ion channels opening and closing were expressed in terms of their conductance. The probability of opening or closing of these ionic channels was governed by three gating parameters m, n and h following the first order kinematics.

( )
In eq. (1), α i and β i represent the rate of opening and closing of the m, n, and h gates. These rate constants depended on the transmembrane voltage.
Hodgkin and Huxley wrote down the rate of change of membrane potential as ( ) Following the mathematical formulation of Hodgkin and Huxley, in 1991 Luo and Rudy developed a model of the ventricular cardiac action potential. They introduced and analyzed more ion channels then other contemporary physiological models. They replaced the three ionic current terms in Hodgkin and Huxley's model with six other currents. They were a fast inward sodium Na + current I Na , a slow inward calcium Ca ++ current I Si , a time dependent delayed rectifier potassium K + current I k , and time-independent potassium inwardly rectifying current I k1 , the plateau K + current I kp , and a time independent background current I b . The membrane current depended on [K + ] 0 as a parameter in this model. Considering a membrane patch of area 1 cm 2 with capacitance 1 µF/cm 2 equation (2) can be rewritten as

Abstract
We report simulation of experimental data for heterogeneous distribution of extracellular potassium, by developing and applying a mathematical model formulated by incorporating a Luo-Rudy representation of the cell membrane current into the bidomain model of the cardiac tissue. Our model serves to examine propagation of extra cellular potassium from a region of low potassium (normal) to a region of high potassium (ischemic). Starting with computer simulation of unipolar stimulation of cardiac tissue, we reproduced results of the passive bidomain model, and also simulated the space-clamped Luo Rudy model. Finally, we combined the Luo-Rudy model with the bidomain model. The diastolic threshold for anodal stimulation first decreased and then increased with increasing potassium concentration ([K] 0 ), reaching a minimum value at [K] 0 =12mM.
Where, E kp = E k1 1. 7.488 1 exp 5.98 The background current is formulated as where, Then the total time independent current is ( )  (19) The complete listing of all the governing equations, parameters and their corresponding values can be found in Appendix.

Bidomain model
The two or three dimensional cable model, more commonly known as the bidomain model, describes the interaction of electric current with cardiac tissue. It was developed in the late 1970's (Muler and Markin, Miller and Geselowitz, Tung) [9][10][11] and was applied to more problems in cardiac electrophysiology in the 1980s (Plonsey and Barr, Barr and Plonsey, Roth and Wikswo) [12][13][14].
Here we will solve the coupled parabolic and elliptic equations using a finite difference technique (Roth) and discuss some consequences of unipolar stimulation of cardiac tissue [15]. The first analysis of this kind of work was done by Sepulveda et al. [16]. The starting point here is to solve the set of two coupled partial differential equations governing the transmembrane potential, V m , and extracellular potential, V e To solve the parabolic type of equation, we use the simple Euler method. The idea here is to solve eq. (20) for the new value of V m , and then solve the elliptic eq. (21) for the new value of V e using the new value for V m as a source term, by over-relaxation.
By, rearranging, we get, In eq. (22), let, iy ey g g X y All these ionic currents also depend on the gating variables, which are actually probabilistic values indicating the closing and opening of various ion channels. The time variation is given as where, and, Here, y can be any of the gating variables (m,h,….), τ y is the time constant, and y ∞ is the steady state value of y. α y and β y are voltage dependent rate constants.
The fast inward sodium current in eq. (3) was described using the Ebihara-Johnson model [4] (Ebihara and Johnson) coupled with the slow inactivation gate (j) of the Beeler-Reuter model (Beeler and Reuter) [5], where, Na G is the ,…maximum conductance of the sodium channel (23 mS/cm 2 ); m the activation parameter and h the inactivation parameter; and E Na is the reversal potential of sodium.  (10) and, d and f are two other gating parameters.
The delayed outward rectifier current I k was based on a formulation proposed by Shibasaki and included a time-dependent activation gate X, a time independent (fast) inactivation gate X i capturing the inward rectification properties of the delayed rectifier current, and a conductance modulated by extracellular K + concentration [6]. where, The formulation of I k1 was based on patch clamp data of Sakmann and Trube as well as other data from single ventricular myocytes [7].
( ) where, The voltage-dependent and time independent plateau current that activates at depolarized plateau potentials was identified by Yue and Marban [8], brings about a faster rate of repolarization at higher [K] 0 . Also the peak of I K1(T) occurs earlier (at more positive membrane potential), reflecting an increase of the reversal potential E K1 .
Furthermore, Figures 5-7 indicates that during repolarization the So, we can write an expression for extracellular potential, V e in terms of V m , as follows:

Results and Discussion
By developing a code using the FORTRAN software, we simulated the distribution of extracellular potassium by using the bidomain model with a LR action potential. Sepulveda et al. used a passive bidomain model to study unipolar stimulation of cardiac tissue [16]. They found that near a unipolar cathode, the tissue is strongly depolarized under the cathode, and more weakly hyperpolarized in two regions near the cathode along the fiber direction. Figure 1 shows our model generated transmembrane potential produced around the cathode, located at the origin.
In presence of the cathode the tissue is depolarized (positive transmembrane potential) in a dogbone shape. However, about 2 mm from the cathode, in the direction along the fibers, is a region of hyperpolarization (negative), called a "virtual anode". In Figure 1, only one quarter of the tissue is shown, so two virtual anodes exist one on each side of the cathode. The values in the plots are given in volts. Figure  2 shows anodic stimulation, which gives rise to two virtual cathodes on either side of the origin. Figure 3 shows the calculated transmembrane potential during unipolar cathodal stimulation of an anisotropic twodimensional sheet of passive cardiac tissue in a color plot.
A space clamped calculation of the LR action potential was performed in Luo and Rudy's original paper (Luo and Rudy) [17]. Figure 4 shows simulated action potentials for different extracellular concentrations of potassium.
As [K] 0 decreases, action potential duration increases and resting potential becomes more negative. Figure 5 shows the time course of I K1(T) , I K , and I si during the action potential. A peak in I K1(T) is observed during repolarization. For [K] 0 =7mM the peak is almost 3 µA/cm 2 but for [K] 0 =3mM it is less than 2 µA/cm 2 . This large difference in I K1(T) Figure 1: Transmembrane potential induced in a 2D sheet of tissue by unipolar extracellular cathode, located at the origin. The fibers are along the X axis.  The 1991 Luo-Rudy Phase I model did not describe intracellular Ca 2+ cycling, sarcolemmal Na + -Ca 2+ exchange, Na + -K + pumping and temporal variation of intracellular ion concentrations [17]. Luo and Rudy produced another model called Luo-Rudy phase II model, which is one of the most popular cardiac ventricular cell models even today (Luo and Rudy) [18]. However, because focus is made on extracellular potassium concentration and not on intracellular calcium dynamics, it is not relevant here.
Our next analysis is made by combining the bidomain model with an active Luo-Rudy membrane, which bears similarity with simulations performed previously by Roth [15]. Figure 8 shows stimulation threshold as a function of [K] 0. . The threshold initially decreases as [K] 0 increases, but then rises abruptly for [K] 0 above 12mM. Figure 9 shows the resting potential as a function of potassium ion concentration. The resting potential follows approximately the potassium Nernst potential, with a very negative value at low [K] 0 and then rising to less negative values as [K] 0 increases.

Vr (mV)
[K]0 (mM) Figure 9: Resting potential as a function of extracellular potassium ion concentration. The straight line represents the potassium Nernst potential. a case where propagation was successful in the normal tissue but then failed to enter the ischemic tissue. Also attempt was made by applying more than two stimuli (up to five) with various timings, and varied the distribution of extracellular potassium so that in one case the left was 5.4mM and the right was 8mM (Figure 13), and in another case the left was 8mM and the right was 12mM ( Figure 14). Nevertheless, after extensive simulations none of these simulations resulted in propagation on the left but then failure on the right. We validated our calculations by combining both models by verifying previous calculations, in that we reproduced results reported by Sepulveda et al. [16], Luo-Rudy action potential (Luo and Rudy, 1991) and by Roth and Chen [17,19]. This exercise brought about the conclusion that resting potential rises with [K] 0 , which tends to make propagation easier (lower threshold) initially as [K] 0 increases, but makes propagation harder (higher threshold) for very high [K] 0 (Figures 8 and 9).
Finally the Sidorov et al. experiment was simulated. In their case, the wave front propagated from a region of [K] 0 =4mM to a region of [K] 0 =10mM in a Langendorff-perfused rabbit heart [2]. They often found cases in which propagation was successful in the normal (4mM) tissue but then failed at the boundary between the normal and ischemic (10mM) tissue. Here this behavior was never observed. The reason appears to be the dependence of action potential duration on [K] 0. In the Luo-Rudy model, the action potential duration decreases as [K] 0 increases, so that the refractory period in the ischemic region is short. Any action potential that can propagate through the normal tissue (with a long refractory period) should have no problem also propagating through the ischemic tissue (with a short refractory period). Interestingly, in Sidorov et al. data, the action potential duration (and therefore the refractory period) appeared insensitive to extracellular potassium ion concentration. This appears to be the primary difference between our simulations (reported herein) and their experiments. Future simulations of Sidorov et al.'s data would need to include a model of the membrane kinetics that properly captures the [K] 0 dependence of the action potential duration.

Conclusions
We simulated Sidorov's experiment, in rabit heart, using Having verified that the model gives results similar to previous calculations, we can now apply the model to simulating the experiments of Sidorov et al. [2]. We can consider a tissue that has [K] 0 =4mM (normal) on the left half of the tissue, and [K] 0 =10mM (ischemic) on the right (Figure 10). The tissue sheet is 10 mm by 10 mm (100 × 100 nodes), with a space step in each direction of 0.1 mm.
We can stimulate the tissue by raising the extracellular potential along the left edge to 20 mV at time t=0. An action potential then propagates across the tissue from left to right. Plots of the transmembrane potential versus time at three locations are shown in

V (mV)
Time (milli-seconds) Figure 12: The transmembrane potential when the second stimulus is applied in the interval 542<t ≤ 543 ms. It can be observed that both, the first and second action potential waves, pass through the normal and ischemic tissue.
a mathematical model, wherein we incorporated Luo Rudy's representation of the cell membrane current with a bidomain model of the cardaic tissue.
In addition to matching experimental results by Sidorov, reliability of our model is also demonstrated by reproduction results of previously reported models that predicted dependence of extracellular potassium on action potential duration.
Using our model, we examined the propagation of potassium from low to high concentration and found that [K] 0 is dependent on the action potential duration. Further research is necessary to refine our modelling approach by making it more inclusive of physiological parameters.