Effect of Ion Concentration Changes in the Limited Extracellular Spaces on Sarcolemmal Ion Transport and Ca2+ Turnover in a Model of Human Ventricular Cardiomyocyte

We have developed a computer model of human cardiac ventricular myocyte (CVM), including t-tubular and cleft spaces with the aim of evaluating the impact of accumulation-depletion of ions in restricted extracellular spaces on transmembrane ion transport and ionic homeostasis in human CVM. The model was based on available data from human CVMs. Under steady state, the effect of ion concentration changes in extracellular spaces on [Ca2+]i-transient was explored as a function of critical fractions of ion transporters in t-tubular membrane (not documented for human CVM). Depletion of Ca2+ and accumulation of K+ occurring in extracellular spaces slightly affected the transmembrane Ca2+ flux, but not the action potential duration (APD90). The [Ca2+]i-transient was reduced (by 2%–9%), depending on the stimulation frequency, the rate of ion exchange between t-tubules and clefts and fractions of ion-transfer proteins in the t-tubular membrane. Under non-steady state, the responses of the model to changes of stimulation frequency were analyzed. A sudden increase of frequency (1–2.5 Hz) caused a temporal decrease of [Ca2+] in both extracellular spaces, a reduction of [Ca2+]i-transient (by 15%) and APD90 (by 13 ms). The results reveal different effects of activity-related ion concentration changes in human cardiac t-tubules (steady-state effects) and intercellular clefts (transient effects) in the modulation of membrane ion transport and Ca2+ turnover.


Introduction
The increased availability of experimental data obtained from human isolated cardiomyocytes in recent years provided sufficient material to mathematically describe the transmembrane ion transport and electrical activity in human atrial [1,2] and ventricular cells [3][4][5][6][7][8]. None of these models, however, included a description of the effects of simultaneous changes of ion concentrations in the intercellular cleft space and t-tubular system as functional compartments that may affect cellular activity. The present work is the first attempting at filling this blank.
The first model of cardiac cells that included changes of Ca 2+ concentration in the intercellular clefts was developed by Hilgemann and Noble [9]. Focused on the description of Ca 2+ movements and excitation-contraction (EC) coupling in rabbit atrial cells, the model accurately reproduced the variations of extracellular Ca 2+ concentration measured in rabbit atrium [10,11]. Later, intercellular cleft space was also introduced into a model of human atrial cell by Nygren et al. [2].
In 2003, we published a pioneering work [12] indicating that activity-induced ion concentration changes in t-tubules of cardiac ventricular cells may be sufficiently high to modulate membrane currents and cellular electrical activity. Later, works on species-specific models for rat [13] and guinea-pig [14] demonstrated that, during a single beat at steady-state stimulation, the transient Ca 2+ depletion in t-tubules may range between 7% and 20%, while transient K + accumulation may reach 4%. This depended on the characteristics of the ion transporters, on their distribution between the surface and the t-tubular membrane, on the rate of ion diffusion between the t-tubular lumen and the extracellular space and on stimulation frequency. The resulting effect was a decrease of transmembrane Ca 2+ influx and a reduction of sarcoplasmic reticulum (SR) Ca 2+ load and intracellular Ca 2+ -transient. The frequency-dependent profiles of ion concentration transients during a cycle were different between these two models, stressing the need for a dedicated human ventricular cell model.
The present model of human ventricular myocyte has been formulated to take into account the up-to-date knowledge of membrane ionic currents and intracellular homeostasis in human cells. The description of the t-tubular system is based on current knowledge of morphological features in human cardiomyocytes [15]. To increase the accuracy of the model, the concentration changes in the layer adjacent to the tubular membrane were introduced similarly as in our latest work [16]. The t-tubular space was partitioned into concentric subspaces to describe the radial concentration gradients within the t-tubular lumen. In addition, to be consistent with the presence of the limited size of the pericellular interstitial space, a cleft space has been formulated, as done by Nygren et al. [2].
The aim of this study was to explore the functional significance of the presence of these extracellular spaces. Because the values of some parameters that are related to the features of t-tubules in human ventricular myocytes are still unknown (fractions of ionic currents in the t-tubular membrane, the rate of ion exchange between the t-tubular and extracellular space), we explored the physiological consequences of their various combinations, in particular those related to cellular electro-mechanical activity. Figure 1 illustrates the APs, three dominant ionic currents involved in membrane repolarization (I CaL , I Kto and I K1 ) and currents controlling intracellular ionic homeostasis (I NaCa , I NaK ) at 1 and 2.5 Hz steady-state stimulation. While the differences between APs at surface and t-tubular membranes are minimal [17], the differences between surface and tubular components of ionic currents are marked. They are caused by uneven distribution of corresponding ion transporters between both membrane systems (see Table S1 in the Supplementary Information), but also partly by differences in extracellular ion concentrations, to which both membranes are exposed during the stimulation cycle ( Figure 2).  The depletion of t-tubular Ca 2+ is predominantly induced by the activation of tubular I CaL (I CaL,t ). The two peaks in the tubular K + accumulation result from the activation of tubular I Kto (I Kto,t ) at the beginning of the AP and from the increase of tubular I K1 (I K1,t ) during the course of AP repolarization. The traces of tubular [Ca 2+ ] and [K + ] in Figure 2 show how activity-induced depletion of tubular Ca 2+ and accumulation of tubular K + decrease from the layer adjacent to the t-tubular membrane (first segment; [Ca 2+ ] t1 , [K + ] t1 ) to the central part of the t-tubule lumen. Their differences demonstrate luminal ion concentration gradients in the radial direction. The alterations of [Ca 2+ ] in clefts ([Ca 2+ ] c ) were induced by the activation of the main Ca 2+ currents in the surface membrane (I CaL,s , I pCa,s , I NaCa,s ) and by Ca 2+ diffusion between t-tubules and clefts. Small changes of [K + ] c were induced by activation of I Kto,s , I K1,s and I NaK,s and by K + diffusion between t-tubules and clefts. In addition, concentrations of all ions in the cleft space were affected by diffusion of ions between the cleft space and the external bulk space.  2). The model was run when ion concentration changes in extracellular spaces were either allowed to change (red lines) or fixed at external bulk levels (black lines).

