Little is known about how small variations in ionic currents and and diffusion coefficients impact action potential and 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 peaks are highly sensitive to 10% increase in L-type current; moderately influenced by 10% increase in - exchanger, - pump, rapid delayed and slow transient outward currents, and background current; insensitive to 10% increases in all other ionic currents and sarcoplasmic reticulum fluxes. Cell electrical activity is strongly affected by 5% shift of L-type channels and - exchanger in between junctional and submembrane spaces while -activated -channel redistribution has the modest effect. Small changes in submembrane and cytosolic diffusion coefficients for , but not in 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 signaling.

1. Introduction

Several ionic models have been developed to investigate the subcellular mechanisms regulating excitation-contraction coupling (ECC) in rabbit ventricular cardiomyocytes [17]. In 2004, Shannon and colleagues published a detailed model for handling and ionic current that accurately represents sarcoplasmic reticulum (SR) -dependent release and simulates basic ECC phenomena. This model was the first to include (1) the subsarcolemmal compartment to the other two commonly formulated cytosolic compartments (junctional and bulk, [8]), (2) the variations in the locations of ion transporters throughout the cell surface membrane, (3) and transport between the subcellular compartments, and (4) buffering inside. Latest studies extended further the Shannon et al. ionic model in rabbit ventricular cells. Mahajan and colleagues modified L-type current and 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 and examined the relative contributions of voltage- and -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 signal regulation in two rabbit cell models [1, 4]. Systematically characterizing AP properties, and 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 levels are strongly dependent on , , and 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 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 [1019]. Recently, the contributions of structural heterogeneities to local and global 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 diffusion and buffering on ECC were not investigated.

Taken together the published computational models have provided important insights into the underlying mechanisms regulating intracellular 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 and transport between cellular subdomains affect important cellular biomarkers (i.e., AP shape, , , and ). 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 and 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 fluxes included in the Shannon et al. electrophysiological model is shown in Figure 1.

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 concentrations ( ) 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, 2628].

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 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 signals ( , , ) at 9-10 s of every simulation. In all numerical experiments (  min) remained unchanged while the variations in were small (ranging from 20 μM to 100 μM) with insignificant effects on calculated AP and traces at 6–8 s. The and delta ( , , , ) quantities were used to analyze the model sensitivity. is defined as the difference between peak concentration in the particular compartment and the diastolic 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 13). 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, . 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 and the output is 0.5, then when is one standard deviation greater than the mean, will increase by half a standard deviation. Conversely, if the value is −0.5, then when is one standard deviation greater than the mean, 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.

3. Results

3.1. Effects of 10% Changes in Channels Conductances on AP Morphology and Intracellular Ca2+ 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 (APD60, , , ).

The increase in L-type channel permeability ( ) had the most pronounced effects on and signals, by increasing action potential duration and enhancing maximum peaks. The effects of , , and , while less pronounced than , were still significant. Figure 2 also shows that 10% increases in and decreased and peaks while 10% increase in decreased peaks but increased . The increase in maximal pump rate ( ) had still detectable but minor effect, by decreasing both and peaks. The changes in SR parameters ( , , and ) by 10% had negligible effects on and whereas the submembrane and junctional peaks (e.g., and ) were most affected. Figure 2 shows that 10% increases in and increased and while had the opposite effect. Interestingly, 10% increase in background conductance ( ) showed similar magnitude as and and demonstrated the same effects in decreasing and 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 and KNCX combined, as well all other two-parameter group combinations (not shown in Figure 2), had slight or no effect on the selected cellular biomarkers ( , , , and ). The plots in Figure 3 confirm further the specific, quantitative predictions made by PLS and Nimrod/E for the effects of 10% increases in , , , , , and on AP morphology and transients.

