Multiscale stochastic modeling for calcium dynamics in cardiac electrophysiology: assessing whole-cell model reliability under phosphorylation and LCC downregulation

Intracellular calcium (Ca 2+) dynamics drives contractile function in cardiac myocytes. In particular, L-type Calcium Channels (LCCs) and Ryanodine Receptors (RyRs) are organized in microdomains, where LCCs trigger substantial Ca 2+ release from the Sarcoplasmic Reticulum (SR) via RyRs. Different microdomains can be coupled at different length scales by calcium diffusion or common external activation. We present a Scalable Aggregate Calcium Release Unit (SA-CaRU) model for human ventricular myocytes that integrates a recently developed Markov Chain (MC)-based description of LCCs, replacing classical Hodgkin–Huxley gates. Our approach is based on previously published MC-based frameworks for the human heart, enabling stochastic gating and reproducing evoked local Ca 2+ release statistics across different effective levels of microdomain aggregation. Our single-SA-CaRU system captures, within a unified framework, key features of microscale and macroscale Ca 2+ cycling and allows, for the first time, systematic exploration of variability in SR Ca 2+ release as a function of effective microdomain size and coupling. Simulations with increasing numbers of channels reveal that the transition from stochastic to deterministic-like Ca 2+ behavior is typically sharp at a specific cluster size. Under normal (healthy) conditions, this occurs at O (1 0 2) LCCs (with mild sensitivity to the RyR:LCC scaling). However, under high phosphorylation or LCC upregulation, stochasticity persists and convergence to deterministic-like behavior is absent or markedly delayed even for total LCC numbers as large as 20,000. In these conditions, whole-cell deterministic models become doubtful, since their behavior can be qualitatively different from that arising from any plausibly mediated coordination of subcellular calcium release units.

The heart is a mechanical pump and contraction is regulated by an electrical signal that travels from cell to cell, known as the action potential (Malmivuo and Plonsey, 1995). The disruption in the synchronization of these electrical impulses can lead to various arrhythmias, including atrial flutter, atrial fibrillation, ventricular fibrillation, and ventricular tachycardia (Zipes et al., 2017). Some of such arrhythmias are fatal, making cardiovascular disease the leading cause of death in developed countries (Mendis et al., 2011; Timmis et al., 2020).

Clinical studies have shown that the trigger for this arrhythmia can be localized excitation, which can lead to premature ectopic beat. This abnormal excitation can appear due to micro-reentries (Oliveira et al., 2018; Alonso et al., 2016) or due to the interaction of abnormal calcium cycling and sodium currents in cardiac cells (Clusin, 2003; Wehrens et al., 2005; Xie et al., 2015). Such suggestions have been further studied using mathematical models of cardiac tissue (Alonso et al., 2025; Niederer et al., 2019), which basically describe the electrophysiological properties of cardiac cells and, in particular, the dynamics of the action potential. Such models have suggested that malfunctions in calcium cycling can trigger ectopic beat action at the single-cell level, both as early afterdepolarizations (EADs) and as delayed afterdepolarizations (DADs) (Colman et al., 2022; Chen et al., 2009).

Cardiac contraction is based on highly regulated intracellular Ca 2+ cycling. At the onset of each heartbeat, upon cell depolarization, L-type Calcium Channels (LCCs) open and allow moderate Ca 2+ influx through the membrane, triggering the opening of RyR2 in clusters in front of the LCC initiating a dramatic release of Ca 2+ from the sarcoplasmic reticulum (SR) in a process known as Calcium-Induced Calcium Release (CICR). The calcium released must be later reuptaken into the SR via the SR Ca pump (SERCA) and outside the cell, mainly via the sodium-calcium exchanger (NCX). The joint structure of the RYR2 receptor cluster and the LCC in front of them, while each is in a different membrane, is tightly packed with a small dyadic volume between them. They represent a coherent calcium release unit (CaRU) at the core of heart contraction. These microdomains are distributed along the Z-planes of the cardiomyocyte and associated with invagination of the cell membrane, t-tubules, with a characteristic distance between them of 300–500 nm within the plane and a wide distribution of sizes. Between the Z-planes, roughly 2 μ m apart, some RyR2 clusters without associated LCC are also present, but they do not constitute the core of the normal calcium transient and are very sparse.

Although the core of cardiac Ca 2 + cycling originates in discrete microdomains, modeling efforts have often operated at the whole-cell scale using deterministic models that treat the whole cell as single set of typically five compartments, such as the ones developed by Pandit et al. and ten Tusscher and Panfilov (Pandit et al., 2001; Ten Tusscher and Panfilov, 2006), to simulate the global behavior of action potentials and calcium transients. Necessarily, they treat ionic currents through population-level averages using Hodgkin-Huxley-type formulations, and the stochastic opening and closing of individual ion channels and the resulting local fluctuations are not considered or studied. Subcellular behavior such as the rate of stochastic local releases (sparks) during diastole or the possible presence of diastolic calcium waves (Sato et al., 2013; Colman et al., 2022), where the probabilistic interaction between a small number of L-type Ca 2 + channels (LCCs) and neighboring ryanodine receptors (RyRs) within the dyadic cleft plays a crucial role, is not accounted for. Any relevance they might have to the global average dynamics of a normal heartbeat is ignored. This means that other possible sources of pro-arrhythmic behavior may be outside of the scope of these models (Colman et al., 2022).

To address these local phenomena, subcellular models have emerged that treat the calcium release unit (CaRU) as the fundamental computational element, incorporating thousands of discrete CaRUs where a small number of LCCs and RyRs are present (Marchena and Echebarria, 2018). These units are coupled directly by calcium through diffusion and indirectly via a common membrane voltage. Their stochastic dynamics require the use of either Monte Carlo techniques or Markov Chain (MC) state diagrams. Such approaches can be used to reproduce critical features such as Ca 2 + sparks, brief localized calcium releases from the SR, as well as spark-mediated Ca 2 + waves that can propagate across the cell.