Effect of Ion Concentration Changes in Extracellular Spaces on Intracellular Ca 2+ SR Load and Ca 2+ -Transient at Steady-State Stimulation
Though small, this effect shows a potential of these variations to slightly affect the Ca 2+ turnover in cardiac cells at steady state. The time integral of intracellular Ca 2+ -transients was reduced by 2.5% and 3.4% at 1 and 2.5 Hz, respectively. A deeper analysis also showed that this effect is predominantly induced by depletions of Ca 2+ concentration in the vicinity of the t-tubular membrane ([Ca 2+ ] t1 ); the relative contribution of changes of [Ca 2+ ] c , [K + ] c and [K + ] t1 in total was below 0.3%.
To explore the effect of Ca 2+ concentration changes in t-tubules on cytosolic Ca 2+ -transient for potentially different combinations of critical fractions f CaL,t , f NaCa,t and f pCa,t in human ventricular cells, we performed a series of simulations enabling us to construct these relations in 3D graphs ( Figure 4). The results displayed for t-tubular fraction of L-type Ca 2+ channels of 0.64 {based on data from guinea-pig [18] (basic setting)} showed that, if the fraction, f pCa,t , was small (<0.2), the repeated depletions of Ca 2+ in human t-tubules might result in a reduction of Ca 2+ -transient by 3% at resting heart rate (1 Hz) and by 4% during exercise (2.5 Hz) ( Figure 4, left panel). However, increasing the t-tubular fraction of L-type Ca 2+ channels to 0.8 (consistent with data from rat cardiomyocytes [19]) would lead to a reduction of Ca 2+ -transient by more than 4% and 5.4%, respectively, under the same conditions ( Figure 4, right panel).
If, in addition, the rate of Ca 2+ diffusion between human t-tubules and the cleft space was made slower (τ Ca,ct increased from the basic value of 240 to 480 ms), this effect would further increase to more than 9% at 2.5 Hz (not illustrated). While the dependency of this effect on the value of f NaCa,t appeared to be small, high values of f pCa,t would even reverse it, as they would lead to accumulation of Ca 2+ in the tubular space and, consequently, to a slight increase of the Ca 2+ -transient.  Figure 5 shows the responses of the model to a sudden change of stimulation frequency from 1 to 2.5 Hz followed during 10 s (left panels) and 600 s (right panels). It demonstrates the potency of ion concentration changes in both extracellular spaces to modulate NSR Ca 2+ load and intracellular [Ca 2+ ], [Na + ] and [K + ] during transient events after a sudden change of stimulation frequency. The results from the basic model (blue dots) are compared with the results after fixing ion concentrations in the cleft space (red dots) or in both t-tubular and cleft spaces (black dots) at external bulk levels. The right panels illustrate progression toward the steady state in an expanded time scale.

