Testing the limits of cardiac electrophysiology models through systematic variation of current

1.   Introduction

Mathematical models of cardiac electrophysiology are being developed in the hope that they might one day be used to guide personalized treatment of cardiac diseases ^[29]. Unfortunately, detailed models can be problematic ^[35]. The large number of free parameters in complex models can make parameter identifiability difficult ^[10], particularly in a time frame that is clinically useful. It is also intuitively difficult to understand how changes in one parameter alter model predictions and whether estimated parameter values fall within a biologically realistic range. Finally, complex models have problems with numerical stability ^[9,31] and uniqueness of the solutions to the model ^[47,67]. This limits the models' usefulness in making predictions about scenarios outside the "normal" range, such as patients with genetic mutations changing their electrophysiology ^[94] or patients being treated with antiarrhythmics ^[33].

Mathematical models of cardiac cell ionic currents developed as extensions of the Hodgkin-Huxley model of the squid giant axon ^[34], with the first simply including a second potassium current ^[55] along with Hodgkin and Huxley's sodium, potassium, and leakage currents. Since then, increasingly complicated ionic current models have been developed for different animal species and different types of cardiac cells within the heart ^[30,56]. New currents are added based on experimental studies that measure individual ionic currents ^[13,19] and the resulting new models are often validated by comparing to experimental data consisting of multi-current processes such as an action potential ^[12,70]. Beyond this broad validation of a model, there is not much study of how the addition of a new current alters the dynamics of existing model currents and there is usually no experimental data on how changes in one current alter the dynamics of the other currents or the entire system ^[20,80].

This limited validation and testing of new models has consequences for the attempt to personalize models for clinical use. When personalizing a model, parameters of the model are chosen to fit some aspect of a patient's electrophysiological response. These parameters are used in the model to simulate a patient's response to treatment, meaning that some of these parameters are further perturbed from their base values. Since most models are so finely tuned to the experimental conditions from which they derive, it is not clear that perturbing the model parameters in this way will lead to accurate predictions. Researchers have started to examine this by studying populations of models with parameters randomly drawn from specified distributions ^{[46,52,54,78,86]}. An alternative approach is a sensitivity analysis to examine how different parameters affect a certain outcome ^[79]. Both of these techniques have been useful for correlating subsets of parameter values with certain desired (or undesired) outcomes.

In this paper, we use a more systematic approach to examining the limitations of a cardiac model by continuously varying the value of a single parameter. We focus on the effect of currents, specifically the sodium, potassium, and calcium currents and examine how variation in the conductivities alters the predicted model outcome. Changes in current conductivities can often be reproduced experimentally through the use of antiarrhythmic drugs ^[75], so model predictions of the effect of changes in current can be verified experimentally. By systematically varying parameters, we can identify specific values at which there might be problems with stability or with realistic representation of the electrophysiological response. Once anomalous points have been identified, comparison with experiment can help pinpoint which currents might be inaccurately or incompletely formulated, leading to further improvement of the model and an extension of its range of validity.

2.   Materials and method

We use the Bernus human ventricular epicardial cell model ^[7] and the Fox canine ventricular cell model ^[24] to compute the total ionic current through the cell membrane at a given time. These models are not the most complex cardiac models, but do have several currents that have been added to or modified from earlier models without examining the interactions with existing currents. All currents in both models are Hodgkin-Huxley type currents controlled by gating variables. The simplicity of the models allows for faster numerical computation and avoids the computational instability seen in more complex models ^[88], particularly when model parameters are varied ^[47]. We also use the Fox model since this model exhibits steady-state alternans ^[24], a precursor to more complex arrhythmias ^{[53,69,76,81]}. Using two models that exhibit different dynamics allows us to assess whether our results are model-dependent.