Bridging the gap between these two modeling paradigms—whole-cell deterministic and subcellular stochastic-is a major challenge. The behavior of individual CaRUs does not scale trivially to the whole cell level, as spatial organization, local calcium diffusion, and the nonlinearities of CICR introduce emergent behaviors not predictable from averaged dynamics. The computational burden of simulating thousands of stochastically coupled release units with spatial resolution is significant, making the simulation of whole organ behavior using subcellular models numerically impossible. Proper whole-cell organ models require an understanding of the relationship between the two approaches. More specifically, understanding when the average behavior of calcium release units can be properly translated into whole-cell dynamics is crucial. When local and global behaviors diverge under specific conditions, whole-cell deterministic models must be carefully adapted for studies at the whole-heart scale.

In this manuscript, to address the stochastic behavior of the calcium-induced calcium release (CICR) process, we define a Scalable Aggregate Calcium Release Unit (SA-CaRU) in which both LCC and RyR channels are modeled as MCs. We simulate the dynamics of a single cluster and systematically analyze the variability of Ca 2+ release as a function of the number of channels within the SA-CaRU (i.e., cluster size), systematically varying the number of LCC and RyR channels from a single LCC channel and 5 RyRs to thousands of channels within the same SA-CaRU. An SA-CaRU with a low number of LCC reproduces the behavior of a real cluster of RyR, a real Calcium Release Unit. For slightly larger values of LCC, it captures the behavior of macrodomains coupled by diffusion, and for a large number of LCC, it reproduces a mean-field analysis with unrealistically high diffusion, allowing us to check whether it captures average whole-cell behavior. We want to stress that the SA-CaRU does not represent a single enlarged microdomain when the number of channels is very large. Rather, it is a mean-field statistical description of many independent CRUs, analogous to population-density approaches. No assumption of instantaneous or perfect local coupling is made; CRUs interact only through the global cytosolic Ca concentration. This formulation preserves known spatial constraints of dyadic microdomains while capturing the collective stochastic behavior of CICR, allowing the study of the system at different scales.

We consider both the original parameters and two sets of key parameters that emulate the common adjustments made in whole-cell models. The first set mimics conditions observed in patients with marked downregulation of LCCs, while the second represents increased pumping activity mediated by phosphorylation. The robustness of the relationship between deterministic whole-cell models and the subcellular coupling among a few microdomains under these conditions provides important insights into the adequacy and limitations of whole-cell approaches.

Our results reveal a sharp behavioral transition around 100 L-type calcium channels (LCCs) in healthy cardiomyocytes, marking a critical threshold where calcium release shifts from a stochastic to a more deterministic regime. Below this threshold, release events display high variability, primarily driven by random fluctuations in the opening and closing of individual channels. In contrast, under LCC downregulation, and especially under conditions of strong phosphorylation, the distinction between deterministic and stochastic regimes persists even when more than 10,000 channels are involved. In these cases, deterministic whole-cell models fail to reproduce physiologically realistic behavior. Instead, the subcellular system exhibits intrinsically stochastic dynamics characterized by a bimodal distribution, alternating between large and very large calcium release events that can only be captured by an explicitly stochastic subcellular model.

We begin with a baseline deterministic cell model and introduce stochastic Markov Chain representations for L-type Calcium Channels (LCCs) and Ryanodine Receptors (RyRs). This approach enables us to adapt subcellular calcium compartments and fluxes to exhibit both deterministic and probabilistic behaviors depending on the channel densities. We can define in this way a Scalable Aggregate calcium Release Unit (SA-CaRU) that, depending on he number of LCC channels and RyRs receptors, will move from a purely stochastic behavior at a very low number of channels towards a full deterministic model when the number of channels goes to infinity. We proceed to describe the details of this general idea.

Electrophysiological model

We take as our original model, which we then modify, the ten Tuscher model for the human ventricle cell (Ten Tusscher and Panfilov, 2006). We replace the original Hodgkin-Huxley (HH)-type gating formulations with a novel representation based on Markov Chain (MC), which produces an equivalent response to experimental conditions (Novaes et al., 2022).

More specifically, in the base model used in this work, since cardiac cells are separated from the neighboring extracellular space by a phospholipid bilayer membrane, the cell is considered a capacitor that stores electric charge. The total electrical current through the cell membrane is equal to the change in the charge of the membrane, giving rise to the standard single-cell models characterized by the change in the membrane potential: where V = V i − V e in Equation 1 is the electrical potential of the membrane, C m is the capacitance of the membrane and the term I s t consists of an external stimulus or pacing. In our model, we adopt a popular formulation of the total ionic current (Ten Tusscher and Panfilov, 2006), which includes contributions from sodium, potassium, and calcium currents. Mathematically, it is defined as: where the terms in Equation 2 represent ionic fluxes through voltage-gated channels, pumps, and exchangers, see Figure 1. To provide a full cell model, sodium, potassium, and calcium cycling have to be included to track the concentration of ions inside the cell. So, in order to complete the single-cell description, intracellular ionic concentrations are dynamically computed, with particular emphasis on the calcium-handling system, which is central to the excitation–contraction coupling mechanism and, therefore, to this study.

We change the I CaL and the calcium cycling of the original ten Tusscher model (Ten Tusscher and Panfilov, 2006) in terms of a Scalable Aggregate Calcium Release Unit (SA-CaRU). More specifically, we modify the L-type calcium current (I CaL) and the calcium cycling subsystem of the original ten Tusscher model to incorporate a Markov Chain–based representation leading to a stochastic framework (Novaes et al., 2022) which allows us to describe the probabilistic opening and closing of L-type Calcium Channels (LCCs) and their interaction with Ryanodine Receptors (RyRs) within a Scalable Aggregate Calcium Release Unit (SA-CaRU). The SA-CaRU provides a unified structure that encompasses the possibility to simulate channel-level stochasticity and whole-cell-level calcium signaling, allowing us to investigate the gap between local microdomain behavior and global calcium cycling. The general structure of the model is provided in Figure 1 and we proceed to provide details.

Scalable aggregate CRU setup

Our focus is to study an Scalable Aggregate Calcium Release Unit (SA-CaRU), composed of a cluster of LCCs and RyRs governed by a Markov Chain process within a dyadic domain. In our framework, the SA-CaRU does not represent only a single anatomical microdomain; rather, it serves as a flexible abstraction that can range from an individual subcellular Ca 2+ release unit to a coordinated cluster of microdomains, and even up to a cell-wide release structure. This allows us to investigate calcium dynamics on different spatial and organizational scales using a unified modeling approach.

