Cardiac strength-interval curves calculated using a bidomain tissue with a parsimonious ionic current

The strength-interval curve plays a major role in understanding how cardiac tissue responds to an electrical stimulus. This complex behavior has been studied previously using the bidomain formulation incorporating the Beeler-Reuter and Luo-Rudy dynamic ionic current models. The complexity of these models renders the interpretation and extrapolation of simulation results problematic. Here we utilize a recently developed parsimonious ionic current model with only two currents—a sodium current that activates rapidly upon depolarization INa and a time-independent inwardly rectifying repolarization current IK—which reproduces many experimentally measured action potential waveforms. Bidomain tissue simulations with this ionic current model reproduce the distinctive dip in the anodal (but not cathodal) strength-interval curve. Studying model variants elucidates the necessary and sufficient physiological conditions to predict the polarity dependent dip: a voltage and time dependent INa, a nonlinear rectifying repolarization current, and bidomain tissue with unequal anisotropy ratios.


Introduction
An in-depth understanding of cardiac excitability is critical to the development of pacemakers and defibrillators. Modeling the dynamics of wave excitation in cardiac tissue is a way to predict the complex behavior of the tissue in response to an applied stimulus. Many models have been used to predict this behavior; typically these models have numerous parameters and ion currents, and each model has its own strengths and weaknesses [1]. The complexity of these models make it difficult to identify the relative role of parameters and variables on cell behavior such as the action potential duration (APD) and the excitation threshold. Other "phenomenological" models are simpler with fewer parameters, but are less realistic, and do not allow a link to detailed physiology such as voltage clamp experiments.
The strength-interval (SI) curve is a useful way to characterize the behavior of cardiac tissue in response to a stimulus [2]. Two stimuli are applied separated by an interval of time. The first (S 1 ) triggers an action potential, and the second (S 2 ) probes the tissue excitability as the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 tissue recovers from S 1 refractoriness. Previous experiments [3], as well as simulations using the bidomain model [4] with the Beeler-Reuter model [5], have shown the mechanisms of excitation underlying cathodal and anodal SI curves are cathode make, cathode break, anode make, and anode break. "Cathode" and "anode" refer to the polarity of the stimulus electrode, and "make" and "break" to the mechanism of tissue excitation; make is when excitation begins directly from depolarized tissue following the onset of the stimulus pulse, and break is when it begins (usually indirectly through hyperpolarized tissue) following the termination of the pulse.
Similar modeling studies [6], including one using the Luo-Rudy dynamic model [7,8], analyzed the "dip" in the anodal SI curve. Typically, the SI curve has negative slope, indicating that the S 2 threshold is lower at longer intervals when the tissue has recovered from S 1 refractoriness. For anodal stimulation, however, part of the SI curve has a positive slope, indicating that S 2 excitation is easier when the tissue is less recovered. This counterintuitive behavior is prominent in the break section of the anodal SI curve, and is suppressed or absent from the cathodal SI curve. One proposed mechanism responsible for the dip is that following the S 2 stimulus, the tissue surrounding the anode is hyperpolarized and excitable. The limiting factor determining the threshold is how well this excitable tissue interacts with nearby depolarization that can be produced by either the S 2 stimulus or the previous S 1 action potential. As the interval is shortened, more S 1 depolarization is present so less S 2 depolarization is required, lowering the threshold S 2 stimulus strength. Thus, the dip arises from the electrotonic interaction of excitable tissue with surrounding depolarization [2,6,7]. Other studies have claimed that the mechanism of the dip is caused by calcium dynamics [9] or is an intrinsic property of the membrane current [10,11].
A recently developed parsimonious ionic model [12,13] reproduces many experimentally measured characteristics of ventricular action potentials using only three variables, two ion currents, and eleven parameters. The advantage of a parsimonious model is that it highlights essential features of the membrane kinetics while not getting bogged down in nonessential details. The purpose of this manuscript is to combine the parsimonious ionic model with the bidomain formulation to study the events responsible for the shape of the strength interval curve. Specifically, we explore which features of the parsimonious model are essential for predicting the anomalous dip in the anodal SI curve. Our hypothesis is that a parsimonious model coupled with the bidomain formulation is sufficient to predict the shapes of the anodal and cathodal SI curves, and that the nonlinearity of the repolarization current is necessary for the dip in the anodal SI curve. unequal anisotropy ratios predict a complicated transmembrane potential distribution with adjacent regions of depolarization and hyperpolarization [3,16] for both anodal and cathodal stimulation.
Cathodal SI curve. The SI curve for cathodal stimulation (Fig 1, blue curve) shows a decrease in threshold S 2 stimulus strength as the tissue is recovering from the S 1 action potential. At long intervals (rheobase), the cathodal threshold is 30 mA/cm. At an interval of 155 ms the tissue is nearly recovered from S 1 refractoriness and cathode make excitation originates from the depolarization surrounding the cathode. Details of this behavior are shown in Fig 2 and S1 Movie (interval = 155 ms, S 2 stimulus = 33 mA/cm; see supplementary information). At shorter intervals cathode break excitation occurs. The tissue is still refractory at 143 ms when the S 2 pulse begins. During and immediately after the S 2 stimulus the depolarization is extensive and hyperpolarization along the fiber direction is evident (the "virtual anode"). Break excitation occurs following the end of the S 2 pulse, when the depolarization surrounding the cathode diffuses into the hyperpolarized virtual anode, exciting a wave front that initially propagates through the hyperpolarized and therefore excitable tissue along the fiber axis (Fig 3  and S2 Movie, 143 ms, 435 mA/cm). At less than 140 ms, the tissue is unexcitable for any S 2 strength (absolute refractory). Anodal SI curve. The anodal SI curve (Fig 1, red curve) has a prominent dip (a local minimum created by an abrupt drop at 127 ms followed by a region of positive slope), as seen in other models [2,6,7]. The threshold value of the anodal stimulus at the minimum of the dip is 587 mA/cm. For S 1 -S 2 intervals below 127 ms the tissue is unexcitable. The maximum The transmembrane potential is shown as a function of x (horizontal, parallel to the fibers) and y (vertical, perpendicular to the fibers). The black rectangle is the electrode and the number in each frame is the time in ms. S 1 -S 2 interval = 155 ms, S 2 stimulus = 33 mA/cm. These frames are taken from the S1 Movie (see supplementary files). The color scale for this and other plots of transmembrane potential is shown in the movie; purple is rest, yellow is depolarized, and blue is hyperpolarized.  threshold value is 1017 mA/cm recorded at 147 ms. The range from 155 to 160 ms contains a supernormal period similar to that in the space-clamped simulations, during which the threshold dips to at most 3% below its asymptotic long-interval (rheobase) value of 598 mA/cm. Anode make excitation occurs at long intervals. The hyperpolarization surrounding the anode plays little role in excitation, and the wave front originates from the depolarized region (the "virtual cathode") along the fiber direction (Fig 4 and S3 Movie, 158 ms, 640 mA/cm). The supernormal period (155-160 ms) was restricted to the make section of the SI curve. Anode break excitation occurs at shorter intervals. The depolarization along the fiber axis diffuses into the hyperpolarized region surrounding the anode, exciting a wave front that propagates initially perpendicular to the fiber direction ( Fig 5 and S4 Movie, 130 ms, 710 mA/cm). At intervals shorter than 127 ms, anode break initiates an action potential but it can only propagate locally near the electrode because the surrounding tissue is too refractory and therefore is incapable of generating a wave that propagates through the heart (Fig 6 and S5 Movie, 125 ms, 600 mA/cm).