The Bernus model has 9 different ionic currents that include a rapid sodium current $I_{Na}$, a slow calcium current $I_{Ca}$, a transient outward current $I_{to}$, a delayed rectifier potassium current $I_{K}$, an inward rectifier potassium current $I_{K1}$, background sodium and calcium currents $I_{Na, b}$, $I_{Ca, b}$, a sodium-potassium pump $I_{NaK}$, and a sodium-calcium exchanger $I_{NaCa}$. The Fox model includes 13 ionic currents. In addition to the currents present in the Bernus model, it has a plateau potassium current $I_{Kp}$, a potassium current through calcium channel $I_{CaK}$ and a sarcolemmal pump current $I_{pCa}$. Also, the delayed rectifier potassium current $I_{K}$ includes a rapid component $I_{Kr}$ and a slow component $I_{Ks}$.

The cell membrane potential $V$ in both models follows an ordinary differential equation,

where $I_\mathrm{stim}$ is an external stimulus current, $C_m$ is the cell membrane capacitance per unit surface area and $I_\mathrm{ion}$ is the sum of all the membrane currents.

Our implementation of both models uses the $\mathtt{Odeint}$ function from SciPy (open source Python library) to solve the system of ordinary differential equation. $\mathtt{Odeint}$ dynamically switches from the Adams method to the backward differentiation (BDF) method as needed with a tolerance of 1.5$\times 10^{-8}$. We use the same values for the various parameters as in the original derivation of the models ^[7,24], except as noted below.

2.1.   Pacing protocols

We study the effect of the sodium, potassium or calcium ionic currents by multiplying the corresponding conductance by $(1-\varepsilon)$, where $\varepsilon$, which ranges from 0 to 1, is the relative decrease of the current. Specifically, we apply $\varepsilon$ to $I_{Na}$ and $I_{Ca}$ in both models. While for potassium currents, we apply $\varepsilon$ to $I_K$ in the Bernus model and $I_{Ks}$ in the Fox model. While this type of systematic variation can be used on any parameter, we choose to vary conductances since there is some limited experimental data where specific currents have been reduced or blocked ^[8,26,27,38].

For the Bernus model, action potentials are produced with an external stimulus current of 19.6 pA/pF of 2 ms duration. For the Fox model, action potentials are produced with an external stimulus current of 80.0 pA/pF of 1 ms duration. Action potential durations (APDs) for both models are measured using a threshold membrane potential corresponding to 90% repolarization.

Since we are interested in studying how different drugs alter cardiac cell dynamics, we generate bifurcation diagrams for a variety of changes in conductance. Bifurcation diagrams show the relationship between APD and basic cycle length (BCL) ^[32], which is the time interval between consecutive stimuli. The diagram is obtained by recording APD for different values of BCL, first with BCL decreasing (BCL downsweep) to a certain fixed value and then increasing the BCL (BCL upsweep) from its lowest value to its original value. For the Bernus model, we begin at a BCL of 400 ms and step down to a BCL of 200 ms with a step size of 1 ms, then return to our original BCL using the same step size. We recorded 20 APD for each BCL. Note that the Bernus model reaches steady-state within two or three stimuli. For the Fox model, we begin at a BCL of 300 ms and step down to a BCL of 90 ms with a step size of 1 ms, then return to our original BCL using the same step size. For the Fox model, we pace for 40 stimuli and then record the next 20 APD as the Fox model requires at least 30 stimuli to attain the steady state.

Finally, we examine the effect of reducing currents on restitution properties of the cells. Restitution curves are relationships between the APD and diastolic interval (DI), which is the time between the end of one action potential and the beginning of the next. Several different restitution curves have been defined ^[42] based on different pacing protocols, and the resulting curves are thought to be related to the stability of the 1:1 cardiac response ^[43]. Here we examine the two most commonly measured restitution curves, the dynamic restitution curve (DRC) and the S1S2 restitution curve (SRC). The DRC consists of the steady-state values of APD and previous DI, while the SRC consists of the APD response to a perturbation in DI from steady-state values. We simulate the DRC by pacing for 20 stimuli in the Bernus model and for 40 stimuli in the Fox model to reach steady state. The last APD and preceding DI are recorded and used as one point of the DRC. The BCL is decreased and the procedure is repeated until we no longer see a 1:1 response. The SRC is generated by pacing for 20 stimuli at a BCL of 600 ms in the Bernus model and 40 stimuli at a BCL of 400 ms in the Fox model. A single stimulus of a different BCL is given and the resulting DI and APD are recorded as part of the SRC. The procedure is repeated using the same base BCL, but different BCLs for the perturbation.