Let N LCC and N RyR denote the number of L-type calcium channels and ryanodine receptors in the cluster, respectively. A typical single microdomain configuration might use (N LCC, N RyR) = (4,20) (Bers and Stiffel, 1993; Qu et al., 2022). However, changing and increasing the number of LCC and RyR, we can explore how stronger coordination between microdomains works. Our model offers flexibility in fully exploring up to very large values of N LCC and N RyR, which should reproduce the actual models of single-cell behavior.

In this single unit, local calcium concentrations, such as [ C a ] s s (calcium in the subsarcolemma), [ C a ] S R (calcium in the Sarcoplasmic Reticulum), and [ C a ] i (cytosolic calcium), are determined by fluxes resulting from I CaL, I rel, and I u p, similar to standard deterministic modeling approaches for the total cell. The governing equations, see Equations 3–9, for these variables read:

Fluxes related to NCX activity I NaCa, calcium pump I pCa, calcium background current into the cell I bCa, and the leak of calcium from the SR into the cytosol I leak follow the same formulation as in the original ten Tusscher model. They read:

The different parameters in Equations 10–13 also come from the ten Tusscher model. For more details, see the original model (Ten Tusscher and Panfilov, 2006). We proceed now to describe the other key currents.

L-type calcium channel

The original ten Tusscher model uses a configuration of this channel based on HH-gated variables, which does not have an equivalent stochastic counterpart. Therefore, we define the L-type current as in (Novaes et al., 2022), where each channel can be in seven different configurations represented by a seven-state Markov Chain, see Figure 2, and was designed to represent the dynamics of the original ten Tusscher model. Therefore, the seven-state L-type Ca channel Markov model employed in this study was previously introduced and fully validated (Novaes et al., 2022). That work demonstrated that the model reproduces canonical I–V relationships and ion channel dynamics. When embedded in the ten Tusscher ventricular framework, it also preserved action potential morphology, Ca transient amplitude, restitution behavior, and Ca cycling stability across pacing rates. Because the present manuscript focuses specifically on the stochastic-to-deterministic transition, we use the same validated formulation without modification. There are also different types of models of the cardiac L-type calcium current (Agrawal et al., 2023).

The channel transitions between states in a stochastic manner. The key state responsible for generating the current is the open state O L C C, as depicted in Figure 2. When the channel occupies this state, the resulting current is given by: where Δ V = V − 15 in Equation 14. When the channel is in any other state different from the open state, it does not generate current. When in the open state current is produced.

We must emphasize that the transitions between the LCC states are stochastic and that, in the limit of a very large number of channels in the SA-CaRU we revisit the deterministic whole-cell model. It will be important to track the fluctuations when we consider a low number of channels. For this, we use Tau leaping to efficiently model systems with discrete events. This is achieved by exploiting the binomial distribution to represent the number of successful transitions of states (Tian and Burrage, 2004). For each transition with rate k, the number of channels that changed state during a time step Δ t was drawn from a binomial distribution with probability of success p = 1 − e − k Δ t. This approach guarantees non-negative populations and avoids excessive-leap instabilities associated with classical tau-leaping. All simulations used Δ t = τ = 0.001 ms, which satisfies k Δ t ≪ 1 over the full range of channel rates. Using this approach, the model captures the stochastic dynamics of the system while significantly reducing the computation time compared to simulating each event individually.

Transitions implement the voltage-dependent opening, calcium-dependent inactivation, and closed or inactivated states discussed in other models [as in rabbit (Mahajan et al., 2008)]. This formulation has been shown to adequately reproduce the deterministic evolution of a healthy human ventricular cardiomyocyte (Novaes et al., 2022).

Section 2.1 in the Supplementary Material provides a description of the different rates that define the stochastic process in the L-Type calcium channel. Most of the parameters were taken obtained using a genetic algorithm to search for a proper reproduction of the whole-cell model in the human ventricle (Novaes et al., 2019).

The basic structure of these parameters was designed to provide a more biophysically faithful representation of how a modest trigger current can initiate a large-scale release of calcium from the SR while keeping the basic properties already adequately described in whole-cell models.

Finally, we will compare the stochastic behavior with the exact deterministic behavior coming from this set-up. In this case, each of the Markov Chains becomes an occupation number that tracks a set of differential equations. These differential equations are reproduced numerically.

Calcium release from the SR into the cytosol. The RyR2 model

Calcium release flux through RyR channels is activated when the receptor is in the open state and is modeled proportional to the calcium concentration gradient. The flux is given by: where the variable O R y R in Equation 15 indicates whether the RyR is in the open state. We used a RyR model with four different states as shown in Figure 3, one closed, one open, and two inactivated. This configuration of the calcium release flux model through RyR was already used in the original model (Ten Tusscher and Panfilov, 2006).

The four-state model of RyR has cytosolic calcium regulation corresponding to the calcium-induced calcium-release mechanism (CICR) and inactivation is also directly mediated by cytosolic calcium, mimicking the regulation of the Ca-calmodulin complex (Wei et al., 2021). Furthermore, in this RyR model, both the transitions from the resting state to the open (O) state or the resting-inactivated (RI) state are also regulated by the calcium concentration of the sarcoplasmic reticulum (SR), capturing luminal control of channel gating (Shannon et al., 2004).

Table 1 provides the list of parameters needed to perform the stochastic simulation of the states of the Markov chain. Again, to simulate the number of RyR in the open state, we stochastically computed the number of multiple events that occur within a specified time interval (Δ t) using a Tau leaping approach as in the previous section. We will also need to compare with the full set of equations for calcium cycling that is provided in the Appendix.

Modeling of calcium cycling alterations

We implemented two modifications in the parameters of the SA-CaRU model by targeting two sets of key parameters which aim to reproduce two characteristic features of, first, calcium regulation via LLC regulation due to genetic diseases and, second, calcium regulation via specific phosphorylation.

We use a single parameter indicator ϕ to indicate changes in phosphorylation and another single parameter η for up or downregulation of LCC channels. We describe how we propose to implement these changes.

