Sensitivity of Rabbit Ventricular Action Potential and Ca²⁺ Dynamics to Small Variations in Membrane Currents and Ion Diffusion Coefficients

Abstract

Little is known about how small variations in ionic currents and Ca²⁺ and Na⁺ diffusion coefficients impact action potential and Ca²⁺ dynamics in rabbit ventricular myocytes. We applied sensitivity analysis to quantify the sensitivity of Shannon et al. model (Biophys. J., 2004) to 5%–10% changes in currents conductance, channels distribution, and ion diffusion in rabbit ventricular cells. We found that action potential duration and Ca²⁺ peaks are highly sensitive to 10% increase in L-type Ca²⁺ current; moderately influenced by 10% increase in Na⁺-Ca²⁺ exchanger, Na⁺-K⁺ pump, rapid delayed and slow transient outward K⁺ currents, and Cl⁻ background current; insensitive to 10% increases in all other ionic currents and sarcoplasmic reticulum Ca²⁺ fluxes. Cell electrical activity is strongly affected by 5% shift of L-type Ca²⁺ channels and Na⁺-Ca²⁺ exchanger in between junctional and submembrane spaces while Ca²⁺-activated Cl⁻-channel redistribution has the modest effect. Small changes in submembrane and cytosolic diffusion coefficients for Ca²⁺, but not in Na⁺ transfer, may alter notably myocyte contraction. Our studies highlight the need for more precise measurements and further extending and testing of the Shannon et al. model. Our results demonstrate usefulness of sensitivity analysis to identify specific knowledge gaps and controversies related to ventricular cell electrophysiology and Ca²⁺ signaling.

1. Introduction

Several ionic models have been developed to investigate the subcellular mechanisms regulating excitation-contraction coupling (ECC) in rabbit ventricular cardiomyocytes [1–7]. In 2004, Shannon and colleagues published a detailed model for Ca²⁺ handling and ionic current that accurately represents sarcoplasmic reticulum (SR) Ca²⁺-dependent release and simulates basic ECC phenomena. This model was the first to include (1) the subsarcolemmal Ca²⁺ compartment to the other two commonly formulated cytosolic Ca²⁺ compartments (junctional and bulk, [8]), (2) the variations in the locations of ion transporters throughout the cell surface membrane, (3) Ca²⁺ and Na⁺ transport between the subcellular compartments, and (4) Na⁺ buffering inside. Latest studies extended further the Shannon et al. ionic model in rabbit ventricular cells. Mahajan and colleagues modified L-type Ca²⁺ current and Ca²⁺ cycling formulations, based on new experimental patch-clamp data, and used the updated model to investigate the mechanisms regulating ventricular tachycardia and fibrillation [4, 5]. Morotti et al. improved Mahajan’s et al. model of rabbit ventricular 𝐼_Ca and examined the relative contributions of voltage- and Ca²⁺-dependent inactivation to total current inactivation [7]. Aslanidi et al. determined the functional role of each ionic channel current in modulating the action potential (AP) shape from different locations in rabbit ventricular wall [6]. Recently, using parameter sensitive analysis, Romero et al. [9] identified also key subcellular factors involved in the electrical signal generation and Ca²⁺ signal regulation in two rabbit cell models [1, 4]. Systematically characterizing AP properties, Ca²⁺ and Na⁺ dynamics, and their rate dependence on ±15%, ±30, and ±45% variations in the main transmembrane currents’ conductance they found that APD is significantly modified by most repolarization currents and that steady-state Ca²⁺ levels are strongly dependent on 𝐼_Ca, 𝐼_NCX, and 𝐼_NaK currents. Important limitation of these “common pool models”, however, is that the subcellular compartment geometries (junctional cleft, submembrane space, cytosol, and SR) were simplified to allow prediction of total Ca²⁺ concentration ([Ca](t)) in the compartment of interest only.

A few reaction-diffusion models that seek to incorporate the effects of idealistic or more realistic subcellular and whole-cell geometries and to compute the spatial distributions of [Ca] in cardiac, skeletal or smooth muscle cells have been published also [10–19]. Recently, the contributions of structural heterogeneities to local and global Ca²⁺ signals in rabbit ventricular myocytes with realistic 𝑡-tubule microanatomy have been investigated [20]. Limitations of these spatial models, due to complexities, are that the cell electrophysiology was simplified, the AP dynamics was not reproduced, and the effects of intracellular Na²⁺ diffusion and buffering on ECC were not investigated.

Taken together the published computational models have provided important insights into the underlying mechanisms regulating intracellular Ca²⁺ dynamics and electrophysiological behavior in rabbit ventricular cells under control and certain pathological conditions. It remains, however, poorly understood how small variations in ion currents conductance and distribution and small changes in Ca²⁺ and Na⁺ transport between cellular subdomains affect important cellular biomarkers (i.e., AP shape, [Ca]_𝑖, [Ca]_SL, and [Ca]_jct). In this study, we applied sensitivity analysis and used parameter estimation tools to investigate these questions. The sensitivity of each biomarker to 5%–10% changes in ionic currents properties and Ca²⁺ and Na⁺ diffusion coefficients were quantified using the Shannon et al. model. Validation of simulation results was performed by a comparison to experimental data in rabbit ventricular myocytes. This study identifies new experimental data required for better understanding of rabbit’s electrophysiology and demonstrates the importance of sensitivity analysis as a powerful method for systematic and in-depth validation of AP models.

2. Methods

2.1. Electrophysiological Model in Rabbit Ventricular Myocytes

A diagram describing the arrangement of subcellular compartments, ionic currents and pumps, and intracellular Ca²⁺ fluxes included in the Shannon et al. electrophysiological model is shown in Figure 1.

Figure 1 

Schematic diagram illustrating the four subcellular compartments, electrophysiology, and Ca2+ dynamics in the Shannon et al. model in rabbit ventricular myocytes [1]. See Glossary and Tables 1–3 for the notations of the parameters used throughout the study.

