Sensitivity of Rabbit Ventricular Action Potential and Ca2+ Dynamics to Small Variations in Membrane Currents and Ion Diffusion Coefficients

Little is known about how small variations in ionic currents and Ca2+ and Na+ diffusion coefficients impact action potential and Ca2+ 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 Ca2+ peaks are highly sensitive to 10% increase in L-type Ca2+ current; moderately influenced by 10% increase in Na+-Ca2+ 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 Ca2+ fluxes. Cell electrical activity is strongly affected by 5% shift of L-type Ca2+ channels and Na+-Ca2+ exchanger in between junctional and submembrane spaces while Ca2+-activated Cl−-channel redistribution has the modest effect. Small changes in submembrane and cytosolic diffusion coefficients for Ca2+, 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 Ca2+ signaling.


Introduction
Several ionic models have been developed to investigate the subcellular mechanisms regulating excitation-contraction coupling (ECC) in rabbit ventricular cardiomyocytes [1][2][3][4][5][6][7]. In 2004, Shannon and colleagues published a detailed model for Ca 2+ handling and ionic current that accurately represents sarcoplasmic reticulum (SR) Ca 2+ -dependent release and simulates basic ECC phenomena. This model was the first to include (1) the subsarcolemmal Ca 2+ compartment to the other two commonly formulated cytosolic Ca 2+ compartments (junctional and bulk, [8]), (2) the variations in the locations of ion transporters throughout the cell surface membrane, (3) Ca 2+ 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 Ltype Ca 2+ current and Ca 2+ 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 2+ -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 2+ signal regulation in two rabbit cell models [1,4]. Systematically characterizing AP properties, Ca 2+ 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 steadystate Ca 2+ 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 2+ 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][11][12][13][14][15][16][17][18][19]. Recently, the contributions of structural heterogeneities to local and global Ca 2+ 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 2+ 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 2+ 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 2+ 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 2+ 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.

Electrophysiological Model in Rabbit Ventricular Myocytes.
A diagram describing the arrangement of subcellular compartments, ionic currents and pumps, and intracellular Ca 2+ 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 Ca 2+ concentrations ([Ca] ) changes under control conditions [1].

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.

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.

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.  Notes: fraction of total flux in SL ( (SL) = 1 − (jct) ). In the Shannon et al. article [1] Cl(jct) corresponds to Cl(Ca)(jct) . The baseline total " " flux fraction values used in this study are the same as in [1].  [1] Ca(junction-SL) corresponds to Ca(jct-SL) , Ca(SL-cytosol) to Ca(SL-) , Na(junction-SL) to Na(jct-SL) , and Na(SL-cytosol) to Na(SL-) . The baseline diffusion coefficients values used in this study are the same as in [1].
Nimrod/O performs multiple searches for an algorithm and multiple algorithms, all in parallel.

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.  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 twolevel 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 BioMed Research International 5 instance, if the regression coefficient for the input CaL and the output APD 60 is 0.5, then when CaL is one standard deviation greater than the mean, APD 60 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 60 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 2 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.