LCC regulation

The lower expression of LCCs, corresponding to downregulation, decreases the inward calcium current, which is crucial to trigger calcium release from the sarcoplasmic reticulum (Shaw and Colecraft, 2013). The decrease in calcium influx has important consequences for cardiac cell excitability and results in altered cell contractility (Triggle, 1998) and disturbed electrical activity. Finally, such downregulation alters the configuration of the action potential, influencing the overall cardiac rhythm (Benitah et al., 2010).

Increased LCC function, corresponding to upregulation, induces an elevated inward calcium current in cardiac cells, which enhances calcium-triggered calcium release from the sarcoplasmic reticulum. This increase in channel open probability can increase the probability of atrial fibrillation (Klein et al., 2003). The increase in cytosolic calcium transients can improve the strength of cardiac muscle contraction (Zhao et al., 2009). Finally, high LCC activity prolongs the plateau phase of the cardiac action potential and can increase the probability of arrhythmia (Satin and Schroder, 2009).

We model the regulation of the LCC channels by the modification of the inactivation rates, see Figure 2: where we have included the extra parameter η. We have simulated the decrease η = 0.5 and the increase η = 1.5 of the inactivation rates, i.e., up and downregulation of the LCC channels, respectively, in Equations 16, 17. The modification of these two rates is also equivalent to the increase and decrease of the rate d +, see Figure 2, as implemented in Sato et al. (2018) to model the cooperative gating of L-type Ca 2 + channels (Sato et al., 2018).

Effect of phosporylation

Protein kinase A (PKA) regulates the cardiac excitation–contraction coupling because it can be activated by cAMP and phosphorylate critical proteins in the channels LCC, RyR, and, phosphorylates also Phospholamban (PLB), aftecting SERCA. Therefore, PKA has a important influence on calcium dynamics and the contractility of the myocytes. Alterations in RyR, LCC and SERCA phosphorylations play a critical role in various cardiac diseases, including heart failure.

The effects of PKA are threefold. First, PKA prolongs the opening time of LCC and increases calcium influx into the cell (Ahern et al., 2019; Ahern et al., 2021). This modification of LCC activity increases SR Ca 2 + release, promotes Ca 2 + -dependent remodeling and prometes spontaneous SR Ca 2 + -release events, To model such increase, we multiplied the transition rate from the C ′ to O state - t C ′ → O - of the LCC MC (see Figure 2) by a factor ϕ similar to some downregulation approaches (Sato et al., 2018). So, the updated version of the transition rate t C ′ → O reads

Second, PKA also affect the behavior of the RyR2 channels, since the RyR2 has some important phosphorylation sites (Dobrev and Wehrens, 2014) that change the stability of the different states upon the binding of the phosphate group. The key effect of this phosphorylation is to increase the opening rate probability of the RyR2 increasing the probability of reaching the opening state as described in (Warren, 2021; Wehrens et al., 2006; Respress et al., 2012).

We implement this effect by changing parameter k 1, (see Table 1), by ϕ. This modifies the rate of opening, and when ϕ > 1 it increases the open probability. So, the updated version of the parameter k 1 reads:

And, finally, Phospholamban (PLB) is a regulatory protein that inhibits the activity of the Sarcoplasmatic reticulum calcium ATPase (SERCA). PKA phosphorylation of PLB decreases SERCA inhibition, which results in faster calcium recapture (Tomek and Zaccolo, 2023). The Sarcoendoplasmic reticulum calcium ATPase (SERCA) uptake reads where the effective term, K u p / ϕ, takes into account different possible affinities resulting from PKB regulation.

Notice that, in this way, all the changes produced by phosphorylation are introduced in ϕ in Equations 18–20 as a single parameter, changing the three key targets of phosphorylation. More complex relations could be established by changing the effect in each case, but we propose that a common value is physiologically relevant and a very clear way to investigate different phosphorylation levels.

In this manuscript, we consider two phopshorylation level scenarios relative to control (ϕ = 1). In scenario 1, parameter ϕ was reduced by 50% (i.e., ϕ = 1 ↦ ϕ = 0.5) representing a drop in phosphorylation. In scenario 2, the same parameter was increased by 100% (i.e., doubled: ϕ = 1 ↦ ϕ = 2) representing an increase in phosphorylation.

Summary of computational experiments

As described in the previous sections, this study introduces a stochastic Markov Chain (MC) model in which the cell is represented by a single Scalable Aggregate Calcium Release Unit (SA-CaRU). In the limit of a large number of L-type calcium channels (LCCs) and ryanodine receptors (RyRs), the model behavior converges to that of the deterministic system. However, in finite and physiologically relevant channel numbers, the stochastic formulation captures variability in calcium fluxes, particularly in the L-type calcium current (I CaL) and the release current from the sarcoplasmic reticulum (I rel).

To investigate how many channels are required for the system to approximate the deterministic limit and to quantify how calcium dynamics vary between stochastic and deterministic-like regimes, we conducted numerical simulations in which the number of LCCs in the SA-CaRU was systematically varied. Specifically, we explore N LCC ∈ { 1,2,4, ⋯, 32768 }, increasing the count by powers of two (i.e., up to 2 15).

The purpose of this scaling was to evaluate how the degree of channel clustering affects the emergence of coordinated calcium signaling. To isolate this effect, we maintained a constant ratio between LCCs and RyRs throughout all simulations. Based on anatomical and functional estimates from the literature (Bers and Stiffel, 1993; Qu et al., 2022), we fixed N RyR = 5 ⋅ N LCC for each experiment.

For each combination of N LCC and N RyR, we simulated 1100 pacing cycles at a frequency of 1 Hz. The initial 100 beats were discarded to eliminate transient dynamics and ensure that the system reached steady-state conditions. The remaining 1000 beats were used for analysis. This setup enabled us to capture the beat-to-beat variability intrinsic to the stochastic formulation and to study its dependence on the channel population size.