Model Responses to Sudden Changes of Stimulation Rate
The fast transient increase of Ca 2+ flux into the cell decaying within several seconds (integral of total trans-sarcolemmal Ca 2+ flux during a cycle; panel a, left) induced a marked transient reduction of mean Ca 2+ concentration in the cleft space (panel b, left) that was also reflected by a similar change in the mean Ca 2+ concentration in the first tubular segment (panel c, left). Note that this reduction of [Ca 2+ ] in both extracellular spaces was substantially larger (blue dots) than the reduction of [Ca 2+ ] t1 in the model after fixing ion concentrations in clefts space (red dots). This initial rapid reduction in extracellular [Ca 2+ ] decays with a time constant of approximately 50 s due to equilibration between the clefts and bulk space (panels b, right, and c, right).
The consequence of such a profound transient depletion of Ca 2+ in both cleft and t-tubular spaces was a parallel decrease of Ca 2+ concentration in NSR and a significant reduction of cytosolic These results suggest that, while the role of ion concentration changes in cleft spaces in modulating cellular electrical activity and the Ca 2+ -transient is relatively small at steady state, it becomes important during perturbations, due to a change in stimulation frequency. This is illustrated in Figure 6 showing APs, Ca 2+ flow through sarcolemma and [Ca 2+ ] i -transient at a stimulation frequency of 1 and 2.5 Hz under conditions of free running and fixed extracellular concentrations. Steady-state differences in the course of AP were negligible, and changes in [Ca 2+ ] i -transient were small at any frequency. In contrast, both quantities were markedly affected 10 s after a sudden change of stimulation rate from 1 to 2.5 Hz. APD 90 (the duration of the action potential at 90% repolarization) was reduced by 13 ms in accordance with a decreased influx of Ca 2+ ions in the course of AP. [Ca 2+ ] i -transient dropped by 15%. If the stimulation frequency returned from 2.5 to 1 Hz, the opposite deviations in AP and [Ca 2+ ] i -transient appeared (not shown).
To verify the ability of our model to reproduce the changes in [K + ] and [Ca 2+ ] in the clefts that were documented in the literature, we simulated, from steady-state stimulation at 1 Hz, the application of a train of stimuli at 3 Hz for durations from 34 up to 700 s and a subsequent return to 1 Hz for 700 s. Figure 7 shows the end-diastolic values of K + and Ca 2+ concentrations in the cleft space ([K + ] c,end and [Ca 2+ ] c,end ) and the end-diastolic values of resting transmembrane voltage (V m,end ) at the peripheral membrane. During the train of stimuli applied at 3 Hz, K + accumulation was taking place initially with a time constant of 10 s and peaked at 0.8 mM above the initial level of 5.41 mM at 1 Hz. While the stimulation train was continued, it thereafter decayed with a time constant of 214 s toward the initial level at 1 Hz. A similar time course was observed for the K + depletion upon returning to 1 Hz stimulation. In addition, the peak level of [K + ] c,end during the depletion depended on the duration of the stimulation train at high frequency. The left and right columns show the responses of the model during 10 and 600 s, respectively, after the change of stimulation rate. The model was run when t-tubular, and cleft ion concentrations were allowed to change (blue), when ion concentrations in the cleft space were fixed (red) and when both t-tubular and cleft ion concentrations were fixed at external bulk levels (black). Figure 6. The effect of ion concentration changes in extracellular spaces on the action potential (AP) (V m ), total sarcolemmal Ca 2+ flux (J Ca,tot ) and cytosolic Ca 2+ -transient after a sudden increase of stimulation frequency from 1 to 2.5 Hz. The left, middle and right columns show the results obtained, respectively, at 1 Hz steady-state stimulation, at 10 and at 600 s (steady state) after a sudden increase of simulation rate to 2.5 Hz. The model was run when t-tubular and cleft ion concentrations were allowed to change (blue lines), when ion concentrations in the cleft space were fixed (red lines) and when both t-tubular and cleft ion concentrations were fixed at external bulk levels (black, practically identical to red traces).
These features compare well with the measurements of interstitial [K + ] done with a K + -selective microelectrode by Kunze [20]. Another feature of Kunze's data that was well reproduced in our model is that the peak magnitude of the K + depletion upon returning to 1 Hz depended on the duration of the stimulation at 3 Hz. The insets show the corresponding experimental records digitized from Kunze [20] and Attwell et al. [21].
Changes in the membrane voltage closely paralleled the kinetics of [K + ] changes in the clefts: the membrane was depolarized during K + accumulation and hyperpolarized during K + depletion, by about 4 mV. The onset of the depolarization at the application of 3 Hz stimulation and that of the hyperpolarization upon returning to 1 Hz stimulation were as fast as from the microelectrode recordings of Attwell et al. [21]. After the peak, the time courses of the subsequent decay of the depolarization or of the hyperpolarization were comparable to those of the membrane voltage signal in [21] and of the [K + ] signal in [20]. Thus, the cumulative K + accumulation and membrane depolarization measured at the end of diastole in the model during a high frequency stimulation train are comparable to published data, given the differences with our modelling: (i) unlike Kunze's and Attwell et al.'s data, which were recorded in the atrial tissues of rabbit and guinea-pig, respectively, our model is designed for human ventricular cells; and (ii) the tissues were not perfused, whereas our model assumes an infinite volume of the bulk solution, which is equivalent to a perfect perfusion.
In our simulations, the reduction of [Ca 2+ ] c,end at the beginning of the 3 Hz stimulation train was fast and amounted to 0.25 mM (Figure 7). Subsequently, cleft [Ca 2+ ] decayed with a time constant of 46 s, while the model cell was still stimulated at 3 Hz. In addition, our simulations show that, upon returning to 1 Hz, a transient accumulation of Ca 2+ took place in the cleft space with a very similar time course as the former depletion. This compares well with Ca 2+ depletion measured in the perfused right ventricular tissue of rabbits (amounting to 0.3 to 0.5 mM) reported by Hilgemann and Langer [22] in a train of eight stimulations at 0.5-s intervals applied from rest.