3.   Results

3.1.   Changes in the action potential

A basic element of cardiac electrophysiology is the action potential, so we first use our methodology to contrast and compare the effects of different currents on the action potential. Initial conditions from the original derivations of the models are used to initiate the simulation which is allowed to run for 100 ms to ensure that the model is at the resting membrane potential. A single stimulus pulse is then applied to either model, and we look at changes to the subsequent action potential as $\varepsilon$ varies. The results are shown in Figure 1. Both models show similar responses for all three currents. When the sodium current decreases (Figure 1, A & D), the models predict that the APD does not change for most $\varepsilon$. For the Bernus model, the action potential loses its spike and the AP disappears at high $\varepsilon$, while for the Fox model, the action potential is substantially reduced at high $\varepsilon$, but does not completely disappear. For a reduction in potassium current (Figure 1, B & E), the models predict that the APD increases with increase in $\varepsilon$. The Fox model suggests that changes in the potassium current also alter the resting potential, although we do not see a similar effect in the Bernus model. In the case of calcium current (Figure 1, C & F), the models predict that the action potential duration decreases with increase in efficacy. The effect of reduction in calcium current is to shorten the duration of the plateau phase of the action potential. Even when the calcium channel is completely blocked, there is still an action potential, although it does not have a plateau.

Figure 1.  Effect of decreasing Na, K and Ca currents on action potentials of the Bernus (top row) and the Fox model (center row). (A & D) For Na current, APD remains almost constant for a wide range of $\varepsilon$, markedly decreasing at high values. (B & E) For K current, APD increases with $\varepsilon$, preventing the transmembrane voltage from returning to the rest state when $\varepsilon$ exceeds $\sim$0.6. (C & F) Decrease in the Ca current decreases the APD by decreasing the plateau phase. Arrows indicate the trend in the direction of decreasing current. The three currents have very different effects on APD, as shown in the summary diagrams for (G) the Bernus model and (H) the Fox model, with reductions in Na current (black line) showing little effect, reductions in K current (red line) causing increase in APD, and reductions in Ca current (blue line) causing almost linear decrease in APD.

DownLoad: Full-Size Img PowerPoint

Each of the three currents under study has a distinct effect on APD. These findings are summarized in Figure 1 G & H, where we have plotted the action potential duration as a function of $\varepsilon$. These summarizing figures allow us to clearly contrast the effects of the currents over their full range of effect. Reduction in sodium current has no effect on APD until $\varepsilon\sim0.94$ for both models. For slightly higher values of $\varepsilon$, the APD increases slightly before the action potential is completely blocked in the Bernus model. However, in the Fox model, the APD still exists as the action potential is not completely blocked. Reduction in the potassium current causes a rapid increase in the APD as $\varepsilon$ increases, while reduction in the calcium current causes a slower decrease in the APD as $\varepsilon$ increases in both models.

3.2.   Rate-dependent changes

Action potentials change as the pacing rate changes and interactions between currents can change as gates do not have time to open or close completely. We examined the rate dependence of APD as various currents were changed by pacing at a variety of BCLs until a steady-state response was reached. We constrained our study to 1:1 responses and for this range of BCLs, the Bernus model reached a steady-state APD within two or three action potentials, as can be seen in the bifurcation diagrams presented later, whereas the Fox model reached a steady state after about 30 action potentials ^[22]. The last action potential was used to calculate the APD for each BCL. The lengthening or shortening of APD, APD_Diff is calculated by taking the difference between APD without changes in the currents and APD with one of the currents reduced.