To quantify the degree of stochasticity in our simulations, we focused on the calcium release current from the sarcoplasmic reticulum, I rel, representing the flux of calcium ions into the cytosol. For each beat, we computed the total calcium released over a 1-s interval by summing I rel (t), yielding a scalar quantity that characterizes the amplitude of release. Mathematically, it reads where N in Equation 21 is the number of Δ t steps used to simulate one period.

This summation was performed for each of the 1000 beats simulated per channel configuration, producing a distribution of S ̄ (I rel) values. We used descriptive statistics to characterize the variability of this distribution and to compare the degree of stochasticity between simulations with different numbers of LCCs and RyRs.

To visually assess the variability in calcium release under different channel configurations, we generated histograms of the summed release signal S ̄ (I rel) for each stochastic simulation. These histograms capture the distribution of calcium release amplitudes across 1000 steady-state beats for each N LCC, and reveal how stochastic fluctuations evolve as the number of channels increases. These visualizations offer qualitative insight into the degree of stochasticity present in the system and complement the statistical descriptors of variability.

To contrast this behavior with the idealized deterministic regime, we also analyzed the last 10 pacing cycles of the deterministic version of the model. For each of these cycles, we extracted the temporal evolution of key variables, including membrane potential V, intracellular calcium concentration [ C a ] i, subspace calcium [ C a ] s s, L-type calcium current I CaL, and calcium release current I rel. These traces reflect the steady-state behavior of the system in the absence of stochastic fluctuations and serve as a visual reference against which the effects of noise can be qualitatively assessed.

The output data were produced for each of the five simulation scenarios investigated: the baseline physiological condition, the high and low phosphorylated state (representing β -adrenergic stimulation) and the LCC up and downregulation mimicking voltage-dependent modulation. This allowed for direct comparisons across pathophysiological conditions and cluster sizes.

The results are organized as follows. First, we present the dynamics of the pulses for the standard configuration of the model as a function of the stochasticity, which is controlled by the size of the SA-CaRU.

Next, we introduce two modifications to the model, which enable the extraction of information not only from the deterministic dynamics, as is standard in scientific research, but also from the dependence on stochasticity under various scenarios of LCC malfunction and phosphorylation.

Dependence of the dynamics on the number of channels in our baseline model

Large number of channels–deterministic limit

As explained in the methodology, we numerically pace the deterministic single-cell model every one second through the channel I s t. After the transitory behavior, we find a repetitive pattern of pulses as seen in Figure 4. The action potential has the usual shape of a ventricle cell (Ten Tusscher and Panfilov, 2006) and the free calcium and the two fluxes determined by I CaL and I rel take the usual shape from both descriptions based in HH-gates, see (Ten Tusscher and Panfilov, 2006) and the corresponding description in terms of Markov chains (Novaes et al., 2022). In every beat, the concentration of ions that leave RyR channels (I rel) remains constant, reaching a steady state.

Finite number of channels–stochastic effects

We also consider the stochastic behavior of the ion channels and employ the description in terms of a stochastic algorithm as explained in the methodology. In the limit of a large number of channels, the behavior of the system has to coincide with the deterministic simulation. For a lower number of LCC, we want to observe the type of behavior obtained. We have measured the total release of calcium through the RyR channels to visualize the different dynamics of the model as a function of the number of channels considered; see Figure 5. As expected, the release converges to the deterministic value when the number of channels is large.

However, the release decreases as the number of channels considered is reduced. For a small number of channels, the average release is smaller and has some variability, as shown by the intervals in Figure 5. Therefore, as expected for a stochastic process, the release of calcium through RyR channels varies from beat to beat. However, the high dispersion of transients at a low number of LCC gradually converges to a Gaussian-type dispersion of values as the number of LCC increases. A remarkable result is that this dispersion is drastically reduced after N LCC = 128, where the values converge very rapidly to the deterministic case. The same reduction in the variability is observed in the I CaL traces plotted in Figure 6. This indicates that if calcium diffusion can coordinate microdomains of up to 200 channels, or around 50 microdomains if we consider that there are 2-8 LCC per cluster, the local coupled-behavior does not exhibit strong stochasticity.

Similar tendencies to the results shown in Figure 5 have been obtained when studying the AP traces (Figure 7). One should be careful here, since the action potential in an SA-CaRU does not accurately reflect a realistic action potential. For an SA-CaRU with a large number of channels, it indeed approximates the mean-field AP, allowing us to assess whether the stochastic behavior converges to the deterministic one. For low (2–16) or intermediate (32–2024) numbers of LCC, however, the action potential is merely a representation of the local potential in the model and must not be understood as a realistic, measurable local potential. Its variability and how it differs from the deterministic gives us a clear picture of when the stochastic effect diminishes. Action Potential Duration at 50% (A P D 50) and 90% (A P D 90) are plotted as functions of the number of LCC channels, see, respectively, Supplementary Figures S1, S2. The complete set of action potential shapes and their corresponding average shapes is shown in Supplementary Figure S3 for several values of the number of LCC channels. The variability of APD and its convergence to the deterministic case are clearly observed.

This is one of the most interesting features we unveil in this manuscript. The transition towards deterministic values is not a smooth transition but a sudden change in the level of noise-induced dispersion at rather low values of LCC channels. In order to better visualize this point and observe the distribution of release for the different pulses, we plot the histograms of this value for the 1000 pulses considered in each simulation for some selected examples of LCC channels in Figure 8 (baseline N RyR = 5 N LCC). The histograms allow for a clearer understanding of the dependence of the distribution on the number of channels. For a very low number of channels (N LCC = 1 − 2), the majority of pulses do not arrive to release a substantial amount of ions, and the value is frequently zero. For a slightly high number of channels, the release increases, although it does not arrive at the deterministic limit. The distribution shows a Gaussian-like distribution around the average number, typical for a stochastic algorithm. For N LCC = 16 a slightly bimodal distribution is observed, with a small second maximum above the value of the deterministic limit. The first maximum drift to the right as a function of the number of channels in a continuous displacement towards the value in the deterministic case. Crucially, around O (1 0 2) channels (in the baseline scaling N RyR = 5 N LCC), the distribution becomes very sharp and centered on the value of the deterministic case. The implications of this finding will be discussed in the discussion.