Equal anisotropy ratios
The SI curves were recalculated using equal anisotropy ratios, so no virtual anode was produced during cathodal stimulation and no virtual cathode occurred during anodal stimulation. Anodal stimulation did not excite the tissue. Cathodal stimulation was always by the make mechanism and not the break mechanism.

Effect of I K current on strength-interval curves
The repolarization current I K contains three parameters (see Methods section), each of which modifies the I K current-voltage curve in a unique way. The conductance g K scales the magnitude of the I K current, E K is the reversal potential, and k r controls the nonlinearity of the current-voltage curve.
The variation of I K with transmembrane potential for different k r values is shown in Fig 7a. The values for k r are chosen to show the effect of the I K current when it is linear (k r = 1) and nonlinear (k r = 21.28 and 40 mV). The value k r = 21.28 mV corresponds to that used in Gray's original model [12,13]. When g K and E K are fixed (0.3 mS cm -2 and -83 mV), the current-voltage curves have the same intercept and slope for I K = 0.   almost constant with an increase in S 1 -S 2 interval and the rheobase strength is highest of all k r values examined. The APD and the amplitude of the action potential depend upon k r . With a nonlinear I K current the anodal SI curve exhibits a transition from break to make excitation; for k r = 21.28 mV the transition occurs at 155 ms, and for k r = 40 mV at 38 ms. The longer APD (k r = 21.28 mV) results in a SI curve with a smaller threshold dip than the short APD (k r = 40 mV). When k r = 40 mV, the anodal SI curve has a dip at 19 ms with a stimulus strength of 981 mA/cm. The maximum threshold current is 1565 mA/cm at 31 ms. For k r = 21.28 mV, anodal SI curve dip is at 127 ms with a strength of 587 mA/cm. Rheobase is highest for k r = 1.