Figure 2 shows the lengthening or shortening of APD as a function of the BCL for several values of $\varepsilon$ for the Bernus and the Fox models, respectively. For sodium current (Figure 2 A & D), our models predict that there is no appreciable change in APD at different BCLs in the Bernus model, no matter how much sodium current is reduced. This is consistent with the single AP results which indicate that changes in sodium current cause little change in the APD (Figure 1 A & D). However, we notice a slight increase of APD_Diff with increase in $\varepsilon$ in the Fox model at lower BCLs. Potassium current increases the APD (Figure 1 B & E), so we expect a positive APD_Diff as seen in (Figure 2 B & E). Here, both models predict that there is a larger effect on APD when potassium current is reduced as BCL increases. APD_Diff is somewhat higher in the Fox model as compared to the Bernus model, suggesting that the changes in potassium current have a larger effect in the Fox model. Calcium current decreases the APD (Figure 1 C & F), so we have a negative APD_Diff (Figure 2 C & F). Again, the models predict an increased effect at longer BCLs in both models.

Figure 2.  Changes in the APD, APD_Diff induced by reduction in Na, K and Ca currents for the Bernus (top row) and Fox model (center row). (A & D) There is little change in APD_Diff by reduction in Na current at any pacing rate. (B & E) APD is lengthened by reduction in potassium current with greater changes at longer BCLs. (C & F) APD is shortened by reduction in calcium current with greater changes at longer BCLs. Arrows indicate the trend in the direction of decreasing current. The rate of change of the drug-induced APD_Diff with BCL is plotted as a function of $\varepsilon$ for the three currents for (G) the Bernus model and (H) the Fox model. There is little change in the slope for changes in Na current (black line). For K current (red line), the slope grows with increase in the $\varepsilon$. For Ca current (blue line), the slope decreases linearly at slower rate and then decreases at higher rates with increase in $\varepsilon$.

DownLoad: Full-Size Img PowerPoint

We summarize the changes in rate-dependent behavior by calculating the slope of APD_Diff vs. BCL curves (Figure 2 G & H). While the curves are not all linear, as a simple approximation of the rate-dependent changes, we fit a line to the curve to calculate the slope. Interestingly, the general pattern here is the same as for changes in APD — blocking sodium channels has almost no effect on the slope; blocking potassium channels causes a rapid increase in the slope; and blocking calcium channels causes a more gradual decrease in the slope. Thus blocking potassium channels results in larger rate-dependent effects than changes in other currents.

3.3.   Changes in the bifurcation diagram

Bifurcation diagrams provide a detailed picture of the electrophysiological response of cardiac cells by recording APD at various BCLs. The bifurcation diagram for the Bernus model in the absence of any changes is shown in Figure 3 (left). During the downsweep (right to left) the APD decreases but remains in a 1:1 pattern until a BCL of 275 ms when the cell reaches a point at which it cannot recover fast enough to respond to every stimulus. There is an abrupt increase in the APD as the cell transitions to the 2:1 response and remains in the same pattern until the BCL decreases to 200 ms. As we begin the upsweep (left to right), the cell remains on the 2:1 branch until BCL = 333 ms after which the cell makes an abrupt transition back to 1:1 behavior. The bifurcation diagram for the Fox model is shown in Figure 3 (right). As BCL is decreased it remains in a 1:1 pattern until the BCL reaches 200 ms. At this BCL, alternans occur showing the bifurcation from a 1:1 to a 2:2 response. The cell continues with 2:2 behavior until the BCL reaches 150 ms. The cell returns to 1:1 again and maintains this pattern until the BCL reaches 95 ms when there is an abrupt increase in the APD as the cell transitions to the 2:1 response and remains in the same pattern until the BCL decreases to 90 ms. In the Fox model, there is little difference in behavior when increasing or decreasing BCL, so we do not observe a hysteresis window in this model.

Figure 3.  Bifurcation diagram for the Bernus model (left) and the Fox model (right) without changes in currents. The red triangles correspond to APD during BCL downsweep (right to left) and the black diamonds correspond to APD during BCL upsweep (left to right).

DownLoad: Full-Size Img PowerPoint