The Shannon et al. model was adopted without changes and recent MatLab version of the code has been uploaded (http://somapp.ucdavis.edu/Pharmacology/bers/). The model behavior at 1 Hz pacing frequency is consistent with experimental measurements for AP-shape and global cytosolic Ca²⁺ concentrations ([Ca]_𝑖) changes under control conditions [1].

2.2. Sensitivity Analysis Tools

Recently, several sensitivity analysis and parameter estimation tools and methodologies have been developed with a specific aim to assess the robustness of the cardiac cell models’ parameters, including (1) the parameter variability analysis which examine the effects of varying one parameter in time [21], (2) techniques which assess the significance of parameters using multivariable regression over randomized set of input parameters values [22], and (3) global sensitivity analysis tool which implements the Morris screening algorithm [23] within the Nimrod platform [24]. In this study, we explored two methodologies [22, 23] and compared each method potential in evaluating parameter-output correlations in the Shannon et al. electrophysiological model.

2.2.1. Partial Least Squares (PLS) Regression

PLS is a statistical method that creates a linear regression model to predict a set of dependent variables (i.e., model outputs) from a set of independent variables (i.e., model inputs) [22, 25]. PLS is particularly useful when inputs are highly correlated (i.e., multicollinear) and when there are more inputs than outputs. In regression analysis, if two highly correlated effects are measured, the more dominant effect will falsely enhance the weaker effect, leading to inaccurate sensitivity information, and misinform truly significant parameter effects [25]. Since PLS retains good predictive power when working with multicollinear data, myocyte models containing multiple correlated variables can be analyzed using PLS.

2.2.2. Nimrod Tool Family

Nimrod uses a simple declarative modeling language to define parametric experiments, and it automates the tasks of formulating, running, monitoring, and collecting results from multiple individual experiments. Nimrod is not a single tool. It incorporates Nimrod/G, component that distributes computations to the recourses [24, 26]; Nimrod/O, a component that searches for “good” solutions using nonlinear optimization algorithms [27]; Nimrod/E, component that helps evaluating which parameter settings are important using the experimental design [28]. Nimrod/G is designed to assist scientists in performing studies, using concurrent execution on a cluster of processors or the resources of a computational grid. The user prepares a “plan file” which specifies the experiment [24]. A plan file is expanded into a “run file” which lists the appropriate parameter combinations required for execution. The Nimrod/G core runs on a processor called the “root node” whereas individual jobs are executed on “remote nodes.” Nimrod/G handles the file transfers required for the remote nodes, execution of computational tasks, and transfer of results back to the root node. The number of concurrent jobs is limited only by the number of processors available. Thus the user may achieve high concurrency without modifying the executables and without concern for grid specific details. Nimrod/O optimizes the numerical outputs of computational models. It provides a range of optimization algorithms and leverages Nimrod/G to perform batches of concurrent evaluations. The user prepares a “schedule file,” which, like the Nimrod/G plan file, specifies the parameters, their ranges, and the tasks required for execution of the model. But it also specifies the optimization algorithms to be used and the settings for the algorithms. Nimrod/O performs multiple searches for an algorithm and multiple algorithms, all in parallel.

2.2.3. Nimrod/E and Fractional-Factorial Design (FFD)

It is often of interest for researchers to identify variable interactions because they may provide physiological insights into cell dynamics. While colinearity in multiple linear regressions is problematic and can lead to inaccurate predictions, Nimrod/E is able to circumvent this problem and successfully identify both main effects and multifactor interactions [28]. Similar to PLSR, FFD attempts to create a linear model to explain the data generated from the model; it measures and ranks main effects and two-level interactions of input variables with outputs via experimental design. However, the method chooses only those jobs that generate the most significant results by ignoring high-level interactions with little influence, thereby saving a lot of computation time. From the resulting linear model effect coefficients can also be extracted that tells the correlation between parameters and model outputs. For more complete view of successes and failures in modeling pursuits of Nimrod toolsets, refer to Abramson and colleagues papers [24, 26–28].

2.3. Parameter Sensitivity Analysis

To perform sensitivity analysis, the Shannon et al. MatLab code was modified to extract the desired model outputs. Under basal conditions or by varying one parameter at a time AP and Ca²⁺ transients appeared stable (after 6–8 s until 8 min) at stimulus frequency of 1 Hz. For this reason we recorded and analyzed AP and Ca²⁺ signals ([Ca]_𝑖, [Ca]_SL, [Ca]_jct) at 9-10 s of every simulation. In all numerical experiments (0≤𝑡≤8 min) [K]_𝑖 remained unchanged while the variations in [Na]_𝑖 were small (ranging from 20 μM to 100 μM) with insignificant effects on calculated AP and Ca²⁺ traces at 6–8 s. The APD₆₀ and delta [Ca]_𝑗 (Δ[Ca]_𝑗, 𝑗=𝑖, SL, jct) quantities were used to analyze the model sensitivity. Δ[Ca]_𝑗 is defined as the difference between peak Ca²⁺ concentration in the particular compartment and the diastolic Ca²⁺ concentration of 0.1 μM. The modified code was uploaded onto Nimrod/E toolset and was executed repeatedly using permutations of parameter sets that were ±10% or ±5% of default values (see the appendix, Tables 1–3). Then Nimrod/E fractional-factorial analysis and PLS were used to analyze the effects of single parameter changes. In addition, Nimrod/E allowed examining the two-level parameter interactions on model outputs. In the Results section the sensitivities of selected biomarkers to variations of parameters are displayed graphically as either “Lenth plots” in Nimrod/E, or as “bar plots” in PLS. Before performing PLS regression analysis to obtain sensitivity values, 𝑧-scores of input and output matrices were calculated using MatLab’s 𝑧-score function; that is, 𝑧=(𝑥−mean(𝑥))/std(𝑥). Each 𝑧-score value was computed using the mean and standard deviations along each column of the matrices. The columns of matrices have mean zero and standard deviation one. The PLS regression coefficients, or sensitivity values (see Figures 2, 4, and 6), indicate how changes in input parameters lead to changes in outputs. Examining these numbers allows for an assessment of the relative contributions of the various parameters. The sensitivity values in the “bar plots” could be interpreted quantitatively as follows. Because input and output matrices are mean-centered and normalized to standard deviations computed column by column, each sensitivity value is defined relatively to the relevant deviations. For instance, if the regression coefficient for the input 𝑃_CaL and the output APD₆₀ is 0.5, then when 𝑃_CaL is one standard deviation greater than the mean, APD₆₀ will increase by half a standard deviation. Conversely, if the value is −0.5, then when 𝑃_CaL is one standard deviation greater than the mean, APD₆₀ will decrease by half a standard deviation, and vice versa. The estimate values in the “Lenth plots” (𝑦-axis) provide a qualitative overview of the inputs’ relative effects on outputs and offer a comparison to the PLS-bar graphs. Changing the number of parameters studied slightly varied the magnitude of “regression coefficients” of bar plot and “estimate” of Lenth plot, but the relative effects and overall parameter-to-output relationships remained constant. The effects of sample sizes on model predictive efficacy by randomly picking subsets from sample pool, ranging from 15, 20, 50, 100, 200, and 500 to 1000 trials was examined. The adjusted 𝑅² values (coefficient of determination, quantifying the explanatory capacity of the regression analysis) were calculated and averaged. At the low end of 15 samples, the regression model explained 90.8 ± 2.4% of the variance. At the high end of 1000 samples, the model explained 96.4 ± 0.3% of the variance. This relatively low decrease in prediction efficacy suggests that sample sizes do not significantly change the qualitative information obtained from the experiments. Considering this statistical analysis, the results were not included.

Table 1 

Maximal conductance varied in the model.

Table 2 

Fractions of total “𝑗” flux varied in the model.

Table 3 

Diffusion coefficients varied in the model.

Figure 2 

Parameter sensitivity analysis of APD60, Δ[Ca]𝑖, Δ[Ca]SL, and Δ[Ca]jct sensitivities toward 10% increases in maximal transporters’ conductance and permeability. Positive value indicates direct relationship between parameter and output (e.g., enhanced L-type Ca2+ permeability prolongs APD60), and negative value indicates indirect relationship between parameter and output. Nimrod/E-Lenth plots (left column), PLS-bar plots (right column). Nimrod/E analysis produced similar qualitative results as PLS analysis.

3. Results

3.1. Effects of 10% Changes in Channels Conductances on AP Morphology and Intracellular Ca²⁺ Dynamics

Figure 2 shows the effects of 10% changes in 19 default ion channel conductances and permeabilities (see the appendix, Table 1). Both PLS and Nimrod/E methods yielded consistent results with respect to single parameter’s changes on the model outputs (APD₆₀, Δ[Ca]_𝑖, Δ[Ca]_SL, Δ[Ca]_jct).

The increase in L-type Ca²⁺ channel permeability (𝑃_CaL) had the most pronounced effects on APD₆₀ and Ca²⁺ signals, by increasing action potential duration and enhancing maximum Ca²⁺ peaks. The effects of 𝐺_Kr, 𝐺_to, and 𝐾_NCX, while less pronounced than 𝑃_CaL, were still significant. Figure 2 also shows that 10% increases in 𝐺_Kr and 𝐺_to decreased APD₆₀ and Ca²⁺ peaks while 10% increase in 𝐾_NCX decreased Ca²⁺ peaks but increased APD₆₀. The increase in maximal NaK pump rate (𝐾_NaK) had still detectable but minor effect, by decreasing both APD₆₀ and Ca²⁺ peaks. The changes in SR parameters (𝐾_rel, 𝐾_up, and 𝐾_leak) by 10% had negligible effects on APD₆₀ and Δ[Ca]_𝑖 whereas the submembrane and junctional Ca²⁺ peaks (e.g., Δ[Ca]_SL and Δ[Ca]_jct) were most affected. Figure 2 shows that 10% increases in 𝐾_rel and 𝐾_up increased Δ[Ca]_SL and Δ[Ca]_jct while 𝐾_leak had the opposite effect. Interestingly, 10% increase in background Cl⁻ conductance (𝐺_ClBk) showed similar magnitude as 𝐺_Kr and 𝐺_to and demonstrated the same effects in decreasing APD₆₀ and Ca²⁺ peaks. Furthermore, Nimrod/E had the advantage in predicting how 10% increases in two-parameter group can affect the model outputs. Nimrod/E predicts here that 𝑃_CaL and K_NCX combined, as well all other two-parameter group combinations (not shown in Figure 2), had slight or no effect on the selected cellular biomarkers (APD₆₀, Δ[Ca]_𝑖, Δ[Ca]_SL, and Δ[Ca]_jct). The plots in Figure 3 confirm further the specific, quantitative predictions made by PLS and Nimrod/E for the effects of 10% increases in 𝑃_CaL, 𝐾_NCX, 𝐺_Kr, 𝐺_to, 𝐾_NaK, and 𝐺_ClBk on AP morphology and Ca²⁺ transients.

(a)

(b)

(c)

(d)

Figure 3 

Verification of ion conductance sensitivity predictions. Steady-state AP and Ca2+ transients (recorded 9-10 s) in each of the three non-SR compartments in response to 1 Hz stimulus. The predicted changes in APD60 and Ca2+ peaks (see Insets also) are consistent with the generated by Nimrod/E and PLS outputs with respect to 10% increases in transporters’ conductance. Control conductance values (black solid line) and 10% increases in 𝐾NCX (black dashed line), 𝑃CaL (black dash-dot line), 𝐺Kr (gray dash-dot line), 𝐺to (gray solid line), 𝐾NaK(gray dotted line), and 𝐺ClBk (black dotted line).

Figure 4 

Parameter sensitivity analysis of APD60, Δ[Ca]𝑖, Δ[Ca]SL, and Δ[Ca]jct sensitivities toward 5% increases in ion transporter distributions in the junctional cleft. Transporters’ distribution between sarcolemma and junctional cleft summed up to 100%. Increasing fraction of total 𝑗-current in the junction by 5% is analogous to decreasing fraction of total 𝑗-current in the subsarcolemma by 45% for 𝐼Ca and by 40-41% for all other currents. Positive value indicates direct relationship between distribution and output (e.g., increased density of 𝐼Ca in the junctional cleft shortens APD60) and negative value indicates indirect relationship between parameter and output. Nimrod/E-Lenth plots (left column), PLS-bar plots (right column). Nimrod/E analysis produced similar qualitative results as PLS analysis.

3.2. Effects of Changes in Membrane Transporter Distributions on AP Morphology and Intracellular Ca²⁺ Dynamics

In Shannon et al. model, cell was separated into four lumped compartments (see Figure 1): the junction junctional cleft (0.077% of total cell volume assuming 11% of the surface membrane junctional), the subsarcolemmal space (2% of total cell volume assuming 89% of the membrane nonjunctional) the bulk cytosol space (65% of total cell volume with remainder of the volume accounted for by mitochondria) and the SR (3.5% of total cell volume). The L-type Ca²⁺ channels were assumed concentrated within the junctional membrane such that 90% of total number were located there while all other membrane transporters (Na⁺ channels, Na⁺ leak, Na⁺/K⁺ pump, slow delayed rectified K⁺ channel, Ca²⁺ activated Cl⁻ channel, Ca²⁺ leak, Na⁺/Ca²⁺ exchanger, and sarcolemmal Ca²⁺ pump) were considered evenly distributed across the sarcolemma with 11% in the junctional cleft and 89% in the subsarcolemmal compartment [1].

To examine how the changes in ion transporter distributions affect the selected cellular biomarkers (APD₆₀, Δ[Ca]_𝑖, Δ[Ca]_SL, and Δ[Ca]_jct), we performed systematic analysis, increasing by 5% the basic 𝐹𝑥_𝑗(jct) while decreasing 𝐹𝑥_𝑗(SL) to keep total number of transporter protein complexes within the sarcolemma unchanged; that is, 𝐹𝑥_𝑗(jct)+𝐹𝑥_𝑗(SL)=1 (the appendix, Table 2). Figure 4 shows that PLS and Nimrod/E analysis again showed consistent results. The 5% increase in L-type Ca²⁺ channels fraction in the junctional cleft (𝐹𝑥_CaL(jct)=0.945) had the most pronounced effect on the model outputs, decreasing APD₆₀ and Ca²⁺ peaks. The results also suggest that 5% increase in the junctional Na⁺-Ca²⁺ current (𝐹𝑥_NCX(jct)=0.1155) prolonged APD₆₀ while slightly affected Ca²⁺ peaks in the junctional cleft and subsarcolemmal and bulk cytosol compartments. The 5% increase in the junctional 𝐼_Cl(Ca) (𝐹𝑥_Cl(Ca)(jct)=0.1155) had similar effect on all model outputs as the increase in 𝐹𝑥_CaL(jct) but with much smaller magnitude. Results in Lenth plot (Figure 4) for the four model outputs reveal that “𝐹𝑥_NCX(jct)𝐹𝑥_CaL(jct) combination” as well all other two-parameter combinations had insignificant impact (not shown in Figure 4). Additional analysis of the data has been done by increasing each 𝐹𝑥_𝑗(jct) by 5% and decreasing 𝐹𝑥_𝑗(SL) (𝑗 = CaL, NCX, Cl(Ca)) while keeping all other model parameters constant. The results in Figure 5 confirm the predictions made by both PLS and Nimrod/E methods.

(a)

(b)

(c)

(d)

Figure 5 

Verification of transporter distribution sensitivity predictions. Steady-state model outputs (recorded 9-10 s) in response to 1 Hz periodic pulse. The predicted alterations in APD60 and Ca2+ peaks (see Insets also) are consistent with the generated by Nimrod/E and PLS results with respect to 5% increases in transporters’ fraction in the junctional cleft. Default 𝑗-currents’ fractions in junctional cleft (black solid line) and 5% increases in 𝐹𝑥CaL(jct) (black dashed line), 𝐹𝑥Cl(Ca)(jct) (grey dotted line), 𝐹𝑥NCX(jct) (black dash-dot line).

Figure 6 

Parameter sensitivity analysis of APD60, Δ[Ca]𝑖, Δ[Ca]SL, and Δ[Ca]jct sensitivities toward 10% increases in Ca2+ and Na+ diffusion coefficients between subcellular compartments. Positive value indicates direct relationship between diffusion coefficient and output (e.g., increased diffusion coefficient between junctional cleft and sarcolemma prolongs APD60) and negative sensitivity value indicates indirect relationship between parameter and output. Nimrod/E-Lenth plots (left column), PLS-bar plots (right column). Nimrod/E analysis produced similar qualitative results as PLS analysis.

3.3. Effects of 10% Changes in Diffusion Constants on AP Morphology and Intracellular Ca²⁺ Dynamics

The effects of changes in default parameters describing Ca²⁺ and Na⁺ diffusion between junctional and subsarcolemmal compartments and between subsarcolemmal and bulk compartments were also analyzed in this study (the appendix, Table 3). Figure 6 shows that 10% increases in basal 𝐷_Ca(jct-SL) altered all model outputs, prolonging APD₆₀ and increasing Δ[Ca]_𝑖 and Δ[Ca]_SL but decreasing Δ[Ca]_jct. Moreover 10% increases in basal 𝐷_Ca(SL-𝑖) value shortened APD₆₀ and increased Δ[Ca]_𝑖, Δ[Ca]_jct, and Δ[Ca]_SL. Interestingly, 10% increases in basal 𝐷_Na(jct-SL) and 𝐷_Na(SL-𝑖) values in all three compartments had no effect on APD₆₀ and Ca²⁺ peaks. Nimrod/E analysis also suggests that two-level interaction between 𝐷_Ca(jct-SL) and 𝐷_Ca(SL-𝑖) had a minor effect on increasing Δ[Ca]_𝑖 while APD₆₀, Δ[Ca]_SL, and Δ[Ca]_jct were unaffected. Figure 7 additionally confirms the PLS and Nimrod/E predictions for the effects of 10% increases in 𝐷_Ca(jct-SL) and 𝐷_Ca(SL-𝑖) on AP morphology and Ca²⁺ transients.

(a)

(b)

(c)

(d)

Figure 7 

Verification of Ca2+ and Na+ diffusion sensitivity predictions. Steady-state AP and Ca2+ transients (recorded 9-10 s) in three non-SR compartments in response to 1 Hz stimulus. The relative changes in APD60 and Ca2+ peaks (see Insets also) are consistent with predictions generated by Nimrod/E and PLS with respect to 10% increases in Ca2+ and Na+ diffusion constants. Control diffusion coefficients (black solid line) and 10% increases in 𝐷Ca(jct-SL) (black dashed line), 𝐷Ca(SL-𝑖) (black dash-dot line), 𝐷Na(jct-SL) (gray dotted lines), 𝐷Na(SL-𝑖) (light gray dotted lines).

4. Discussion

The main goal of this study was to perform a systematic investigation on how small variations in ionic currents and Ca²⁺ and Na⁺ diffusion coefficients modulate ventricular myocyte electrophysiology and Ca²⁺ dynamics in rabbit ventricular myocytes. To determine the most influential parameters controlling underlying subcellular mechanisms, we applied and compared the outcomes from two sensitivity analysis toolkits (PLS and Nimrod). The sensitivity of selected biomarkers (APD₆₀, Δ[Ca]_𝑖, Δ[Ca]_SL, and Δ[Ca]_jct) to the alterations in ionic currents’ properties and intracellular ionic transport was quantified using the Shannon et al. ionic model in rabbit ventricular cells, the predictive power of which has been demonstrated previously [1, 2].

4.1. Effects of Small Variations in Ionic Currents on AP Duration and Ca²⁺ Transients

Using parameter sensitive analysis, Romero et al. quantified the sensitivity of AP properties, Ca²⁺ and Na⁺ dynamics, and their rate dependence with respect to ±15%, ±30%, ±45%, and ±100% variations in the main transmembrane currents’ conductance [9]. It remains, however, poorly understood how small variations (±5, ±10%) in ion currents conductance affect cellular biomarkers. Our next studies contribute to filling up this knowledge gap further. Analysis suggests that the selected biomarkers were strongly dependent on 𝐼_Ca current properties. Results indicate that 10% increase in maximal permeability of L-type Ca²⁺ channel (𝑃_CaL), with respect to other currents’ conductance changes, most notably affected APD₆₀ and Ca²⁺ transients’ peaks at 1 Hz. Our predictions for the effects of 𝑃_CaL on APD and cytosolic Ca²⁺ peak are also consistent with reported data in rabbit ventricular myocytes after applications of 𝐼_Ca blockers [5, 9, 29] or in cells with high L-type Ca²⁺ channel density [30]. The effects of 10% increase in 𝐺_Kr, 𝐺_to, 𝐾_NCX, or 𝐾_NaK on cellular biomarkers, while less pronounced than 𝑃_CaL, were still detectible. Our predictions for 𝐾_NCX increase at 1 Hz are in qualitative agreement with data in ferret myocytes [31] and data supporting APD prolongation following 𝐼_NCX block reported in some experiments in rabbit cardiomyocytes [32, 33]. Several studies in rabbit ventricular cells report also a strong effect of 𝐼_NCX in modulating systolic [Ca]_𝑖 transient [31–34] while others report only a moderate effect [35], similar to our predictions for Δ[Ca]_𝑖. Results demonstrate also that 10% increase in 𝐾_NCX slightly decreased [Ca]_jct and [Ca]_SL peaks. Increase in 𝐺_Kr and 𝐺_to at 1 Hz resulted in APD shortening, which is consistent with experiments in rabbit ventricular myocytes [36, 37]. Important role of 𝐼_to in modulating Ca²⁺ transients through its effects on 𝐼_Ca and 𝐼_NCX during AP repolarization period has been suggested [38, 39]. Slight effects, however, due to 10% increase in 𝐺_to at 1 Hz on Ca²⁺ peaks are predicted here. Furthermore, at stimulus frequency of 2.5 Hz, Romero et al. report 𝐼_NaK (±30% changes or severe 𝐼_NaK block; 0–8 min) as important influent current in most simulated cellular biomarkers, including APD₆₀ [9]. Our studies quantitatively confirmed Romero et al. findings at 2.5 Hz (data not shown). At 1 Hz, however, the effect of ±30% changes in 𝐾_NaK (increment 5%, 0–8 min) on APD₆₀ was less pronounced. The main reason for this observation was that the predicted changes in [Na]_𝑖 were smaller at 1 Hz (20–100 μM) versus 1–1.5 mM at 2.5 Hz, in this way affecting differently 𝐼_NaK and 𝐼_NCX activities. These results indicate that the stimulus frequency is important factor regulating the electrophysiological biomarkers. Accordingly, because the Shannon et al. model mimics rabbit electrophysiology more accurately at normal pacing rates of 1 Hz [9], the model predictions for the effects of changes of ion current conductivities at faster rates should be considered with caution. Interestingly, 10% increase in Cl⁻ background current (𝐼_ClBk), which is a new feature of the Shannon et al. model, showed a negative relationships with four outputs (APD₆₀, Δ[Ca]_𝑖, Δ[Ca]_SL, and Δ[Ca]_jct) of approximately similar magnitude as for 𝐺_to or 𝐺_Kr. Experimental data reporting the effects of 𝐼_ClBk on rabbit ECC are missing in the literature however. Another interesting finding is the relatively small sensitivity of the biomarkers to 10% increases in 𝐺_Na, 𝐺_Na,b, 𝑃_CaK, 𝑃_CaNa, 𝐺_Ca,b, 𝐺_to,f, 𝐺_Ks, 𝐺_K1, 𝐺_Cl(Ca), 𝐾_rel, 𝐾_leak, 𝐾_up, and 𝐾_pCa factors at 1 Hz. Thus, to further test the predictive power of the Shannon et al. model we examined the effects of severe currents’ and SR-fluxes’ block (up to 100%) (data not shown). In qualitative agreement with experiment [9] our studies show that (1) APD₆₀ was highly sensitive to 100% 𝐼_Na, 𝐼_K1, 𝐾_rel, and 𝐾_up and block and less sensitive to 𝐼_CaK, 𝐼_Cl(Ca), 𝐼_to,f, and 𝐼_pCa inhibition while 100% block of 𝐼_Na,b, 𝐼_CaNa, 𝐼_Ca,b, and 𝐼_Ks currents and SR Ca²⁺ leak (𝐾_leak) had small or no effect; (2) steady-state Ca²⁺ peaks were strongly affected when 𝐼_Na, 𝐼_Ca,b, and SR Ca²⁺ fluxes (𝐾_rel, 𝐾_leak, and 𝐾_up) were inhibited whereas 100% block 𝐼_pCa, 𝐼_Cl(Ca), 𝐼_CaK, 𝐼_K1, and 𝐼_to,f had less pronounced effects. New measurements in rabbit ventricular myocytes are needed to define the relative importance of 𝐺_Na, 𝐺_Na,b, 𝑃_CaK, 𝑃_CaNa, 𝐺_Ca,b, 𝐺_to,f, 𝐺_Ks, 𝐺_K1, 𝐺_Cl(Ca), 𝐾_rel, 𝐾_leak, 𝐾_up, and 𝐾_pCa in determining physiological responses due to the small changes in these factors at 1 Hz and to test further the capabilities and limitations of the Shannon et al. model. Using Nimrod/E, we also examined how simultaneous changes in two transmembrane currents affect APD and Ca²⁺ dynamics. Our analysis suggests that 10% increase in both 𝑃_CaL and 𝐾_NCX had minor effect at 1 Hz. Additionally, Nimrod/E showed that simultaneous increases in control 𝐼_Ca and K⁺ related currents had slight or no effect (data not shown). Furthermore, the predicted effects of ±5% changes in channels conductance were identical to those with ±10% variations (data not shown). Such insights, coming from mathematical models, might be useful in searching new avenues for development of methodologies to predict drug action effects on behavior of cardiac cells. For example, most 𝐼_Ca inhibitors suppress partially K⁺-related currents also; therefore the balance of these currents is critical for predicting 𝐼_Ca-related drug action.

4.2. Effects of Variations in Membrane Transporters Distribution on AP Duration and Ca²⁺ Transients

Experimental studies have demonstrated that in rabbit ventricular myocytes marked variations in the distribution of ion transports along the surface membrane probably exist [1, 40–42]. Using reaction-diffusion models, recently we demonstrated that the subcellular Ca²⁺ signals are highly sensitive to 𝐼_Ca and 𝐼_NCX fluxes distributions via the sarcolemma [13, 16, 19, 20]. In these spatial models, however, the cell electrophysiology was simplified; equations describing 𝐼_Ca, 𝐼_NCX, 𝐼_pCa, and 𝐼_Ca,b properties were only included; and AP dynamics were not reproduced. Additionally, small variations in L-type Ca²⁺ channels density in junctional cleft have been suggested to lead to large changes in control SR Ca²⁺ release, subcellular Ca²⁺ transients, and AP waveform [4, 5]. Thus, we used PLS and Nimrod/E to gain further insights into how ±5% changes in L-type and other ion transporters’ distributions influence control AP and Ca²⁺ transients. We evaluated the biomarkers’ sensitivity to 5% increases in the fraction of total ionic current in the junctional cleft (𝐹𝑥_𝑗(jct)) assuming total number of transporters within the sarcolemma unchanged; that is, 𝐹𝑥_𝑗(jct)+𝐹𝑥_𝑗(SL)=1. Our results yielded three channel distributions (𝐹𝑥_CaL(jct), 𝑥_NCX(jct), and 𝐹𝑥_Cl(Ca)(jct)) that significantly affect the model outcomes of interest (i.e., APD₆₀, Δ[Ca]_𝑖, Δ[Ca]_SL, and Δ[Ca]_jct). Increasing by 5% the density of L-type Ca²⁺ channels in the junctional cleft (while decreasing by 45% 𝐼_Ca density in SL compartment to keep 𝐹𝑥_CaL(jct)+𝐹𝑥_CaL(SL)=1) had the strongest effect on both APD₆₀ and Ca²⁺peaks, decreasing APD₆₀ and [Ca]_𝑖, [Ca]_SL, and [Ca]_jct peaks. The reason for the predicted decrease in [Ca]_jct was the strongly increased rate of Ca²⁺ transfer from the junctional cleft to the submembrane space evoked by the 5% increase of 𝐹𝑥_CaL(jct) and reduced 𝐼_Ca (~45%) into the SL compartment. The predicted decreases into [Ca]_SL peak and APD₆₀ were also consequence of 𝐹𝑥_CaL(SL) reduction to 45% though Ca²⁺ flux from the cleft was increased. The drop in [Ca]_SL evoked decrease in the rate of Ca²⁺ transfer from the SL to cytosolic compartment thus causing drop in [Ca]_𝑖 peak as well. Our studies also showed that the effects of 𝐼_NCX redistribution were quite different. The 5% increase in 𝐹𝑥_NCX(jct) (while decreasing by 40-41% 𝐼_NCX density in SL compartment to keep 𝐹𝑥_NCX(jct)+𝐹𝑥_NCX(SL)=1) prolonged APD₆₀to approximately the same magnitude as for 𝐼_Ca and slightly decreased the control Ca²⁺ peaks. Significant effects due to 𝐼_Ca redistribution and negligible effects due to 𝐼_NCX redistribution on the local and global Ca²⁺ transients’ peaks were also found using our reaction-diffusion model in rabbit ventricular cells [20]. Calcium activated chloride current (𝐼_Cl(Ca)) has been extensively studied [43–45]. Experimental data in rabbit ventricular myocytes suggest that 𝐼_Cl(Ca) is strongly temperature dependent (very small at room temperature but substantial at 35°C–37°C); 𝐼_Cl(Ca) can be eliminated by blocking 𝐼_Ca; and 𝐼_Cl(Ca) may normally facilitate rapid repolarization when SR Ca²⁺ load and release are high. To bring further insights into the effects of 𝐼_Cl(Ca) properties on ECC, we examined how 5% increase in 𝐼_Cl(Ca) density in the junctional cleft (assuming 𝐹𝑥_Cl(Ca)(jct)+𝐹𝑥_Cl(Ca)(SL)=1) may affect the cell electrical activity and Ca²⁺ signaling. Interestingly, although much less pronounced than 𝐹𝑥_CaL(jct) variations, changes in 𝐹𝑥_Cl(Ca)(jct) were able also to shorten APD₆₀and provide negative feedback on cell Ca²⁺ load. Negligible effects due to 5% increases in 𝐼_Na, 𝐼_Na,b, 𝐼_NaK, 𝐼_Ks, 𝐼_Ca,b, and 𝐼_pCa currents’ fractions in the junctional cleft were found. Finally, Nimrod/E predicted slight or no effect on model outputs when 𝐹𝑥_CaL(jct) and 𝐹𝑥_NCX(jct), or any other two-parameter 𝐹𝑥_j(jct) groups, were combined. These model studies imply that more accurate experimental knowledge of 𝐼_Ca, 𝐼_Cl(Ca) and 𝐼_NCX, fluxes distributions is needed to better understand the physiological role of 𝐼_Cl(Ca) current, the local Ca²⁺ signaling, and the electrical cell activity. Furthermore, the estimates of % total flux in junctional cleft versus submembrane space under control conditions are still under debate. Additional quantification of the model sensitivity to ±10% changes in 𝐹𝑥_CaL(jct) or ±10% and ±15% variations in other transporters’ densities in the cleft (not done in this study) may provide a more convincing conclusions.

4.3. Effects of Small Variations in Intracellular Ion Transport on AP Duration and Ca²⁺ Transients

In the Shannon et al. model the diffusion of Ca²⁺ and Na⁺ ions was considered where (1) the junctional compartment communicates only with the SL compartment, (2) the SL compartment with the junctional and the bulk cytosolic compartments, and (3) the bulk cytosolic compartment only with the SL compartment. The values of the effective diffusion coefficients for Ca²⁺ and Na²⁺ were average values measured in different cardiac tissues and species under control conditions. In the literature, however, these estimates are contradictable [1, 8, 10–20]. In this study, we tested how small changes (±5, ±10%) in suggested diffusion coefficients may influence the model outputs. We used PLS and Nimrod toolkits to determine the most influential ion diffusion parameter(s) governing the properties of the electrical signal (such as APD₆₀) and of the Ca²⁺ signal (such as peak level of Ca²⁺ transient in the junctional cleft, subsarcolemmal space, and bulk cytosol). While 10% increases in basal Ca²⁺ diffusion coefficients (𝐷_Ca(jct-SL) and 𝐷_Ca(SL-𝑖)) had pronounced effects on the model outputs, no visible changes in all selected biomarkers were detected due to 10% increases in Na⁺ diffusion coefficients (𝐷_Na(jct-SL) and 𝐷_Na(SL-𝑖)). Interestingly marked effects on both APD₆₀ and Ca²⁺ peaks were predicted when Na⁺ diffusion was blocked, assuming 𝐷_Na(jct-SL)and 𝐷_Na(SL-𝑖) zeros shorten APD₆₀ and decreased [Ca]_𝑖, [Ca]_SL, and [Ca]_jct peaks (data not shown). Our results also demonstrate that accelerating Ca²⁺ transfer from junctional to submembrane space (e.g., increasing 𝐷_Ca(jct-SL)), (1) sensitively decreased Ca²⁺ peak in the junctional cleft and prolonged [Ca]_jct time to peak, (2) enhanced [Ca]_SL and [Ca]_𝑖 peaks, and (3) prolonged APD₆₀. The effects of 𝐷_Ca(SL-𝑖) increase were quite different. Accelerating Ca²⁺ transfer from the SL compartment to the bulk cytosol (e.g., increasing 𝐷_Ca(SL-𝑖)) shortened APD₆₀ and enhanced all Ca²⁺ transients’ peaks. In addition, 𝐷_Ca(SL-𝑖) increase most significantly enhanced [Ca]_jct peak while the control [Ca]_SL peak was less affected. Possible reason for the predicted increase in [Ca]_jct peak, due to 𝐷_Ca(SL-𝑖) increase, was enhanced [Ca]_𝑖 which promoted increase in the SR load (i.e., [Ca]_SR) in this way causing additional Ca²⁺ release from the SR and consequent elevations in [Ca]_jct and [Ca]_SL, respectively. The enhanced [Ca]_SL transient produced faster Ca²⁺-dependent inactivation of 𝐼_Ca, thus leading to APD shortening. Nimrod/E showed also that simultaneous increases in control 𝐷_Ca(jct-SL) and 𝐷_Ca(SL-𝑖) slightly increased Δ[Ca]_𝑖 while all other biomarkers remained insensitive to the changes in any other two-level parameter combinations. Also predicted relative effects of ±5% changes in intracellular ion transport were identical to those of ±10% variations (data not shown). In summary, our results reveal that small variations in diffusion coefficients for Ca²⁺ may impact sensitively the cell electrical activity and Ca²⁺ homeostasis. New and more precise estimates of the effective Ca²⁺ and Na⁺ diffusion constants are needed to better understand ECC in rabbit ventricular myocytes under control and certain pathologies.

We need to acknowledge that although we demonstrated good correlation with experiment, the results of sensitivity analysis strongly depend on the model details and simplifications. For example, in the lumped Shannon et al. model the ionic diffusion from junctional cleft to cytosolic compartment goes solely through subsarcolemmal space which is questionable with regard to the contact areas between these compartments [19]. In addition, in all studied cases the variations of intracellular potassium concentration were disabled (i.e., 𝑑[K]_𝑖/𝑑𝑡=0 in Shannon’s MatLab code) and whether or not this drawback was eliminated in the Romero’s code or in recent modifications of the Shannon’s code is unclear [4–7, 9]. This limitation leaded to the following: (1) fast stabilisation of the model for [Ca]_𝑖 and AP (after 6–8 s until 8 min); (2) [Na]_𝑖 appeared stable after ~200 s; (3) diastolic [K]_𝑖 remained unchanged in the interval 0≤𝑡≤8 min. Thus, models neglecting slow intracellular concentration changes in [K]_𝑖 produced mainly by Na/K pump need to be revised further. Finally, the 𝑡-tubule compartment, the ion diffusion between the 𝑡-tubule lumen and bulk extracellular space, and the mitochondrial Ca²⁺ fluxes were omitted [46–50]. These are likely to have important impact on predicted AP and Ca²⁺ profiles as well.

5. Conclusions

Our studies demonstrate that parameter sensitivity approaches can be used to obtain new insights into relationships between model parameters, model outputs, and experimental data. More specifically, the simulation results provided information on how small variations in the ionic currents and intracellular Ca²⁺ and Na⁺ diffusion modulate ventricular myocyte electrophysiology and Ca²⁺ dynamics in rabbit ventricular myocytes at 1 Hz. We found that APD₆₀ and Ca²⁺ peaks are (1) highly sensitive to 10% increase in 𝐼_Ca, while the effects of 𝐼_NCX, 𝐼_NaK, 𝐼_Kr, 𝐼_to, and 𝐼_ClBk were moderate, and (2) insensitive to 10% increases in all other Ca²⁺, Na⁺, and K⁺ transmembrane currents and SR Ca²⁺ fluxes. The myocyte electrical activity is highly sensitive to 5% variations in 𝐼_Ca and 𝐼_NCX and less sensitive to 𝐼_Cl(Ca) redistribution between junctional cleft and submembrane space. Small changes in submembrane and cytosolic diffusion coefficients for Ca²⁺, but not in Na⁺ diffusion coefficients, may impact notably the myocyte contraction. This study highlights the need for additional and more precise experimental data and further updating and testing of the Shannon et al. model and its recent modifications to fill specific knowledge gaps related to rabbit ventricular electrophysiology and intracellular Ca²⁺ signaling.

Appendix

See Tables 1, 2, and 3.

Glossary

Abbreviations

ECC:	Excitation-contraction coupling
AP:	Action potential
APD60:	Action potential duration at 60% repolarization
SL:	Sarcolemmal membrane
SR:	Sarcoplasmic reticulum
RyR:	Ryanodine receptor
SERCA2:	SR Ca2+ ATPase pump
TnC:	Troponin C
NIMROD:	Parameter modeling toolkit for sensitivity analysis
PLS:	Partial least square tool.

Membrane Currents

𝐼Na:	Na+ current
𝐼Na,b (𝐼NaBk):	Na+ background current
𝐼Ca (𝐼CaL):	L-type Ca2+ channel current
𝐼Ca,b (𝐼CaBk):	Ca2+ background current
𝐼to,f:	Fast transient outward K+ current
𝐼to (𝐼to,s):	Slow transient outward K+current
𝐼Kr:	Rapid delayed rectifier K+ current
𝐼Ks:	Slow delayed rectifier K+ current
𝐼K1:	Inward rectifier K+ current
𝐼Cl(Ca):	Ca2+ activated Cl− current
𝐼ClBk:	Cl− background current
𝐼NaK:	Na+-K+ pump current
𝐼NCX:	Na+-Ca2+ exchanger current
𝐼pCa (𝐼SLCaP):	Sarcolemmal Ca2+ pump current.

Concentrations

[Ca]𝑖:	Cytosolic Ca2+ concentration
[Ca]SL:	Subsarcolemmal Ca2+ concentration
[Ca]jct:	Junctional cleft Ca2+ concentration
[Ca]SR:	SR Ca2+ concentration
[Na]𝑖:	Cytosolic Na+ concentration
[K]𝑖:	Cytosolic K+ concentration.

Acknowledgments

The authors thank Dr. Jiří Šimurda for helpful suggestions and careful reading of the paper. This work is supported by NIBIB 1R03EB014593-01 (Michailova), National Biomedical Computational Resource (NIH grant 8P41 GM103426-19), NHLBI Systems Biology Collaboration 1R01HL105242-01 (McCulloch), and NSF funded PRIME program (NSF grant IOSE 0710726).