Discussion
This study demonstrates that ion concentration changes in the t-system of human cardiomyocytes, though small compared to those observed in skeletal muscle [23,24], are not negligible and should be considered in a detailed description of a human ventricular myocyte as a factor contributing to the control of Ca 2+ turnover and cardiac contraction. This confirms our previous results from simulations in models of animal cardiac ventricular myocytes of guinea-pig [14] and rat [16]. In addition, [Ca 2+ ] changes occurring in the limited space of t-tubules are potentiated by [Ca 2+ ] changes in the intercellular clefts, primarily during transients caused by changes in stimulation frequency.

Validation of the Model
The ability of the present model to simulate changes of APs recorded in experiments in human cardiomyocytes at different stimulation frequencies [25] is demonstrated in Figure 8. This frequency dependence of APD was shown to be associated primarily with modulation of I CaL [26]. For the sake of comparison, we reconstructed the frequency dependence of AP configuration also by using three other models of human ventricular myocytes recently published by Iyer et al. [5], Fink et al. [6] and O'Hara et al. [8]. To strictly meet the experimental conditions specified in [25], we applied only 15 stimulation pulses, and [Na + ] i and [K + ] i were fixed at 10 and 130 mM, respectively (ion concentrations used in the pipette solution). Surprisingly, the frequency-dependent changes of AP duration in the models ( [6] and [8]) were quite small. Larger changes in AP duration could be observed in these two models only when a substantially higher number of stimulation pulses were applied (e.g., 300 pulses) and when the simulations were done without fixation of [Na + ] i and [K + ] i .
The ability of our model to simulate the frequency-dependent changes of intracellular ion concentrations measured in human preparations by Schmidt et al. [27] ([Ca 2+ ] i ) and Pieske et al. [28] ([Na + ] i ) is illustrated in Figure 9. The data are normalized to the values at 1 Hz. The resting value of [Ca 2+ ] i in our model (63.6 nM) agrees well with the experimental data of Beuckelmann et al. [29], and that of [Na + ] i (9.3 mM) is in the range of the experimentally observed values of Pieske et al. [30].  [25]. The following models were used: (a) the present model; (b) the model of Iyer et al. [5]; (c) the model of Fink et al. [6] and (d) the model of O'Hara et al. [8]; (e) Representative AP waveforms digitized from Li et al. [25]; and (f) Comparison of values of APD 50 and APD 90 (action potential duration at 50% and 90% repolarization, respectively) at 0.5, 1 and 2 Hz obtained from the experimental data of Li et al. [25] and the present model. In the models, 15 APs were elicited at each stimulation rate, as done in [25].  [27]; and (b) [Na 2+ ] i (600 APs elicited at each stimulation rate); data published by Pieske et al. [28].