Bifurcation diagrams were generated and examined for a variety of values of $\varepsilon$ for each of the currents for both models. We adjust the different ranges of BCLs for some values of $\varepsilon$ in order to see the transitions as APD gets lengthened or shortened with changes in $\varepsilon$. The results for the Bernus model are shown in Figure 4. No matter which current is altered, we find that there is a transition from a 1:1 to a 2:1 pattern and that there is a hysteresis window associated with this transition. For sodium current, we can see hysteresis up to an efficacy of 0.95, after which we do not see hysteresis since the action potential vanishes. For potassium current, $\varepsilon$ is limited to 0.6 as higher values lead to very large APDs exceeding the BCL. When potassium is altered, we see an alternating rhythm appear between BCL≈300–400 ms and move towards longer BCL≈525 ms as $\varepsilon$ increases. Note that the usual definition of alternans is two different responses for every two stimuli, but here we observe two different responses for every four stimuli. For calcium current, the transitions shift to lower BCLs due to shortening of the APD.

Figure 4.  The effect of reductions in current on rate-dependent behavior in the Bernus model. We show the bifurcation diagrams in the presence of various values of $\varepsilon$ for Na current (top row), K current (center row) and Ca current (bottom row).

DownLoad: Full-Size Img PowerPoint

Bifurcation diagrams for the Fox model are shown in Figure 5. Compared to the Bernus model, we can see a completely different picture in this model. Reduction in sodium current shifts the APD alternans bubble at lower and intermediate values of $\varepsilon$ and causes the alternans to vanish when applied at high values of $\varepsilon$. With reduction in potassium current, APD alternans exist only at low values of $\varepsilon$, while higher values of $\varepsilon$ cause response patterns not seen in the untreated model. Bifurcation diagrams for reduced calcium current are quite different compared to the two other currents; even at low values of $\varepsilon$, reducing calcium current causes the alternans and other transitions to disappear, leaving a 1:1 response at all BCLs.

Figure 5.  The effect of current reduction on rate-dependent behavior in the Fox model. We show the bifurcation diagrams in the presence of various values of $\varepsilon$ for Na current (top row), K current (center row) and Ca current (bottom row).

DownLoad: Full-Size Img PowerPoint

The hysteresis window has been linked to arrhythmogenesis in both modeling ^[36,59,91] and experimental ^[58] studies. We can study how the hysteresis window changes as we reduce different currents by defining ${\mathrm{BCL}_\mathrm{down}}$ as the BCL at which we see a transition from 1:1 to 2:1 and ${\mathrm{BCL}_\mathrm{up}}$ as the BCL at which we see a transition from 2:1 to 1:1. These transition BCLs are plotted as a function of $\varepsilon$ and the results are shown in Figure 6. The change in the width of the hysteresis window will only be studied for the Bernus model as there is no significant hysteresis window in the Fox model. For sodium current, the transitions, both up and down, take place at almost the same BCLs for a wide range of $\varepsilon$, only increasing slightly once $\varepsilon$ surpasses ~0.9. For potassium current, the transitions, both up and down, occur at longer BCLs as $\varepsilon$ increases. This is a direct consequence of the increase in APD as $\varepsilon$ increases. ${\mathrm{BCL}_\mathrm{up}}$ increases smoothly, but ${\mathrm{BCL}_\mathrm{down}}$ exhibits a sharp discontinuity near $\varepsilon\approx 0.4$, leading to a decrease in the size of the hysteresis window at higher values of $\varepsilon$. For calcium current, the transitions both occur at shorter BCLs as $\varepsilon$ increases. ${\mathrm{BCL}_\mathrm{up}}$ appears to decrease more rapidly than ${\mathrm{BCL}_\mathrm{down}}$ causing a decrease in the size of the hysteresis window as $\varepsilon$ increases with a complete disappearance of the window once $\varepsilon > 0.88$.