Effect of k r and g K parameters on strength-interval curves
One limitation of varying k r by itself is that the action potential duration changes dramatically. In order to look at changes in the nonlinearity of I K without confounding changes in APD, we varied both k r and g K simultaneously to keep the APD approximately constant . Fig 7b shows that k r = 40 mV and g K = 0.06 mS cm -2 produces a triangular shaped action potential, and that k r = 1 mV and g K = 0.016 mS cm -2 produces a decaying exponential action potential, both with a duration similar to the standard action potential (k r = 21.28 mV and g K = 0.3 mS cm -2 ). Fig 7a shows the variation of the I K current with membrane potential. The SI curves of these three cases are shown in Fig 9. The anodal strength-interval curve for the linear case (k r = 1) does not contain a section for anode break excitation, and does not contain a dip.
A nonlinear I K current is necessary for anode break excitation. To examine why, we performed simulations of anodal excitation for an interval of 130 ms, for k r = 21.28 mV and g K = 0.3 mS cm -2 (normal, nonlinear) and for k r = 1 and g K = 0.016 mS cm -2 (linear).  (Fig 12a) and at the virtual cathode (Fig 12b). In the normal case, the hyperpolarization surrounding the anode is relatively weak and decays quickly after the stimulus pulse ends, creating an excitable path for the break wave front [5,6,17]. In the linear case, the hyperpolarization surrounding the anode is strong and decays slowly. The depolarization does not diffuse into hyperpolarized and excitable tissue. Rather, hyperpolarization diffuses into depolarized and refractory tissue, which does not result in excitation [18]. The nonlinearity of the repolarization current is critical to reducing the extent of the hyperpolarization.

Effect of E K on strength-interval curve
At longer time intervals the primary mechanism of anodal excitation is make, whereas for shorter time intervals the primary mechanism is break [19]. These effects are evident in our calculations with the parsimonious ionic model (anode make in Fig 4 and anode break in Fig  5). When we increase the reversal potential (E K ) to -69 mV, which represents an elevation in extracellular potassium ion concentration, at long time intervals the excitation mechanism is a combination of make and break (Fig 13 and S8 Movie, 200 ms, 500 mA/cm). Excitation with a weaker stimulus just above threshold (200 ms, 458 mA/cm) is purely break excitation, so we will categorize this more complex behavior as break. The change from make to break excitation at long intervals is consistent with the experimental work on elevated extracellular potassium ion concentration by Sidorov et al. [20], who showed that in diastole the anodal excitation changes from make to break. Numerical simulations carried out by Roth and Patel [21] using the Luo-Rudy dynamic model [8] suggest that potassium ion concentration influences anodal stimulation of cardiac tissue. Fig 14 shows the anodal and cathodal SI curves for E K = -83 and -69 mV. The action potential duration shortens with an increase in E K (Fig 15), causing the SI curves to shift to the left. Although the cathodal SI curve retains its shape the anodal SI curve for elevated potassium does not contain a dip. Therefore, elevated E K suppresses the dip in the anodal SI curve, consistent with previous studies [20,21]. It also suppresses the supernormal period.

Discussion
Our simulations using the parsimonious ionic model for SI curves are consistent with results from previous studies [2,6,7,20,21]. The presence of the dip using a parsimonious model supports our hypothesis that the shape of the strength-interval curve is strongly influenced by the Cardiac strength-interval curves calculated using a bidomain tissue with a parsimonious ionic current voltage dependent membrane recovery and bidomain properties of the tissue. By changing the parameter k r in the parsimonious ionic model we were able to obtain SI curves with shorter APD yet which have a transition of anode break excitation to anode make excitation (for k r = 21.28, 40 mV) as shown in movies in the supplementary files. The presence of the dip only at non-infinite values of k r suggests that a nonlinear rectifying repolarization current is necessary for anode break excitation and the dip to occur in the anodal SI curve. Moreover, the parsimonious ionic model predicts the effect of elevated extracellular potassium concentration. Roth and Patel [21] found that the disappearance of the dip was caused by post-repolarization refractoriness, so that when the tissue recovered from refractoriness there was no S 1 depolarization to assist with anode break excitation.
Why did we choose to use a parsimonious model rather than a more detailed model for the ionic current? Gray and Pathmanathan [13] make a compelling case for the utility of minimal models that contain few parameters. Most more detailed ionic models are derived from voltage clamp data obtained under non-physiological conditions. Complex models suffer from overparameterization, non-uniqueness, and have limited validity beyond the data used in their derivation. The parsimonious model is mathematically identifiable [12] and exhibits important emergent phenomena such as alternans and spiral wave breakup [13]. This model is attractive for us because we wish to know if the shape of the strength-interval curve can be predicted using even a simple model of cardiac electrical behavior. Our conclusion is that the parsimonious model can predict the shape of the SI curve, and that more complicated and exotic ionic currents are not necessary to explain this behavior. The parsimonious ionic model has its limitations. APD depends upon the degree of linearity of I K . The ideal model is such that we should be able to vary the degree of linearity without a change in APD, so that we can isolate the effect of nonlinearity. There is a discrepancy between the quantitative value of the threshold currents for SI curves compared to previous studies. The currents are relatively large for the parsimonious ionic model. There can be numerous reasons for this discrepancy [22].
In conclusion, a simple ion current model that includes only an excitable sodium current and a nonlinear repolarization current is sufficient to predict the dip in the anodal SI curve. The parsimonious ionic model with the bidomain tissue model can predict shapes of the anodal and cathodal SI curves and break and make excitations for both anodal and cathodal stimuli. Nonlinearity of the repolarization current and electrotonic interaction of depolarization and hyperpolarization are both necessary for a dip in the anodal strength-interval curve. Other researchers have postulated that the dip may be caused by the "funny current" that activates upon hyperpolarization but that has a depolarized reversal potential so it supplies inward current near rest potential [10,11], or by calcium dynamics [9]. Although we cannot rule out a role for such currents with our analysis, this parsimonious model indicates that such factors are not necessary to reproduce the observed strength-interval curves.