Conditions of Model Stability
The model showed long-term stability, as demonstrated in Figure 11. For all types of ions, the conditions of stability required that, at steady state, the net amount of substance transferred across the cell membrane during one period (T) at any frequency equaled zero. Particularly for calcium ions, this condition can be expressed as: where n Ca,net,s and n Ca,net,t denote net amounts of Ca 2+ ions transferred during each stimulation period, T, across the surface and the t-tubular membrane, respectively. They can be expressed as the time integrals of fluxes J Ca,net,s and J Ca,net,t over the stimulation period (T) or as a sum of integrals of all ionic currents taking a share in transmembrane Ca 2+ transport through each membrane pool.

The Effects of Extracellular Ion Concentration Changes Under Steady-State and Non-Steady-State Conditions
Under steady state at any frequency, the mean ion concentrations in the intercellular clefts are equilibrated with concentrations in the bulk. Not surprisingly, inclusion of the clefts into the model only minutely affected the simulated results of intracellular Ca 2+ concentration changes ([Ca 2+ ] i and [Ca 2+ ] NSR ) at steady state ( Figure 3). When exploring steady-state effects in the model, ionic currents and the quantities describing Ca 2+ turnover remained practically unaltered when ionic concentrations in the intercellular clefts were fixed at bulk concentrations (not illustrated). In contrast, in the t-tubules, the mean Ca 2+ concentration differs from the bulk concentration as a consequence of the sustained flux, J Ca,net,t (cycling between the tubular and surface membrane [16]). In addition, activity-induced steady-state oscillations in the clefts are relatively small as compared to those in the t-tubular lumen (Figure 2), due to the larger volume of the intercellular cleft and the smaller Ca 2+ channel density at the surface membrane. Thus, the currents mediating transmembrane Ca 2+ movements in the course of each cycle at steady state are more affected by a depletion of extracellular Ca 2+ (which reaches a maximum within one cycle during activation of I CaL ) at the t-tubular membrane than at the surface membrane.
Although the fractions of the main ion transporting proteins localized in the tubular membrane have been quantified in some experimental animals (for reviews, see [31,32]), similar data related to human cardiomyocytes are still missing. Therefore, we investigated a possible influence of varying the fractions of Ca 2+ transporters in the t-tubular membrane on the modulation of intracellular Ca 2+ -transient (represented by Δ[Ca 2+ ] i,area,rel ). The graphs in Figure 4 demonstrate the steady-state effect for two values of f CaL,t assessed in guinea-pig (0.64) and rat (0.8) cardiomyocytes. As I CaL,t causes Ca 2+ depletion in the t-tubules, it is not surprising that a higher f CaL,t corresponds to a lower amount of Ca 2+ in the intracellular compartments reflected by Δ[Ca 2+ ] i,area,rel . Conversely, raising fractions f NaCa,t and f pCa,t attenuates changes in Ca i -transient, because I NaCa,t and I pCa,t tend to enhance the mean value of [Ca 2+ ] t . The simulations suggest that the changes in intracellular Ca 2+ -transient may be significant if, simultaneously, f CaL,t was high (≥0.6) and f pCa,t low (≤0.2).
In contrast, under non-steady-state conditions, when sudden changes in heart rate were simulated, the intercellular cleft space significantly affected Ca 2+ handling and APs (Figures 5 and 6). This was caused by a significant transient depletion of extracellular Ca 2+ concentrations in the cleft space and, consequently, in t-tubules. This process resulted from the initial fast increase of Ca 2+ input into the cell (within several stimuli), since, during each shortened stimulation period after transition to 2.5 Hz, Ca 2+ input via I CaL transiently exceeded Ca 2+ output via I NaCa and I pCa . Restoration of equilibrium between the clefts and the bulk space required tens of seconds.
Transient accumulation of extracellular K + in response to sudden changes of stimulation frequency (Figure 7) was evidenced as parallel changes in resting membrane voltage in agreement with experimental results [21]. It also played a role in the reduction of Ca 2+ -transient (not shown), but its effect was substantially smaller than the effect of Ca 2+ depletions. The effects of minute relative variations of [Na + ] i were negligible.