Figure 6.  The effect of current reduction on the hysteresis window. For each of the currents, we plot the BCL at which the dynamics transition from 1:1 to 2:1 behavior (${\mathrm{BCL}_\mathrm{down}}$) and the BCL at which the dynamics return from 2:1 to 1:1 behavior (${\mathrm{BCL}_\mathrm{up}}$). (left) For Na current, the size of the hysteresis window (difference between ${\mathrm{BCL}_\mathrm{up}}$ and ${\mathrm{BCL}_\mathrm{down}}$) remains constant for almost the entire range of $\varepsilon$. (center) For K current, the transitions, both up and down, occur at longer BCLs as $\varepsilon$ increases. (right) For Ca current, the transitions both occur at shorter BCLs as $\varepsilon$ increases.

DownLoad: Full-Size Img PowerPoint

3.4.   Restitution

Restitution curves are theoretically linked to the stability of the 1:1 response ^[85] and provide a measure of the cell's response to changes in pacing rate. We simulate the DRC and SRC in both models at different values of $\varepsilon$ and determine the maximum slope. Figure 7 shows the maximum slope of the DRC (top row) and the SRC (bottom row) for the Bernus model (left column) and the Fox model (right column).

Figure 7.  The effect of current reduction on restitution curves. For each of the currents, we plot the maximum slope of either the DRC (top row) or the SRC (bottom row) for the Bernus (left column) or Fox (right column) models.

DownLoad: Full-Size Img PowerPoint

Reducing the sodium current does little to change the maximum slope of either the DRC or the SRC in the Bernus model until a large reduction in current. Near a reduction of 95%, the RCs get steeper before the blockage of sodium current results in a loss of AP. Blocking the potassium current results in steepening of both the DRC and SRC in the Bernus model. The calcium current shows slightly different responses in the DRC and SRC when it is blocked in the Bernus model. The maximum slope of the DRC remains the same until a roughly 40% block of calcium current. The maximum slope of the DRC then decreases as the calcium current is further reduced. The maximum slope of the SRC initially increases as calcium current is blocked, but then decreases as calcium current is increased beyond ~45%.

The slopes of the restitution curves in the Fox model respond differently than in the Bernus model to changes in sodium, potassium, and calcium currents. The maximum slope of the DRC gradually decreases as the sodium current is reduced until the abrupt increase at a reduction of ~95%. The maximum slope of the SRC, on the other hand, remains almost constant as sodium current is initially reduced (up to ~40%). It then gradually rises to a peak value near $\varepsilon\sim 0.6$ followed by a gradual decline until $\varepsilon\sim 0.95$ when it declines more rapidly. A reduction in potassium current in the Fox model results in an initial decline in the maximum slope of the DRC (until $\varepsilon\sim 0.1$), followed by an increase (until $\varepsilon\sim 0.4$) and a rapid decline thereafter. The maximum slope of the SRC initially increases (until $\varepsilon\sim 0.2$) followed by a more rapid reduction in the slope of the SRC. Reduction of the calcium current initially causes a slight increase in the maximum slope of the DRC, followed by a gradual decrease. The maximum slope of the SRC does not show a similar initial increase as calcium current is reduced, but simply decreases as calcium current is reduced.

4.   Discussion

In the present study, we used mathematical models of cardiac cells to systematically study the effect of changing ionic conductances. Our study revealed several features that are not likely to be found by examining populations of models or doing a sensitivity analysis. These features include the slight increase in APD just before total block of the sodium current in the Bernus model, the change in slope of rate-dependence as calcium conductance is decreased, the discontinuity in ${\mathrm{BCL}_\mathrm{up}}$ as potassium conductance is decreased in the Bernus model, or the differences between maximum slopes of DRC and SRC in the Fox model. Whether these features are characteristics of these particular models or whether they are a general feature of cardiac systems, and if so, whether they are clinically relevant, are questions that need to be addressed in future work.