Effects of 10% Changes in Channels Conductances on AP Morphology and Intracellular Ca 2+
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 The increase in L-type Ca 2+ channel permeability ( CaL ) had the most pronounced effects on APD 60 and Ca 2+ signals, by increasing action potential duration and enhancing maximum Ca 2+ 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 60 and Ca 2+ peaks while 10% increase in NCX decreased Ca 2+ peaks but increased APD 60 . The increase in maximal NaK pump rate ( NaK ) had still detectable but minor effect, by decreasing both APD 60 and Ca 2+ peaks. The changes in SR parameters ( rel , up , and leak ) by 10% had negligible effects on APD 60 and Δ[Ca] whereas the submembrane and junctional Ca 2+ 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 60 and Ca 2+ 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 60 , Δ[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 2+ transients.

Effects of Changes in Membrane Transporter
Distributions on AP Morphology and Intracellular Ca 2+ 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 2+ 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 2+ activated Cl − channel, Ca 2+ leak, Na + /Ca 2+ exchanger, and sarcolemmal Ca 2+ 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 60 , Δ[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 2+ channels fraction in the junctional cleft ( CaL(jct) = 0.945) had the most pronounced effect on the model outputs, decreasing APD 60 and Ca 2+ peaks. The results also suggest that 5% increase in the junctional Na + -Ca 2+ current ( NCX(jct) = 0.1155) prolonged APD 60 while slightly affected Ca 2+ 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.

Effects of 10% Changes in Diffusion Constants on AP
Morphology and Intracellular Ca 2+ Dynamics. The effects of changes in default parameters describing Ca 2+ 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 Na(jct-SL) and Na(SL-) values in all three compartments had no effect on APD 60 and Ca 2+ 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 60 , Δ[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 2+ transients.

Discussion
The main goal of this study was to perform a systematic investigation on how small variations in ionic currents and Ca 2+ and Na + diffusion coefficients modulate ventricular myocyte electrophysiology and Ca 2+ dynamics in rabbit ventricular myocytes. To determine the most influential  most notably affected APD 60 and Ca 2+ transients' peaks at 1 Hz. Our predictions for the effects of CaL on APD and cytosolic Ca 2+ 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 2+ 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   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   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 60 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 2+ leak ( leak ) had small or no effect; (2) steady-state Ca 2+ peaks were strongly affected when Na , Ca,b , and SR Ca 2+ 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 2+ 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.

Effects of Variations in Membrane Transporters Distribution on AP Duration and Ca 2+
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][41][42]. Using reaction-diffusion models, recently we demonstrated that the subcellular Ca 2+ 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 2+ channels density in junctional cleft have been suggested to lead to large changes in control SR Ca 2+ release, subcellular Ca 2+ 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 2+ 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 60 , Δ 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 60 to approximately the same magnitude as for Ca and slightly decreased the control Ca 2+ peaks. Significant effects due to Ca redistribution and negligible effects due to NCX redistribution on the local and global Ca 2+ 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][44][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 2+ 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 2+ signaling. Interestingly, although much less pronounced than CaL(jct) variations, changes in Cl(Ca)(jct) were able also to shorten APD 60 and provide negative feedback on cell Ca 2+ 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 2+ 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. The values of the effective diffusion coefficients for Ca 2+ and Na 2+ were average values measured in different cardiac tissues and species under control conditions. In the literature, however, these estimates are contradictable [1,8,[10][11][12][13][14][15][16][17][18][19][20].

Effects of Small Variations in Intracellular Ion Transport on AP Duration and Ca
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 60 ) and of the Ca 2+ signal (such as peak level of Ca 2+ transient in the junctional cleft, subsarcolemmal space, and bulk cytosol). While 10% increases in basal Ca 2+ 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 60 and Ca 2+ peaks were predicted when Na + diffusion was blocked, assuming Na(jct-SL) and Na(SL-) zeros shorten APD 60 and decreased [Ca] , [Ca] SL , and [Ca] jct peaks (data not shown). Our results also demonstrate that accelerating Ca 2+ transfer from junctional to submembrane space (e.g., increasing Ca(jct-SL) ), (1) sensitively decreased Ca 2+ peak in the junctional cleft and prolonged [Ca] jct time to peak, 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 2+ may impact sensitively the cell electrical activity and Ca 2+ homeostasis. New and more precise estimates of the effective Ca 2+ 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][5][6][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 2+ fluxes were omitted [46][47][48][49][50]. These are likely to have important impact on predicted AP and Ca 2+ profiles as well.

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 2+ and Na + diffusion modulate ventricular myocyte electrophysiology and Ca 2+ dynamics in rabbit ventricular myocytes at 1 Hz. We found that APD 60 and Ca 2+ 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 2+ , Na + , and K + transmembrane currents and SR Ca 2+ 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 2+ , 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 2+ signaling.