Cardiomyopathies Alter the Function and Structure of Tubular System
As documented by numerous studies, cardiac diseases, like ischemic or dilated heart failure, affect the structure and function of the t-tubular system, both in animal [33][34][35] and human cardiomyocytes [36][37][38] (for a review, see [39]). The myocardial remodeling includes a reduced number of t-tubule openings [38,40], altered geometry (diameter and length) [36,41] and redistribution of ion transporting proteins between the surface and tubular membrane, namely L-type Ca-channels and Na/Ca-exchange proteins [42]. As a consequence, alteration in Ca 2+ turnover and EC coupling can be expected. Intercellular clefts can also be modified by, e.g., hypertrophic cardiomyopathy [43].
The purpose of the present model study was to estimate the impact of ionic changes in extracellular spaces on APs, ionic currents and intracellular Ca 2+ handling in normal human cardiomyocytes at different frequencies. Our aim in a future work will be to utilize the present model to investigate the consequences of the above-mentioned changes of ion concentrations in extracellular spaces occurring under pathological conditions.

Inhomogeneities Along the T-Tubules
The distribution of ion channels along the t-tubules is considered uniform in our model. However, the tubular system comprises both transverse and longitudinal extensions, and ion channels are probably unevenly distributed among these two membrane pools [17]. Further, channels are likely clustered in specialized regions [44]. The concentration of Ca 2+ channels, at the limited region of the t-tubules that contribute to the dyads, would result in larger local depletions in luminal Ca 2+ , as modeled using partial differential equations [45].

Deformation of the T-Tubules
Deformation of the t-tubules during activity [46] is likely to modulate the magnitude of ionic currents through mechano-sensitivity. In addition, the t-tubules' diameter is narrowed under high blood pressure or under contracture [47]. Such a change in volume during a cardiac cycle would cause fluid exchange between the t-tubule lumen and the interstitial space and help in damping concentration changes. The magnitude of this effect still remains to be evaluated for human ventricular myocytes.

Cellular Subtypes and Regional Variations
Subtypes of ventricular cardiac myocytes (epicardial, mid-myocardial or endocardial) with different ion channel densities were not considered in our study, as their description in human cardiac ventricular myocytes is not fully documented. It was also shown that the amount of t-tubules-SR junctions (but not the volume fraction of t-tubules) varied among different regions of rabbit ventricular myocardium that displayed a difference in AP configuration [48,49]. For these reasons, this study did not aim to provide a precise quantitative view of the behavior of human cardiac cells, but rather, an average representation that will be refined when a more detailed description of human ventricular cells and tissue becomes available.