We further assessed whether this sharp transition depends on the baseline choice of the RyR:LCC scaling. While the main results use N RyR = 5 N LCC, we repeated the same analysis for several ratios N RyR / N LCC ∈ { 1,2,4,8,16 }. The sharp collapse of dispersion persists across all ratios (Supplementary Figure S4), indicating that the existence of a threshold-like transition is a robust feature of the SA-CaRU dynamics rather than a consequence of a specific fixed ratio. Moreover, when the same results are represented as a function of N RyR (Supplementary Figure S5), the curves collapse more tightly, suggesting that the effective number of RyR channels is a primary determinant of the dispersion of the SR release statistic.

In Figure 9 we show the calcium pulses along time for different values of LCC channels. For intermediate values of the number of channels, see N LCC = 256 in Figure 9, the behavior resembles the deterministic case, a sequence of pulses with a relatively constant shape. For lower values of the number of channels, second pulses appear randomly along the sequence of pulses, giving rise to the bimodal distributions observed for N LCC = 16 − 64 in Figure 8 further decrease of the number of channels produces the alternance of typical pulses with a size similar to that observed in the deterministic case, and very small pulses. Such behavior reminds us of the calcium alternans although without its typical periodicity.

Regulation due to modifications on LCC inactivation rates

Large number of channels–deterministic limit

The up and downregulation of the LCCs will increase and decrease, respectively, the influx of calcium and therefore free calcium; as shown in Figure 10.

For upregulation of the LCC, the deterministic simulation of the whole-cell model, reducing the parameter η by 50%, shows a second peak of flux due to I rel; see Figure 10. The second peak, together with an increase in the release in the first peak compared to the control case, doubles the calcium release due to I rel. This second peak of calcium release produces a secondary peak in the action potential, which lengthens the repolarization process.

The increase in the parameter η in the deterministic whole-cell model corresponds to a downregulation of the LCC channel, and produces a reduction in calcium release and in the free calcium available in the cell; as shown in Figure 10. The release resulting from flux I rel is smaller than in the control case. However, the action potential remains similar to that in the control case.

Finite number of channels–stochastic effects

When the individual behaviors of the SA-CaRU are considered, we observe a much richer phenomenon. As expected, for a very large number of channels, we recover the same dynamics of the deterministic limit. If we incorporate the stochasticity from the GRCU we obtain similar dynamics for the downregulation case; see Figure 11, but with smaller values of the total release; see Figure 11. We also obtain a slight increase in dispersion up to N LCC = 128 in comparison with the control case.

For the upregulation case, we find very broad distributions up to N LCC = 16384 channels; see Figure 11. This is substantially different from the previous case and indicates that, under such circumstances, the stochastic models produce different dynamics than the corresponding deterministic case.

Similar tendencies to the results shown in Figure 11 were obtained for A P D 50 and A P D 90, which are plotted as functions of the number of LCC channels in Supplementary Figures S6, S7. The complete set of different shapes of the action potential and its corresponding average shape can be observed for the upregulated in Supplementary Figure S8 and downregulated cases Supplementary Figure S9, for several values of the number of LCC channels. While the downregulated case converges rapidly to the deterministic case, the upregulated case has broad distributions for a finite but large number of channels, in correspondence with the results for the calcium release observed in Figure 11.

The distributions of the stochastic metric S ̄ shown in Figure 11 can be very broad, depending on the number of channels. In Figure 12 we show the corresponding histograms showing the shape of the distribution of the release values for different numbers of channels, for the two cases under consideration in this section, compared to the control case.

For a small number of channels, the distributions are broad, covering the entire spectrum of values. For N LCC = 16 both regulation cases show a bimodal distribution. For intermediate values, while the downregulation case gives rise to localized standard distributions around a clear value, the upregulation case keeps the bimodal distribution until large values of the number of channels, see the case N LCC = 2048 in Figure 12. The broad distributions observed in Figure 11 for the upregulation case can be associated with the effects of the bimodal distributions obtained in Figure 12. Finally, for more than two thousand channels, in the three cases, the distribution collapses into a single sharp peak, recovering the deterministic cases.

Calcium release depends on the number of channels under consideration, and the resulting free calcium is different for the particular time window and the number of channels. In Figure 13, the free calcium is plotted for ten pulses for both regulatory cases and compared to the control case.

For just one single channel, see Figure 13A, the free calcium in the three cases is very small compared to the deterministic case, see Figure 13F.

For an increasing number of channels, the downregulation case shows an alternation between a small and a large peak of calcium. Finally, for N LCC larger than 32, the simulations of the downregulation condition show a stationary shape of the calcium peak similar to the peak in the deterministic case, however, substantially smaller.

For an increasing number of channels, the upregulation case shows, as expected, much larger calcium pulses, and the shape of the pulses presents several peaks associated with multiple releases of calcium. The number of peaks is not constant for an intermediate number of channels and only collapses to a fixed value for the deterministic case, giving rise to the fixed value of S ̄ shown in Figure 11. The presence of such multiple peaks in the same calcium pulse explains the bimodal distribution found in the results shown in Figure 12 for large values of the number of channels and the corresponding broad distributions shown in Figure 11.

Phosphorylation

Large number of channels–deterministic limit

We consider high and low phosphorylation by increasing and decreasing the parameter ϕ in the deterministic model, respectively.

For high phosphorylation, the deterministic simulation shows an increase in the flux from the RyR channels and an increase in the free calcium after each pulse; see Figure 14. However, the action potential does not change substantially. In this case, the release is the same for all the pulses, and they are identical during the whole simulation.

The decrease in the parameter ϕ produces some unexpected effects in the deterministic dynamics of the action potential and the ion currents. Although the action potential does not seem to change in comparison with the control case, see Figure 14, the flux I rel, and in particular the free calcium concentration, show a strong diversity of pulses with clearly different values. These calcium alternans can have important implications on contractibility. Such an alternation of values in the flux from the RyR gives rise to nonhomogeneous values of release. Therefore, in this case, the release value is not a simple value but a distribution of values. Our results on the effect of phosphorylation on cardiac intracellular calcium handling and the link to cardiac alternans corroborate with previous observations on cardiac tissue (Laurita and Rosenbaum, 2008).

Finite number of channels–stochastic effects

