Source: Cummins MA*, Dalal PJ*, Bugana M, Severi S, Sobie EA. Comprehensive analyses of ventricular myocyte models identify targets exhibiting favorable rate dependence. PLOS Comput Biol. 2014 Mar; 10(3): e1003543
Reverse rate dependence is a problematic property of antiarrhythmic drugs that prolong the cardiac action potential (AP). The prolongation caused by reverse rate dependent agents is greater at slow heart rates, resulting in both reduced arrhythmia suppression at fast rates and increased arrhythmia risk at slow rates. The opposite property, forward rate dependence, would theoretically overcome these parallel problems, yet forward rate dependent (FRD) antiarrhythmics remain elusive. Moreover, there is evidence that reverse rate dependence is an intrinsic property of perturbations to the AP. We have addressed the possibility of forward rate dependence by performing a comprehensive analysis of 13 ventricular myocyte models. By simulating populations of myocytes with varying properties and analyzing population results statistically, we simultaneously predicted the rate-dependent effects of changes in multiple model parameters. An average of 40 parameters were tested in each model, and effects on AP duration were assessed at slow (0.2 Hz) and fast (2 Hz) rates. The analysis identified a variety of FRD ionic current perturbations and generated specific predictions regarding their mechanisms. For instance, an increase in L-type calcium current is FRD when this is accompanied by indirect, rate-dependent changes in slow delayed rectifier potassium current. A comparison of predictions across models identified inward rectifier potassium current and the sodium-potassium pump as the two targets most likely to produce FRD AP prolongation. Finally, a statistical analysis of results from the 13 models demonstrated that models displaying minimal rate-dependent changes in AP shape have little capacity for FRD perturbations, whereas models with large shape changes have considerable FRD potential. This can explain differences between species and between ventricular cell types. Overall, this study provides new insights, both specific and general, into the determinants of AP duration rate dependence, and illustrates a strategy for the design of potentially beneficial antiarrhythmic drugs.
To determine whether antiarrhythmic drugs with desirable rate-dependent properties could be rationally designed, we performed comprehensive and systematic computational analyses of several heart cell models. These analyses calculated the rate-dependent properties of changes in any model parameter, thereby generating simultaneously a large number of model predictions. The analyses showed that targets with favorable rate-dependent properties could indeed be identified from such comprehensive parameter sensitivity analyses. Further simulations uncovered the mechanisms underlying these behaviors. A quantitative comparison of results obtained in different models provided new insight as to why a given drug applied to different species, or to different ventricular regions, might produce different rate-dependent electrical behaviors. Overall this study shows how a comprehensive and systematic approach in computationally analyzing heart cell models can both identify novel drug targets and produce general mechanistic insight into rate-dependent alterations to cardiac electrical activity.
Here we provide the code used for one of the 13 heart cell action potential models and its analysis. This code will allow you to reproduce Figure 4B from the paper, which shows the rate dependence properties of perturbation to ion-channel parameters in this model.
The MATLAB code hund_random.m and run1sim.m is specific to the Hund and Rudy 2004 canine ventricular cell model. The PLS.m, analysis.m, PLS_nipals.m and xticklabelrotate.m scripts/functions are general and were used for all action potential models in the study. The two files, "hund_random.m" and "run1sim.m" must be in the same directory to run the model. Hund_random.m is the master script that defines the model parameters and electrical stimulation protocol. The model parameters are varied randomly from baseline values to create a "population" of models with variable characteristics; the number of variations used in this study was 300 per pacing rate. The master script calls the function run1sim, which uses MATLAB's ODE solvers to integrate the model’s differential equations. Hund_random.m plots the action potential for each model variant and calculates model outputs, e.g. the action potential duration (APD).
Note: In our study we used steady-state conditions, and this model requires 5000 beats to reach a steady state. Each 5000 beat simulation takes about 20 min to run on a desktop PC; therefore, we recommend lowering the number of stimulations (variable n_stimuli in hund_random.m) when simply exploring this model.
The script PLS.m performs a parameter sensitivity analysis on the population of AP models generated by hund_random.m using partial least squares regression with the function PLS_nipals.m (courtesy of Herve Abdi). Parameter sensitivities quantify the magnitude and direction of influence of each parameter on the APD. PLS.m generates 2 figures: a scatter-plot of regression model APD predictions against actual output values, and a bar graph of parameter sensitivities. The function xticklabelrotate.m is used to rotate the parameter names on the x-axis ticks of the bar graph.
The script analysis.m takes the parameter sensitivities from PLS.m, calculated from simulations run at fast and slow pacing rates, and calculates the rate dependence of each parameter sensitivity. In our study we used fast = 500 ms stimulation intervals and slow = 5000 ms stimulation intervals (variable stiminterval in hund_random.m). Analysis.m generates 3 figures: a bar graph of parameter sensitivities at both rates, a scatter-plot in which the y-axis is the parameter sensitivity at fast pacing, and the x-axis is the parameter sensitivity at slow pacing, and finally a bar graph of the rate dependence of each parameter sensitivity, termed BRD. The final figure corresponds to figure 4B in our publication.