3.2. Effects of Changes in Membrane Transporter Distributions on AP Morphology and Intracellular Ca2+ 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 channels were assumed concentrated within the junctional membrane such that 90% of total number were located there while all other membrane transporters ( channels, leak, / pump, slow delayed rectified channel, activated channel, leak, / exchanger, and sarcolemmal 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 ( , , , and ), we performed systematic analysis, increasing by 5% the basic while decreasing to keep total number of transporter protein complexes within the sarcolemma unchanged; that is, (the appendix, Table 2). Figure 4 shows that PLS and Nimrod/E analysis again showed consistent results. The 5% increase in L-type channels fraction in the junctional cleft ( ) had the most pronounced effect on the model outputs, decreasing and peaks. The results also suggest that 5% increase in the junctional - current ( ) prolonged while slightly affected peaks in the junctional cleft and subsarcolemmal and bulk cytosol compartments. The 5% increase in the junctional ( ) had similar effect on all model outputs as the increase in but with much smaller magnitude. Results in Lenth plot (Figure 4) for the four model outputs reveal that “ 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 by 5% and decreasing ( = 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.

3.3. Effects of 10% Changes in Diffusion Constants on AP Morphology and Intracellular Ca2+ Dynamics

The effects of changes in default parameters describing and 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 altered all model outputs, prolonging and increasing and but decreasing . Moreover 10% increases in basal value shortened and increased , , and . Interestingly, 10% increases in basal and values in all three compartments had no effect on and peaks. Nimrod/E analysis also suggests that two-level interaction between and had a minor effect on increasing while , , and were unaffected. Figure 7 additionally confirms the PLS and Nimrod/E predictions for the effects of 10% increases in and on AP morphology and transients.

4. Discussion

The main goal of this study was to perform a systematic investigation on how small variations in ionic currents and and diffusion coefficients modulate ventricular myocyte electrophysiology and 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 ( , , , and ) 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 Ca2+ Transients

Using parameter sensitive analysis, Romero et al. quantified the sensitivity of AP properties, and 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 current properties. Results indicate that 10% increase in maximal permeability of L-type channel ( ), with respect to other currents’ conductance changes, most notably affected and transients’ peaks at 1 Hz. Our predictions for the effects of on APD and cytosolic peak are also consistent with reported data in rabbit ventricular myocytes after applications of blockers [5, 9, 29] or in cells with high L-type channel density [30]. The effects of 10% increase in , , , or on cellular biomarkers, while less pronounced than , were still detectible. Our predictions for increase at 1 Hz are in qualitative agreement with data in ferret myocytes [31] and data supporting APD prolongation following block reported in some experiments in rabbit cardiomyocytes [32, 33]. Several studies in rabbit ventricular cells report also a strong effect of in modulating systolic transient [3134] while others report only a moderate effect [35], similar to our predictions for . Results demonstrate also that 10% increase in slightly decreased and peaks. Increase in and at 1 Hz resulted in APD shortening, which is consistent with experiments in rabbit ventricular myocytes [36, 37]. Important role of in modulating transients through its effects on and during AP repolarization period has been suggested [38, 39]. Slight effects, however, due to 10% increase in at 1 Hz on peaks are predicted here. Furthermore, at stimulus frequency of 2.5 Hz, Romero et al. report (±30% changes or severe block; 0–8 min) as important influent current in most simulated cellular biomarkers, including [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 (increment 5%, 0–8 min) on was less pronounced. The main reason for this observation was that the predicted changes in were smaller at 1 Hz (20–100 μM) versus 1–1.5 mM at 2.5 Hz, in this way affecting differently and 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 background current ( ), which is a new feature of the Shannon et al. model, showed a negative relationships with four outputs ( , , , and ) of approximately similar magnitude as for or . Experimental data reporting the effects of on rabbit ECC are missing in the literature however. Another interesting finding is the relatively small sensitivity of the biomarkers to 10% increases in , , , , , , , , , , , , and 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) was highly sensitive to 100% , , , and and block and less sensitive to , , , and inhibition while 100% block of , , , and currents and SR leak ( ) had small or no effect; (2) steady-state peaks were strongly affected when , , and SR fluxes ( , , and ) were inhibited whereas 100% block , , , , and had less pronounced effects. New measurements in rabbit ventricular myocytes are needed to define the relative importance of , , , , , , , , , , , , and 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 dynamics. Our analysis suggests that 10% increase in both and had minor effect at 1 Hz. Additionally, Nimrod/E showed that simultaneous increases in control and 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 inhibitors suppress partially -related currents also; therefore the balance of these currents is critical for predicting -related drug action.

4.2. Effects of Variations in Membrane Transporters Distribution on AP Duration and Ca2+ 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, 4042]. Using reaction-diffusion models, recently we demonstrated that the subcellular signals are highly sensitive to and fluxes distributions via the sarcolemma [13, 16, 19, 20]. In these spatial models, however, the cell electrophysiology was simplified; equations describing , , , and properties were only included; and AP dynamics were not reproduced. Additionally, small variations in L-type channels density in junctional cleft have been suggested to lead to large changes in control SR release, subcellular 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 transients. We evaluated the biomarkers’ sensitivity to 5% increases in the fraction of total ionic current in the junctional cleft ( ) assuming total number of transporters within the sarcolemma unchanged; that is, . Our results yielded three channel distributions ( , , and ) that significantly affect the model outcomes of interest (i.e., , , , and ). Increasing by 5% the density of L-type channels in the junctional cleft (while decreasing by 45% density in SL compartment to keep ) had the strongest effect on both and peaks, decreasing and , , and peaks. The reason for the predicted decrease in was the strongly increased rate of transfer from the junctional cleft to the submembrane space evoked by the 5% increase of and reduced (~45%) into the SL compartment. The predicted decreases into peak and were also consequence of reduction to 45% though flux from the cleft was increased. The drop in evoked decrease in the rate of transfer from the SL to cytosolic compartment thus causing drop in peak as well. Our studies also showed that the effects of redistribution were quite different. The 5% increase in (while decreasing by 40-41% density in SL compartment to keep ) prolonged to approximately the same magnitude as for and slightly decreased the control peaks. Significant effects due to redistribution and negligible effects due to redistribution on the local and global transients’ peaks were also found using our reaction-diffusion model in rabbit ventricular cells [20]. Calcium activated chloride current ( ) has been extensively studied [4345]. Experimental data in rabbit ventricular myocytes suggest that is strongly temperature dependent (very small at room temperature but substantial at 35°C–37°C); can be eliminated by blocking ; and may normally facilitate rapid repolarization when SR load and release are high. To bring further insights into the effects of properties on ECC, we examined how 5% increase in density in the junctional cleft (assuming ) may affect the cell electrical activity and signaling. Interestingly, although much less pronounced than variations, changes in were able also to shorten and provide negative feedback on cell load. Negligible effects due to 5% increases in , , , , , and currents’ fractions in the junctional cleft were found. Finally, Nimrod/E predicted slight or no effect on model outputs when and , or any other two-parameter groups, were combined. These model studies imply that more accurate experimental knowledge of , and , fluxes distributions is needed to better understand the physiological role of current, the local 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 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 Ca2+ Transients

In the Shannon et al. model the diffusion of and 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 and were average values measured in different cardiac tissues and species under control conditions. In the literature, however, these estimates are contradictable [1, 8, 1020]. 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 ) and of the signal (such as peak level of transient in the junctional cleft, subsarcolemmal space, and bulk cytosol). While 10% increases in basal diffusion coefficients ( and ) had pronounced effects on the model outputs, no visible changes in all selected biomarkers were detected due to 10% increases in diffusion coefficients ( and ). Interestingly marked effects on both and peaks were predicted when diffusion was blocked, assuming and zeros shorten and decreased , , and peaks (data not shown). Our results also demonstrate that accelerating transfer from junctional to submembrane space (e.g., increasing ), (1) sensitively decreased peak in the junctional cleft and prolonged time to peak, (2) enhanced and peaks, and (3) prolonged . The effects of increase were quite different. Accelerating transfer from the SL compartment to the bulk cytosol (e.g., increasing ) shortened and enhanced all transients’ peaks. In addition, increase most significantly enhanced peak while the control peak was less affected. Possible reason for the predicted increase in peak, due to increase, was enhanced which promoted increase in the SR load (i.e., ) in this way causing additional release from the SR and consequent elevations in and , respectively. The enhanced transient produced faster -dependent inactivation of , thus leading to APD shortening. Nimrod/E showed also that simultaneous increases in control and slightly increased 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 may impact sensitively the cell electrical activity and homeostasis. New and more precise estimates of the effective and 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., 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 [47, 9]. This limitation leaded to the following: (1) fast stabilisation of the model for and AP (after 6–8 s until 8 min); (2) appeared stable after ~200 s; (3) diastolic remained unchanged in the interval  min. Thus, models neglecting slow intracellular concentration changes in 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 fluxes were omitted [4650]. These are likely to have important impact on predicted AP and 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 and diffusion modulate ventricular myocyte electrophysiology and dynamics in rabbit ventricular myocytes at 1 Hz. We found that and peaks are (1) highly sensitive to 10% increase in , while the effects of , , , , and were moderate, and (2) insensitive to 10% increases in all other , , and transmembrane currents and SR fluxes. The myocyte electrical activity is highly sensitive to 5% variations in and and less sensitive to redistribution between junctional cleft and submembrane space. Small changes in submembrane and cytosolic diffusion coefficients for , but not in 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 signaling.


See Tables 1, 2, and 3.


ECC: Excitation-contraction coupling
AP: Action potential
: Action potential duration at 60% repolarization
SL: Sarcolemmal membrane
SR: Sarcoplasmic reticulum
RyR: Ryanodine receptor
SERCA2: SR ATPase pump
TnC: Troponin C
NIMROD: Parameter modeling toolkit for sensitivity analysis
PLS: Partial least square tool.
Membrane Currents
: current
( ): background current
( ):L-type channel current
( ): background current
:Fast transient outward current
( ):Slow transient outward current
:Rapid delayed rectifier current
:Slow delayed rectifier current
:Inward rectifier current
: activated current
: background current
: - pump current
: - exchanger current
( ):Sarcolemmal pump current.
: Cytosolic concentration
: Subsarcolemmal concentration
: Junctional cleft concentration
: SR concentration
: Cytosolic concentration
: Cytosolic concentration.


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).