Methods
The two dimensional bidomain model [23] is used to represent the cardiac tissue. The active membrane current is defined by the parsimonious ionic model [12,13], which consists of two currents: a sodium current I Na that activates rapidly upon depolarization, and a time-independent inward-rectifying repolarization current I K I Na ¼ g max Na m 3 hðV À E Na Þ ð1Þ Cardiac strength-interval curves calculated using a bidomain tissue with a parsimonious ionic current where with m 1 ðVÞ ¼ 1 The bidomain model is, ðg ix þ g ex Þ @ 2 V e @x 2 þ ðg iy þ g ey Þ @ 2 V e @y 2 ¼ À g ix where V m and V e are the transmembrane and extracellular potentials, C m is the membrane capacitance per unit area, and β is the surface-to-volume ratio. Values for the intracellular and extracellular conductivities in the directions parallel and perpendicular to the fibers, g ix , g iy , g ex , g ey (Table 1), are chosen such that the tissue has unequal anisotropy ratios [16]. The stimulus current is applied at the center of the tissue using a rectangular unipolar electrode having a length (x-direction) and width (y-direction) of 0.1 cm. The stimulus is applied to a tissue of 1.0 cm by 1.0 cm by applying boundary conditions to the electrode tissue surface: we consider that the total current applied to the electrode is transferred to the extracellular space as the stimulating current, the normal component of the intracellular current density is zero, and the extracellular potential is kept constant over the electrode surface. All stimuli have a duration of 2 ms. Before stimulation the resting potential is equal to E K . In each simulation, the initial stimulus (S 1 ) is cathodal with twice the rheobase stimulus strength. Rheobase depends on the model parameters, so the S 1 strength in mA/cm varies between simulations. All movies were generated using 1 ms per frame. Only one quadrant of the tissue is shown. The transmembrane potential in the other quadrants can be found from the even symmetry. Movies were generated using an S 2 strength about 10% above the threshold strength for that interval.
Space-clamped simulations were performed using Eq 9 with the left-hand-side set to zero, and the stimulus was applied through an additional membrane current added to I K and I Na .
When equal anisotropy ratios were used, the conductivities were chosen to be, g ix = g ex = 2 mS cm -1 , and g iy = g ey = 0.32 mS cm -1 .
The bidomain equations are solved numerically with the finite difference method. Using an initial value of V e (t), we solve Eq 9 for V m (t+Δt). Then we solve Eq 10 for V e (t+Δt) using the value of V m (t+Δt) in the source term, and the successive overrelaxation (SOR) method [4,6]. In order to accelerate convergence, an overrelaxation parameter of w = 1.8 was used [24]. The iterative loop terminates when changes in V e between subsequent time steps is less than 1 μV. All programs are written in Fortran 90 and compiled using the gfortran compiler (www.gcc. gnu.org). The space step is 0.005 cm parallel to the fiber direction and 0.002 cm perpendicular to the fiber direction. The time step is 1 μs. The number of nodes in x-direction is 200 and the number of nodes in y-direction is 500. Because the tissue is stimulated from the center and propagation is symmetrical, we consider only one quadrant of the tissue.
During simulations of anodal excitation, we encountered a computational problem arising from the strong hyperpolarization next to the electrode. The time constant τ h became too small (even smaller than the time step), making the calculations unstable (Fig 16). In order to correct this, we modified the parsimonious ionic model by requiring that the minimum value of τ h be the time step (1 μs). This modification made the simulations stable.