When stochastic simulations are included through the SA-CaRU, the dynamics change strongly with respect to the control case; see Figure 15. In both cases, low and high phosphorylation, the distributions are broad for a wide range of numbers of channels. In the case corresponding to high phosphorylation, the distributions are very broad until N LCC = 32768, and only collapse to a single peak for the deterministic case. It is important to stress that such a large coupling is unrealistic. In order to have thousands of LCC channels coupled, it would require that calcium diffusion strongly couple microdomains at very large distances within the cell. It has been well established that diffusion alone cannot synchronize this behavior. As we discuss later, this indicates that stochasticity is important and that the reduction to the deterministic case is questionable.

On the other hand, for low phosphorylation, the simulations show, for the whole spectrum of numbers of channels, a very dispersed distribution of release values; see Figure 15. Contrary to the previous case, we obtain in the case with N LCC = 32768 the same behavior as in the deterministic case. Here, whole-cell models would seem more appropriate to describe changes at the microdomain level than in the previous case.

For phosphorylation, the tendencies of the results shown in Figure 15 are relatively different from the tendencies obtained for A P D 50 and A P D 90, which are plotted as functions of the number of LCC channels, see, respectively, and Supplementary Figures S10, S11. In both cases, the duration of the action potential remains closer to the control case. The complete set of different shapes of the action potential and the corresponding average shapes can be observed for the low phosphorylation case in Supplementary Figure S12 and for the high phosphorylation case in Supplementary Figure S13, for several values of LCC channels. A second peak in the action potential, eventually associated with an ectopic beat, remains for a large number of simulations, see Supplementary Figure S13.

The assessment of the whole-cell deterministic model can be better understood by the visualization of the histograms associated with the simulations. For stochastic simulations, the release follows a bimodal distribution with two separated maxima compared to the conditions studied above; see Figure 16. See the case with N LCC = 2048, shown in Figure 16, where both phosphorylation cases show clear bimodal distributions in comparison with the single peak of the control case.

For the case of low phosphorylation, in the deterministic limit, a continuous distribution of values is observed; see Figure 16. If the number of channels is decreased, this continuous distribution gives rise to a clear bimodal distribution between two values with a maximal probability, as a clear sign of calcium alternans.

Figure 17 shows the calcium pulses for different numbers of LCC channels. While the bimodal distribution observed under high-phosphorylation conditions may be partly attributed to stochastic fluctuations arising from the numerical simulation, the distribution obtained under low-phosphorylation conditions reflects differences in the amplitudes of calcium pulses induced by external pacing. This behavior confirms the presence of calcium alternans across a broad range of channel numbers.

We have integrated an MC-based L-type current model together with a four-state stochastic RyR model into a single SA-CaRU framework for human ventricular myocytes (Novaes et al., 2022), previously fitted to a well-established model (Ten Tusscher and Panfilov, 2006). Our model maintains deterministic whole cell output at high channel counts while effectively capturing stochastic microdomain subcellular release behavior at lower channel densities. In particular, the transition from random to deterministic-like behavior occurs sharply when the number of L-type calcium channels (LCCs) in the unit exceeds the threshold of O (1 0 2) LCCs.

This has significant implications for our understanding of subcellular calcium dynamics, as it clarifies the importance of effective diffusion coupling between microdomains. The coupling between the LCC and RyR microdomains is influenced by several factors, including the characteristic diffusion length, the distance between microdomains, and the amount of buffering present. Higher diffusion rates, closer microdomain proximity, and minimal buffering pathways lead to greater effective coupling. Since a single microdomain contains only a few LCCs, effective diffusion coupling that can synchronize fewer than 50 microdomains results in a high degree of subcellular variance. In this case, whole cell behavior reflects the average of thousands of highly stochastic and independent units. Conversely, if diffusion coupling can synchronize around 100 microdomains, we begin to see distinct domains of deterministic subcellular behavior. Thus, whole-cell calcium cycling can be understood in terms of the coupling of deterministic systems, both in relation to spark and wave activity.

The degree of coupling between microdomains is not arbitrary; it is physically constrained by the maximum diffusion coefficient of calcium ions in the cytosol (D C a ∼ 800 μ m 2 / s). As the system size increases (i.e., a larger number of CRUs), the time required for diffusion limits the ability of the system to synchronize. Consequently, assumptions of infinite or very high coupling become physically unrealistic for large clusters, as there is a strong finite-size dependence on the effective coupling (Alvarez-Lacalle et al., 2015).

This opens an important avenue for future research. Specifically, the levels of coupling among microdomains within a heterogeneous spatial grid become very important, since the average cell outcome can depend crucially on them. The simulation of spark-induced spark activation and wave propagation in heterogeneous grids should be the next step in our analysis. The key insight from this contribution is that if deterministic behavior emerges spontaneously at around 100 channels, then whether the observed level of coupling crosses this threshold is key to understanding the framework of subcellular calcium cycling and how to transition to proper whole-cell behavior.

A spatially distributed array of CaRUs coupled by cytosolic diffusion represents the most realistic setting for spark-induced spark activation and Ca 2 + wave propagation. The present SA-CaRU framework provides a mean-field baseline that quantifies when stochastic dynamics converge to deterministic behavior as the effective degree of aggregation increases. Extending this work to heterogeneous CRU grids with explicit diffusion coupling is a natural next step, and will allow testing whether physiologically realistic coupling strengths can reach the effective synchronisation regime identified here.

With this in mind, we return to the interpretation of N LCC in the SA-CaRU as an effective proxy for microdomain aggregation. An increase in the number of channels can be understood as an increase in the level of coupling among microdomains. A SA-CaRU with a few LCC reproduces the behavior of a single microdomain. Simulating a SA-CaRU of tens of thousands of LCC is basically a cell that is globally coupled by an effective infinite diffusion of intracellular Ca 2 +. Intermediate values represent realistic diffusion occupying tens of microdomains. By analyzing how the behavior of a single Scalable Aggregate Calcium Release Unit (SA-CaRU) evolves as the number of constituent channels increases, we can directly assess the relationship between stochastic microdomain activity and deterministic whole-cell behavior and track if there is a robust relation between deterministic (whole-cell) and the stochastic behavior obtained at a realistic level of microdomains coupling.

