Experimentally calibrated population of models predicts and explains intersubject variability in cardiac cellular electrophysiology

Abstract

Cellular and ionic causes of variability in the electrophysiological activity of hearts from individuals of the same species are unknown. However, improved understanding of this variability is key to enable prediction of the response of specific hearts to disease and therapies. Limitations of current mathematical modeling and experimental techniques hamper our ability to provide insight into variability. Here, we describe a methodology to unravel the ionic determinants of intersubject variability exhibited in experimental recordings, based on the construction and calibration of populations of models. We illustrate the methodology through its application to rabbit Purkinje preparations, because of their importance in arrhythmias and safety pharmacology assessment. We consider a set of equations describing the biophysical processes underlying rabbit Purkinje electrophysiology, and we construct a population of over 10,000 models by randomly assigning specific parameter values corresponding to ionic current conductances and kinetics. We calibrate the model population by closely comparing simulation output and experimental recordings at three pacing frequencies. We show that 213 of the 10,000 candidate models are fully consistent with the experimental dataset. Ionic properties in the 213 models cover a wide range of values, including differences up to ±100% in several conductances. Partial correlation analysis shows that particular combinations of ionic properties determine the precise shape, amplitude, and rate dependence of specific action potentials. Finally, we demonstrate that the population of models calibrated using data obtained under physiological conditions quantitatively predicts the action potential duration prolongation caused by exposure to four concentrations of the potassium channel blocker dofetilide.

Keywords: cardiac electrophysiology; computational biology; drug; mathematical modeling; systems biology.

Fig. 1.

APs obtained from experimental recordings (red; n = 12), simulations using the models found to be within experimental range (blue; n = 213), and all models considered (black; n = 10,000) at 0.2-, 1-, and 2-Hz pacing frequencies. Plots extend to 400 ms for 2 and 1 Hz and to 800 ms for 0.2 Hz. Each experimental trace shows a representative AP from experiments on isolated female rabbit Purkinje fibers.

Fig. 2.

Scatter plots showing biomarker values for all models when stimulated at 1-Hz pacing frequency. Light gray lines indicate experimental minimum/maximum ranges for each biomarker. White dots correspond to biomarker values for models accepted into the population and, therefore, within experimental range; black dots correspond to rejected models outside of experimental range for at least one biomarker at one or more pacing frequencies. Each plot shows results for a pair of biomarkers.

Fig. 3.

Histograms of the distribution of biomarker values across the population of models for 1-Hz pacing. Dashed lines indicate the experimental range used to calibrate the population of models for each biomarker at this pacing frequency.

Fig. 4.

Scatter plots illustrating the distribution of ionic properties for accepted models in the population. Each panel shows results for a pair of ionic properties (including , , , , , , , , , , and). The scale in all graphs includes ±100% variation with respect to the original value. A representative sample of possible pairings is shown.

Fig. 5.

Correlation plots showing significant PCCs between each parameter that was varied in the population and each biomarker. Coefficient values are represented using the included color bar. For each parameter–biomarker pair, the effects on the biomarker attributable to the unselected parameters are removed as part of the partial correlation method. PCCs are shown for each frequency used in our simulations (2, 1, and 0.2 Hz). Parameter–biomarker pairs with higher coefficient values displayed stronger partial correlation with each other, whereas pairs with a coefficient of 0 did not show statistically significant partial correlation (P < 10⁻⁵). and , positive and negative correlation, respectively.

Fig. 6.

Simulated AP traces obtained for three representative models accepted in the population in control conditions (blue) and following application of 0.01 μM concentration of dofetilide (red) at 2-, 1-, and 0.2-Hz pacing frequencies. This is the concentration closest to the experimentally determined IC₅₀ (therapeutic dose) for dofetilide (0.0124 μM). Plots extend to 500 ms for 2 and 1 Hz and 1,000 ms for 0.2 Hz. Line style indicates which of the control and dofetilide traces correspond to each model.

Fig. 7.

(A) Ranges of prolongation (ΔAPD) caused by four concentrations of dofetilide using the models in the calibrated population. Dots indicate values of ΔAPD independently obtained in five experiments using rabbit Purkinje fiber preparations. (B) Histograms of dofetilide-induced ΔAPD range across sets of 5 models randomly sampled from the calibrated population. A total of 100,000 different sets of five models were used. The range of ΔAPD was calculated as maximum value of ΔAPD − minimum value of ΔAPD, for each set of five models. For each dofetilide concentration, the mean value of the ΔAPD range across the 100,000 samples of five models is shown by a solid blue line, and the ΔAPD range from our experimental data is shown by the dashed red line.

Fig. 8.

(A) PCCs for APD prolongation (ΔAPD) caused by block from a 0.01 μM dofetilide concentration, at each of the three pacing frequencies. ΔAPD was correlated against each of the 12 parameters that were varied to create the population of models, each time controlling for the other 11 parameters as part of the partial correlation process. Three models with outlying values of ΔAPD at 0.2 Hz (ΔAPD = 318, 364, and 395 ms, compared with the rest of the model population’s sample mean of 69 ms ± 26 ms) were excluded from the analysis as they dominate the other models when calculating the PCCs at that frequency. Solid bars denote significant correlations with P < 10⁻⁵, and empty bars indicate correlations with P > 10⁻⁵. (B) Scatter plot of control values of each accepted model at 1 Hz against ΔAPD attributable to block from 0.01 μM concentration of dofetilide.

Fig. 9.

Schematic for the rabbit Purkinje cell model used in our study. Ionic currents included in the model and the calcium handling subsystem are shown. Arrows within the cell represent calcium transport between compartments. PMCA, plasma membrane Ca²⁺-ATPase; SERCA, sarco/endoplasmic reticulum Ca²⁺-ATPase.

Fig. 10.

Representation of each biomarker measured in this study, calculated from a simulation at 1 Hz using a model from the population of models. The maximum upstroke velocity was calculated as the maximum value of the gradient of the membrane potential against time recorded before the point where the Vm Peak occurs.