Since our study examines reduction in current, it is possible to compare our results to experiment through the use of antiarrhythmic drugs ^[75]. While not always exclusive, there are drugs that block sodium channels, potassium channels, or calcium channels. Experiments involving sodium channel blockers have shown a variety of effects, both lengthening ^{[1,27,77,82,92]} and shortening ^{[1,11,57,92,93]} of APD. Our results are consistent with the few studies that have observed no change in APD with the application of sodium channel blockers ^[21,25,61] and also agree with observations of no rate dependence with sodium channel blockers ^[3,5]. Most studies find that potassium channel blockers lengthen APD ^{[2,6,27,28,41,45,72,74,84,87,89]} and calcium channel blockers shorten APD ^[15,27,57], as predicted by our models, but there are experiments that have observed other effects ^[23,27,45]. Our model predictions are also in agreement with the majority of experiments that observed increased changes in APD as BCL decreases in potassium channel blockers ^{[3,4,6,15,28,41,50,72,74,84,89]}. The mixed results of some of these studies could be due to differences in how antiarrhythmics affect different types of cells ^{[1,3,5,27,45,48,68,82]} or they could be a consequence of the fact that most channel blockers do not simply block one channel leaving others unaffected, but rather affect multiple ion channels ^{[18,37,44,49,75,83]}. The role of both of these can be addressed through a combination of modeling and further experiment.

Experiments have also measured the slope of restitution curves in the presence or absence of antiarrhythmic drugs. Some experiments using sodium channel blockers observed a flattening of the dynamic restitution curve with the addition of drug ^[39,71], although another experiment noted an increase in the slope of the DRC ^[63]. Potassium channel blockers have also been observed to lower the slope of the DRC ^[40,60,65], but also to increase the steepness of the DRC in another study ^[62]. Calcium channel blockers were observed to decrease the slope of the SRC ^[66]. Note that these studies only examined the slope of one restitution curve, but both models examined here predict different behaviours as current is reduced for different RCs. One recent study examined the effect of sodium channel blockers and potassium channel blockers on both the DRC and the SRC finding that potassium channel blockers markedly increased the slope of the DRC, but had a much smaller effect on the SRC, and that sodium channel blockers had no effect on the slope of the DRC, but reduced the slope of the SRC ^[64]. While this experiment finds that slopes of the different RCs respond differently to reduction in current, the details of these changes do not match those predicted by our models. Another drawback to the experimental studies is the limited number of doses used in the experiments. Our models predict that some of the restitution curve slopes show increasing and decreasing regions as current is reduced — experiments using more doses might be able to confirm or invalidate these predictions.

This type of study can be helpful in discerning whether observed model behavior is due to model structure or the choice of model parameters. Model structure is the choice of functions and equations used to represent various biological processes. Model structure for cardiac models is often driven by experimental measurements of different currents and typically arises from functional differences in different types of cells (atrial vs. ventricular) or from cells of different species. A large number of computational cardiac cell models now exist to model a variety of cardiac cells from different species ^[16,17]. We used two such models, finding rather good agreement of the effect of changes in current on APD and rate-dependence. This suggests a commonality in the structure of these models that leads to the same outcome. However, the two models examined here showed differences in maximum slopes of RCs as currents were blocked. A systematic study like this highlights the commonalities and differences between models that could help us understand what kind of equations lead to certain behaviors.

5.   Conclusion

This study represents one tool in the effort to develop a mathematical model accurate enough to correctly predict the effect of a new drug ^[14,90]. Note that the Fox and Bernus models in some cases make very different predictions about some dynamics as various currents are blocked. By comparing these to experiment, we can determine which model more accurately predicts experimental observations and perhaps identify which features of current equations produce the correct behaviour. This will help us improve the accuracy of models over a broader range of parameter values, an important consideration when trying to use models for simulation of antiarrhythmics. Once we know that a model can accurately predict behaviour over a range of parameters, it can be used to assess the proarrhythmic potential of a drug before it is tested in patients ^[51,73] and can be used to help optimize treatment regimens. A systematic assessment of model parameters, coupled with experimental confirmation when possible, will help develop more accurate models that can be used for personalized medicine.

Conflict of interest

All authors declare no conflicts of interest in this paper.