Methods
The structure of the model (Figure 10) is similar to that of our recent models incorporating t-tubules. From the two existing species-specific models including t-tubular spaces (rat [16] and guinea pig [14]), the properties of guinea pig cardiomyocytes are nearer to the human one. Therefore, we started with our guinea pig model that we then modified, with respect to the recently published morphological data [5,15], the known properties of ion transfer mechanisms [4,5], the configuration of action potentials (APs) [25] and the changes of intracellular ion concentrations in human ventricular cardiomyocytes [27,28]. , sodium-activated potassium current (I K(Na) ), calcium-activated non-specific current (I ns(Ca) ), sodium-calcium exchange current (I NaCa ), sodium-potassium pump current (I NaK ) and calcium pump current (I pCa ). The intracellular space contains the cytosol (i), the dyadic space (d), network and junctional compartments of sarcoplasmic reticulum (NSR, JSR) and endogenous Ca 2+ buffers [calmodulin (B cm ), troponin (B htrpn , B ltrpn ) and calsequestrin (B cs )]. The t-tubular space is partitioned into nine concentric cylindrical segments (details are given in the Supplementary Information). J up represents Ca 2+ flow via SR Ca 2+ pump and the small filled rectangles in JSR membrane ryanodine receptors. The small bi-directional arrows denote Ca 2+ diffusion. Ion diffusion between the t-tubular and cleft spaces is represented by the dashed arrow and between cleft and external bulk spaces by thick arrows.
The principal modifications include: (i) replacement of the description of I Na , I CaL and I K1 by the formulations used by Iyer et al. [5] and those of I Kr , I Ks and I Kto by the simpler formulations of Ten Tusscher et al. [4]; (ii) reformulation of intracellular Ca 2+ handling mainly on the basis of its description in [5]; (iii) partitioning of t-tubular space along its radius into nine concentric cylindrical segments [16] and readjustment of geometrical parameters of t-tubules to meet the recent experimental data from human hearts [15]; and (iv) formulation of a single-layer cleft space and ion exchange with an external bulk space on the basis of the data in [2].
Because of the paucity of experimental data from human ventricular myocytes, the time constants describing ion exchange between t-tubular and cleft spaces and also t-tubular fractions of ion transporters related to I Na , I CaL , I K1 and I pCa were preliminarily set to be consistent with data obtained from guinea-pig [14]. The t-tubular fractions of other transporters were set to 0.56 (the fraction of membrane area in t-tubules) on the assumption of their uniform density over the whole cell membrane. The model using these settings is referred to as the "basic model" in the text. Some parameters, so far undefined for human ventricular myocytes, such as tubular fractions of L type Ca 2+ -channels (f CaL,t ), Na + -Ca 2+ exchangers (f NaCa,t ) and the sarcolemmal SR Ca 2+ pump (f pCa,t ), were, however, considered as variable in some simulations to evaluate their importance. All details related to the development of the model and its full mathematical description are given in the Supplementary Information (model parameters are summarized in Figure S1 and Tables S2-S10).
The model was implemented in MATLAB 7.2 (MathWorks, Natick, MA, USA), and the numerical computation of the system of 98 nonlinear differential equations was performed using the solver for stiff system ODE-15s. The stability of the model with time is demonstrated in Figure 11. To achieve steady states in Figures 1-4, the model was run for 1200 s at each stimulation rate. To compare the model results with the behavior of other models of human ventricular cells published so far [5,6,8], we used the computational environment for cellular modelling, COR v. 0.9.31.1409 (Dr. Alan Garny), and the available CellML codes of the models.

Conclusions
The available data related to human cardiac ventricular cells were used to design the first mathematical model of human ventricular cardiomyocyte incorporating t-tubular and intercellular cleft spaces. The model was used to study the physiological consequences of activity-related ion concentration changes in these extracellular spaces. Some of them (particularly, the changes of the intracellular Ca 2+ -transient) were explored as a function of a few crucial parameters, still not quantified in human cardiac cells. In the present state, the model visualizes the activity-related accumulation-depletion of Ca 2+ and K + ions in t-tubular and interstitial extracellular spaces. The results suggest that cumulative calcium depletion in the limited extracellular spaces can affect transmembrane Ca 2+ currents, the duration of AP and the cytosolic calcium-transients that govern cellular contraction. Further experimental work is needed to evaluate the missing values describing the distribution of ion transporters in the cell membrane and the rate of ion exchanges between the extracellular compartments of human cardiomyocytes.