This analysis provides a reasonably quantitative basis for evaluating the validity of deterministic whole-cell models. Our units correctly reproduce the behavior of macrodomains when N LCC is low. As the number of channels increases, our mean-field Calcium release Units could accurately capture the behavior of diffusion-coupled units. At very large LCC channels, the results are just a mean-field analysis that allows us to test when the whole-cell model is recovered. If these large N LCC units do not reproduce the deterministic-like behavior, it is highly unlikely that the whole-cell deterministic model can properly reproduce the average behavior of a cluster realistically coupled with diffusion. In a sense, such models rely either on an unphysically high degree of coupling between microdomains, inconsistent with diffusion-limited calcium signaling observed in real cells, or a strongly synergetic behavior between calcium diffusion and synchronization into a single action potential. This would render the model’s outcome fundamentally different from what one would expect. In this situation, whole-cell models would not adequately describe the subcellular level at the whole-heart or tissue level.

Note that the models for RyR and LCC have been generated to reproduce whole-cell properties (Novaes et al., 2022) rather than specific single-channel recordings. Consequently, while the model captures ensemble behavior accurately, the specific dwell-time statistics and microscopic event correlations should be interpreted as a plausible mechanism consistent with global constraints, rather than a uniquely determined solution.

While the quantitative results presented here are influenced by the specific parameterization of the baseline model, the qualitative findings are primarily governed by the Calcium handling subsystem rather than background transmembrane currents. Consequently, we expect these conclusions to generalize to other human ventricular models that share similar balances of LCC, RyR, and SERCA fluxes. However, significant deviations are expected in models representing small mammals, where the much faster SERCA dynamics fundamentally alter the stability and recovery of the calcium transient.

This SA-CaRU-based approach has enabled us to investigate, under varying degrees of stochasticity, the effects of LCC upregulation and high phosphorylation. It allowed us to move beyond the conventional deterministic whole-cell analysis based solely on parameter variation, providing insight into the consequences of increased stochastic behavior. Our results show, for instance, that deterministic simulations under high phosphorylation yield qualitatively different outcomes than their stochastic counterparts, even when a very large number of channels is considered, indicating noise-sustained switching between distinct response regimes. Likewise, stochastic simulations under LCC upregulation display substantial variability even for large channel populations, calling into question the reliability of deterministic whole-cell model predictions in such regimes.

A key mechanistic insight behind this non-convergence is that, under these pathological protocols, the coupled voltage–calcium subsystem operates close to a high-gain threshold, so that stochastic fluctuations can induce beat-to-beat switching between distinct response regimes. In practical terms, individual beats can fall either into a low-release trajectory (weak CICR recruitment) or into a high-release trajectory (strong CICR recruitment), producing non-Gaussian and often bimodal release statistics. This switching is also reflected in the action potential morphology. Under LCC upregulation, the deterministic model already exhibits an EAD/secondary dome, and the stochastic AP ensemble converges only at very large N LCC toward this same EAD-like morphology (Supplementary Figure S8), indicating delayed convergence due to enhanced sensitivity of the repolarization phase to trigger variability. Under high phosphorylation, by contrast, the deterministic AP does not display an EAD-like feature, yet a small late repolarization shoulder persists in a subset of stochastic realizations across a wide range of N LCC (Supplementary Figure S13). This shows that stochasticity can sustain access to alternative repolarization trajectories even when the deterministic solution remains smooth, consistent with a noise-driven switching mechanism near the CICR/repolarization threshold.

There is substantial evidence that LCC dysregulation contributes to the etiology of heart failure (Khalid et al., 2025). Arrhythmic electrical remodeling, often characterized by altered channel density, increases susceptibility to sudden death in patients with heart failure (Erickson et al., 2011). Furthermore, experimental models have shown that such alterations can induce hypertrophy and heart failure in mice (Goonasekera et al., 2012; Van Oort et al., 2010). These pathological conditions, where the number of functional channels may be critical, highlight the importance of stochastic models that can accurately capture the transition between subcellular noise and whole-cell dynamics.

The results presented here may also guide the choice of discretization in mesoscopic subcellular models, depending on the problem under study. Such models do not exactly reproduce all subcellular dynamics but retain aspects of subcellular structure, typically using discretizations of about 2.5 μm (Prudat and Kucera, 2014), 8 μm (Gouvêa de Barros et al., 2015), or 15 μm (Jacquemet and Henriquez, 2009). This approach is analogous to that used in extracellular–membrane–intracellular (EMI) models (Jæger et al., 2021), which employ discretizations around 10 μm, and in Kirchhoff network models, which extend the EMI framework (Jæger and Tveito, 2023).

We should stress that our results indicate a fundamentally different behavior at lower levels of coupling, from intermediate levels up to the large-cell scale. In this regime, whole-cell models would provide a completely inaccurate picture under any realistic assumption of available coupling. This discrepancy is even more severe than other well-known mismatches between subcellular and whole-cell models discussed in the literature. For example, in the case of calcium alternans, as pacing frequencies increase, whole-cell models can reproduce the appearance of an order–disorder phase transition (Alvarez-Lacalle et al., 2015) as a period-doubling bifurcation. In that case, at least, the general behavior remains similar, although important effects may appear at the tissue level. In contrast, in the cases presented here, subcellular and diffusion-limited calcium models would produce fundamentally different outcomes from whole-cell deterministic models. One could argue that whole-cell models should not be used to reproduce whole-heart or tissue behavior under any malfunction unless an assessment of the type presented here has been previously carried out.

In conclusion, integrating the stochastic L-type calcium channel and four-state RyR models within the SA-CaRU framework provides a versatile platform for bridging subcellular stochasticity and deterministic whole-cell calcium dynamics. By explicitly quantifying how the transition from stochastic to deterministic behavior depends on the number of channels, this approach enables a systematic evaluation of the conditions under which deterministic modeling remains valid or stochastic effects become essential. The SA-CaRU-based methodology thus offers a powerful means to elucidate how microscopic fluctuations influence macroscopic calcium cycling, advancing the development of more physiologically realistic models of cardiac excitation–contraction coupling.