Multi-Scale Computational Modeling of Spatial Calcium Handling From Nanodomain to Whole-Heart: Overview and Perspectives

Regulation of intracellular calcium is a critical component of cardiac electrophysiology and excitation-contraction coupling. The calcium spark, the fundamental element of the intracellular calcium transient, is initiated in specialized nanodomains which co-locate the ryanodine receptors and L-type calcium channels. However, calcium homeostasis is ultimately regulated at the cellular scale, by the interaction of spatially separated but diffusively coupled nanodomains with other sub-cellular and surface-membrane calcium transport channels with strong non-linear interactions; and cardiac electrophysiology and arrhythmia mechanisms are ultimately tissue-scale phenomena, regulated by the interaction of a heterogeneous population of coupled myocytes. Recent advances in imaging modalities and image-analysis are enabling the super-resolution reconstruction of the structures responsible for regulating calcium homeostasis, including the internal structure of nanodomains themselves. Extrapolating functional and imaging data from the nanodomain to the whole-heart is non-trivial, yet essential for translational insight into disease mechanisms. Computational modeling has important roles to play in relating structural and functional data at the sub-cellular scale and translating data across the scales. This review covers recent methodological advances that enable image-based modeling of the single nanodomain and whole cardiomyocyte, as well as the development of multi-scale simulation approaches to integrate data from nanometer to whole-heart. Firstly, methods to overcome the computational challenges of simulating spatial calcium dynamics in the nanodomain are discussed, including image-based modeling at this scale. Then, recent whole-cell models, capable of capturing a range of different structures (such as the T-system and mitochondria) and cellular heterogeneity/variability are discussed at two different levels of discretization. Novel methods to integrate the models and data across the scales and simulate stochastic dynamics in tissue-scale models are then discussed, enabling elucidation of the mechanisms by which nanodomain remodeling underlies arrhythmia and contractile dysfunction. Perspectives on model differences and future directions are provided throughout.

Introduction

Intracellular calcium (Ca²⁺) handling is a critical component of cardiac electrophysiology (Cheng et al., 1993; Bers, 2002; Song et al., 2015): it governs excitation-contraction coupling (ECC), is involved in multiple signaling pathways, and its impairment has been causally linked to both mechanical and electrical dysfunction of the heart (Eisner et al., 2009; Voigt et al., 2014; Clarke et al., 2015). Elucidating the fundamental mechanisms of Ca²⁺ homeostasis and the perturbations of the system in disease is therefore vital for understanding the electrophysiology of the heart and identifying better diagnostic and treatment strategies for multiple cardiovascular diseases.

One major challenge of dissecting the specific roles and contributions of the many components of intracellular Ca²⁺ handling to observed (dys)function is the complex, non-linear and multi-scale properties of the system in space and time. Spatially, Ca²⁺ sparks, the fundamental element of Ca²⁺-induced-Ca²⁺-release (CICR; see next sub-section), are controlled at the nanometer-scale in localized nanodomains referred to as dyads or couplons, yet Ca²⁺ homeostasis occurs inherently at the cellular-scale where flux balance through the membrane and subcellular transporters determines the total Ca²⁺ levels in the cell and in the Sarcoplasmic Reticulum (SR, the intracellular Ca²⁺ store). Homeostatic conditions of cardiomyocytes at the cell-level are integrated in the heterogeneous syncytium of cardiac muscle tissue where individual myocytes do not function as isolated entities. Temporally, the gating of Ca²⁺ channels occurs on sub-millisecond time-frames (Zahradníková et al., 1999); the heartbeat itself occurs on the order of a second; signaling and regulation, such as sympathetic stimulation, can occur over minutes (Heijman et al., 2011); and transcription and circadian rhythms can influence dynamics over hours or even days (Black et al., 2019; D’Souza et al., 2021). Thus, structure-function relationships from the nanometer- to the whole-heart-scales and dynamics occurring over nanoseconds to hours all contribute to the macroscopic behavior of the heartbeat.

Recent advances in experimental imaging modalities and image-analysis are enabling the super-resolution reconstruction of the structures responsible for regulating Ca²⁺ homeostasis at the nanometer scale (Baddeley et al., 2009; Crossman et al., 2011; Macquaide et al., 2015; Jayasinghe et al., 2018b; Sheard et al., 2019). Extrapolating functional and imaging data from the dyad to the whole-heart is non-trivial due to multi-scale systems interactions; it is therefore a substantial challenge using experimental techniques alone to employ integrative approaches which aim to understand how macroscopic cardiac function arises from these fundamental building blocks. Computational modeling therefore has important roles to play in helping to dissect these structure-function relationships at multiple scales and elucidate the mechanisms by which cellular phenomena translate to the whole-heart.

Over the last decade in particular there have been substantial advances in the complexity and sophistication of computational models of spatial intracellular Ca²⁺ handling. Due to the variety of independently developed models and range of contexts in which they have been applied, it can be a challenging field to get into and understand, whether one is a computational modeler wanting to use and develop these models, or an experimental researcher hoping to understand the models’ limitations and where they can be used to support one’s research. This review aims to provide an accessible entry-point for those not already familiar with these models and a useful reference for those who are. We focus on methods and approaches, in particular those for image-based and multi-scale modeling, how these differ between models, and the implications of these model differences. Applications of the models will be discussed primarily within this context; the reader is referred to previous reviews for more extensive descriptions of the role of computational modeling in elucidating the Ca²⁺-mediated mechanisms of cardiac (dys)function (Heijman et al., 2016; Maleckar et al., 2017; Vagos et al., 2018; Sutanto et al., 2020).

Structure–Function Relationships in Ca²⁺ Homeostasis: Local Control of Ca²⁺-Induced-Ca²⁺-Release

Excitation-contraction coupling is mediated by CICR (Cheng et al., 1993; Bers, 2002), illustrated in Figure 1: (1) Ca²⁺ enters the cell through the L-type Ca²⁺ channels (LTCC) during electrical excitation (the action potential, AP); (2) This local rise in Ca²⁺ activates the ryanodine receptors (RyRs) to trigger a large release of Ca²⁺ (triggered Ca²⁺ spark) from the SR; (3) Ca²⁺ diffuses throughout the myocyte, binds with the contractile apparatus, and initiates cellular contraction; (4) Peak contraction occurs when Ca²⁺ has diffused sufficiently throughout the cell to permit substantial binding with the contractile apparatus; Ca²⁺ influx has largely terminated at this point; intracellular Ca²⁺ is removed into the extracellular space through the sodium-Ca²⁺ exchanger (NCX) and the plasmalemmal Ca²⁺ pump, and SR-Ca²⁺ is restored through the SR-Ca²⁺ pump (SERCA); (5) As NCX and SERCA reduce the Ca²⁺ concentration in the intracellular volume, myofilaments release Ca²⁺ from their binding sites and cellular relaxation occurs; (6) NCX and SERCA continue to act to restore resting Ca²⁺ levels, ready for the next cycle.

FIGURE 1

Figure 1. Simplified illustration of the intracellular Ca²⁺-handling cycle in ventricular myocytes. Left: illustration of an action potential (voltage), calcium transient (CaT), and developed tension, with the different stages broadly labeled on the CaT. Right: cartoon schematic of the different stages of Ca²⁺ cycling in the intracellular volume. High Ca²⁺ concentration is indicated by the green areas; SL refers to the sarcolemmal membrane and TT to a transverse tubule, illustrated on a real cell image (modified from Gadeberg et al., 2016). Numbered stages correspond to those in the main text and are illustrated at different spatial locations of the schematic for brevity – it should be clarified that all stages occur throughout the cell volume.

The rise and decay of Ca²⁺ in the intracellular volume is referred to as the intracellular Ca²⁺ transient (CaT), and broadly follows the AP (Figure 1). Rather than being a whole-cell, homogeneous event, the CaT is the summation of many thousands of locally controlled Ca²⁺ sparks, conferred in the specialized nanodomains called dyads or couplons which co-localize the RyRs on the SR membrane with the LTCCs on the sarcolemma membrane. Cardiomyocytes therefore feature an intracellular structure which facilitates whole-cell contraction mediated by this local control of CICR: the SR forms a cell-wide network coupling the spatially distributed dyads throughout the intracellular volume; the surface sarcolemma (SL) membrane contains multiple invaginations into the cell interior, consisting of the transverse-tubule (TT) and axial-tubule (AT) system (T-system), harboring LTCCs, NCX, and other ion channels (Dibb et al., 2022) which enables dyads to be formed throughout the cellular volume. The reader is referred to the previous works and reviews by Cannell and Kong (2012, 2017) and Laver et al. (2013) for discussion about the importance of local control to explain the properties of CICR in the heart.

Although the Ca²⁺-handling system is conceptually similar and involves the same machinery in all regions of the heart, there are functional and structural differences between myocytes from the pacemaker regions, cardiac conduction system, atria and ventricles that are important for normal physiology as well as the genesis of cardiac arrhythmias (Sutanto et al., 2020). For example, atrial myocytes do not have as robust and dense a T-system as ventricular myocytes (Richards et al., 2011), featuring more orphaned RyR clusters (those without associated LTCCs); in the pacemaker cells of the sinoatrial and atrioventricular nodes, the Ca²⁺ handling system forms the Ca²⁺-clock which is involved in the generation of APs and does not function primarily to initiate cellular contraction (Maltsev and Lakatta, 2013; Yaniv et al., 2015; Maltsev et al., 2017). Discussing models of pacemaker and conduction system myocytes is beyond the scope of the current review, which will focus on the working myocardium of the ventricles and atria.

Multi-Scale Dynamics of Ca²⁺ Handling

The inherently multi-scale nature of cardiac Ca²⁺-handling is perhaps best illustrated by considering the mechanisms of spontaneous Ca²⁺ release events (SCRE) and their involvement in proarrhythmogenic premature focal excitations (Figure 2).

FIGURE 2

Figure 2. Illustration of the scales involved in cardiac Ca²⁺-handling electrophysiology. Left, upper: illustration of a single RyR and the distribution of RyRs in a dyad, RyR channel function at low and high local Ca²⁺, and propagation of a Ca²⁺-spark between dyads. Cell-structure image re-rendered using data from Jayasinghe et al. (2018a). Left, middle: Illustration of an intracellular Ca²⁺-wave and the impact of the spontaneous Ca²⁺-transient on NCX and the action potential (AP), demonstrating both a delayed-after-depolarization (DAD) and triggered AP (TA). Left, lower: illustration of a focal excitation in tissue which degenerates into a re-entrant excitation pattern. Simulation data from Colman et al. (2013). Right, upper: spatial model of a single dyad (model as presented in Sheard et al., 2019) and illustration of a sub-micron full-cell model with dyads distributed therein. Right, middle: illustration of a coarse-grained 3D cell model with different intracellular and SR compartments labeled. Right, lower: illustration of standard, non-spatial Hodgkin-Huxley electric circuit models of cardiac electrophysiology, their coupling in a 1D strand, and implementation in 3D whole-heart models based on imaging data, reproduced using data from Benson et al. (2011).

One consequence of local control is that each dyad contains only small numbers of channels, with typical values thought to be 5–15 LTCCs and 5–200 RyRs (Baddeley et al., 2009; Jayasinghe et al., 2018a,b), located within a very small volume of order < 10^–3 μm³ (Scriven et al., 2013). Stochastic oscillations of single RyR channels can result in small-scale intracellular Ca²⁺-release (Ca²⁺ quarks), potentially inducing the nanodomain-wide event of a spontaneous Ca²⁺ spark by recruiting further RyRs within the dyad to sustain a release flux. The specific mechanism for this recruitment is either a large flux of Ca²⁺ through the pore of a channel that stochastically opens, which raises the local Ca²⁺ sufficiently to trigger the opening of adjacent RyRs, or a cascade of events where a single channel opening increases the probability of opening of a second round of receptors, which triggers more rounds of receptors opening (Asfaw et al., 2013). Irrespective of the specific mechanism, once a reasonably large spark is originated, spatial-diffuse coupling provides a substrate for the propagation of Ca²⁺-sparks throughout the cell as a spark-induced-spark mediated Ca²⁺-wave (Figure 2).

These SCRE are potentially pro-arrhythmic cellular phenomena: Ca²⁺ release can activate NCX which results in a transient inward current when Ca²⁺ is extruded, depolarizing the cell membrane potential as a delayed-after-depolarization (DAD) or, if of sufficient magnitude, a full triggered AP (TA; Figure 2). Multiple cells must undergo some degree of synchronization of these TA in order for them to overcome electrotonic load and manifest in tissue as a focal excitation (Xie et al., 2010; Campos et al., 2015; Liu et al., 2015; Colman, 2019). Similar considerations apply for many sub-cellular Ca²⁺ handling phenomena, from rate-dependence to arrhythmogenic CaT alternans; the fundamental pumping function itself ultimately depends on these multi-scale interactions and can thusly be potentially perturbed by random, stochastic oscillations at the nanometer-scale. Ca²⁺-dependent regulation of the membrane potential is one of the key factors in understanding arrhythmogenesis. Elucidating these mechanisms is, alongside ECC and contractile function, a primary motivation for the development of biophysically detailed models of intracellular Ca²⁺ handling. For a comprehensive overview of the multi-scale implications of Ca²⁺ handling in normal and abnormal cardiac function, the reader is referred to, for example, the reviews of Eisner et al. (2009, 2017).

Multi-Scale Computational Models of Spatial Ca²⁺ Handling

Due to the importance of local control, common-pool models of the cell – that is, those which treat the intracellular space as single homogenized volumes – fail to properly capture the underlying mechanisms and dynamics of Ca²⁺ handling. For example, Sato et al. (2013) demonstrated that stochasticity in Ca²⁺ cycling is necessary to explain the emergence of discordant alternans, which cannot be reproduced with deterministic, common-pool cell models. Models which explicitly account for the spatial nature of the cardiomyocyte as well as stochastic dynamics of the RyRs and LTCCs are therefore much better suited to detailed analyses of Ca²⁺-handling phenomena. However, these models are also computationally more intensive than common-pool models, by a factor of >10⁵, and less suitable for tissue-scale and especially whole-heart simulations. Therefore, different models need to be considered at different spatial scales (Figure 2).

This review will discuss models of spatial Ca²⁺ handling at the multiple scales of the single nanodomain, the cardiomyocyte, and the whole-heart. Models describing the kinetics of the RyRs will first be discussed, followed by spatial models of the single nanodomain. Approaches to whole-cell modeling will then be discussed, with a particular focus on mechanisms of spatial Ca²⁺ coupling. This discussion will then be expanded to approaches for modeling variable and heterogeneous sub-cellular structure and the integration of experimental imaging data. Finally, approaches to develop simplified, computationally efficient models which still capture important features of spatial and stochastic Ca²⁺ handling will be discussed, in both the context of providing generalizable mechanistic explanations and for performing tissue-scale simulations of many thousands or millions of coupled cells. Overall clarity is prioritized over providing substantial details of all available models and investigations; thus, this review should not be considered exhaustive. The reader is also referred to the extensive overview of multi-scale mathematical and computational modeling methods presented in Qu et al. (2014).

Modeling Ca²⁺-Induced-Ca²⁺-Release: Descriptions of Ryanodine Receptors’ Kinetics

There are multiple descriptions of RyR kinetics which have been used in computational modeling. The simplest form are two-state models, which have only a closed/inactivated state and an open/activated state. Most frequently used are four-state models, which have more details of refractoriness and inactivation (which may or may not physiologically occur, discussed in the later sub-section “Perspectives on model differences”) and can simulate different potential mechanisms of these behaviors. The model structure of the RyR is directly related to the type of behavior that is considered to be behind the appearance of local sparks. In two-state models, refractoriness is not considered relevant in the dynamics of initiation and termination of sparks, whereas refractoriness and more complex gating, that may play critically important roles in the regulation of spark dynamics, can be included in four-state and other models. Previous studies and reviews have compared fundamentally different models of the RyR regarding their ability to reproduce different features of the physiology and/or assess the ability of different proposed mechanisms to explain these features. For example, Stern et al. (1999) evaluated different Markov-chain model constructions for reproducing CICR, and Cannell and Kong (2012, 2017) assessed different mechanisms of Ca²⁺-spark termination.

This section aims to provide a clear indication of the models used in various studies by different groups and explain what the major features of these models’ differences are, in order to guide the reader through the many studies. The focus is primarily on those models which are utilized in whole-cell simulations, rather than those designed specifically to evaluate RyR function in isolation or in bilayers, such as Zahradníková et al. (1999).

Model Structure

The four-state Markov chain RyR model (Figure 3) which forms the basis for the majority of modeling studies, originally presented in Stern et al. (1999), is governed by the following state-equations:

$\frac{dC}{dt}=O\cdot K_{O-C}+C^{{}^{*}}\cdot K_{C^{{}^{*}}-C}-C\cdot (K_{C-O}+K_{C-C^{{}^{*}}})(1)$

$\frac{dO}{dt}=C\cdot K_{C-O}+O^{{}^{*}}\cdot K_{O^{{}^{*}}-O}-O\cdot (K_{O-C}+K_{O-O^{{}^{*}}})(2)$

$\frac{d{C}^{{}^{*}}}{dt}=O^{{}^{*}}\cdot K_{O^{{}^{*}}-C^{{}^{*}}}+C\cdot K_{C-C^{{}^{*}}}-C^{{}^{*}}\cdot (K_{C^{{}^{*}}-O^{{}^{*}}}+K_{C^{{}^{*}}-C})(3)$

FIGURE 3

Figure 3. Illustration of RyR Markov-chain models. The four-state, two-state and three-state models are illustrated with the different features highlighted and a selection of publications in which they are used listed. Block arrows indicate model inheritance.

And by definition/conservation:

$O^{{}^{*}}=1-(C+O+C^{{}^{*}})(4)$

Where O, O*, C, and C* represent the four states of the model and correspond to different physical conditions dependent on the specific model implemented (e.g., active; refractory; inactivated; buffer binding state). The release flux, J_rel, is generally given by:

$J_{rel}=J_{rel}^{\max }\cdot O\cdot ({\left[C{a}^{2+}\right]}_{SR}-{\left[C{a}^{2+}\right]}_i)(5)$

Where J_rel^max is the maximal flux rate and the subscripts SR and i denote the SR and intracellular Ca²⁺ concentrations. The models have the following symmetries (some of which are broken in further model developments):

$K_{C-O}=K_{C^{{}^{*}}-O^{{}^{*}}}(6)$

$K_{O-C}=K_{O^{{}^{*}}-C^{{}^{*}}}(7)$

$K_{C-C^{{}^{*}}}=K_{O-O^{{}^{*}}}(8)$

$K_{C^{{}^{*}}-C}=K_{O^{{}^{*}}-O}(9)$

If all of these symmetries are preserved, the description becomes equivalent to a Hodgkin-Huxley model with two independent gating variables.

The functional form of these transition rates, including the variables that they depend on and which states correspond to release flux, differ between the different implementations. In the original study by Stern et al. (1999), and henceforth referred to as the “Stern-like” models, one single state (O) corresponds to the release flux [equation (5)]; transitions from the closed-to-open states, which are symmetric for C-O and C*–O*, are dependent on the square of local cytosolic Ca²⁺ (“Ca_i” from herein for brevity and to avoid confusion with notation of powers); transitions from the top-to-bottom, also symmetric, are linearly dependent on the local Ca_i:

$K_{C-O}=K_{C^{{}^{*}}-O^{{}^{*}}}=f_0(C{a}_i^2)=k_{0}C{a}_i^2(10)$

$K_{C-C^{{}^{*}}}=K_{O-O^{{}^{*}}}=f_1(C{a}_i)=k_{1}C{a}_i(11)$

$K_{O-C}=K_{O^{{}^{*}}-C^{{}^{*}}}=k_2(12)$

$K_{C^{{}^{*}}-C}=K_{O^{{}^{*}}-O}=k_3(13)$

Where k_0–3 are constants. This model is directly used in, for example, Nivala et al. (2012a,b, 2015).

Shannon et al. (2004) updated this formulation to introduce the SR-Ca²⁺ (Ca_SR for brevity) as a gating variable, which will be referred to here as “direct luminal gating.” The transition rates from the closed-to-open states are now functions of Ca_SR as well as Ca_i², and the transitions from active to inactive (top-to-bottom) are dependent on both Ca_SR and Ca_i:

$K_{C-O}=K_{C^{{}^{*}}-O^{{}^{*}}}=f_4(C{a}_i^2,C{a}_{SR})=\left(\frac{k_4}{k_{CaSR}}\right)C{a}_i^2(14)$

$K_{C-C^{{}^{*}}}=K_{O-O^{{}^{*}}}=f_5(C{a}_i,C{a}_{SR})=k_5{k}_{CaSR}C{a}_i(15)$

Where

$k_{CaSR}=S{R}^{\max }-\left(S{R}^{\max }-S{R}^{\min }\right).{\left[1+{\left(\frac{E{C}_{50}^{SR}}{C{a}_{SR}}\right)}^H\right]}^{-1}(16)$

This formulation is one of the most-commonly used (Alvarez-Lacalle et al., 2013, 2015; Voigt et al., 2014; Marchena and Echebarria, 2018, 2020; Sutanto et al., 2018; Vagos et al., 2020, 2021).

Whereas the above model introduced direct luminal gating, many models include “buffer-mediated luminal gating,” wherein the Ca_SR influences gating of the RyRs not directly but rather mediated through the SR-Ca²⁺ buffer calsequestrin (CASQ), as introduced by Restrepo et al. (2008) and similar to that of Gaur and Rudy (2011). The closed-to-open transition rates have the same form as the original Stern formulation and the top-to-bottom rates are now dependent on CASQ. Two major differences are: (1) There are now two open states (O and O* in the illustration) with the formulation for J_rel [equation (5)] updated accordingly [equation (21)]; (2) The constant for the open transition rate for the lower portion of the model (which now corresponds to CASQ-bound) is smaller than that of the top, breaking one of the symmetries in the previous models [equation (6) is no longer true]:

$K_{C-O}=f_6(C{a}^2)=k_{6}C{a}_i^2(17)$

$K_{C^{{}^{*}}-O^{{}^{*}}}=f_7(C{a}^2)=k_{7}C{a}_i^2(18)$

Where

$k_6>k_7(19)$

And CASQ determines the unbound-bound transition rates:

$K_{C-C^{{}^{*}}}=K_{O-O^{{}^{*}}}=f_8(C{a}_i,CASQ)(20)$

The release flux is now given by:

$J_{rel}=J_{rel}^{\max }\cdot \left(O+O^{{}^{*}}\right).(C{a}_{SR}-C{a}_i)(21)$

The reader is referred to the work of Restrepo et al. (2008) for details of the CASQ buffering and gating equations [equation (20)]. Models which use this formulation include Song et al. (2015, 2017, 2018) and Shiferaw et al. (2017).

In Sato and Bers (2011) and subsequent studies (Sato et al., 2012, 2014, 2016), this model was updated in order to reduce the number of RyRs open during a Ca²⁺ spark:

$K_{C-O}=f_8(C{a}_i^2)=\frac{k_{8}C{a}_i^2}{k_{cp}^2+C{a}_i^2}+w(22)$

$K_{C^{{}^{*}}-O^{{}^{*}}}=f_9(C{a}_i^2)=\frac{k_{9}C{a}_i^2}{k_{cp}^2+C{a}_i^2}+w(23)$

Where k_cp and w are further constants. The majority of the models presented in the field – especially those of whole-cells – implement an RyR model which falls into one of these three broad categories. It is worth noting that the relatively simple functional forms of the transition rates given for the different models above may be modified in studies which aim to fit to experimental data describing RyR open probability, e.g., Voigt et al. (2014), Sutanto et al. (2018), and Vagos et al. (2020); these studies contain further parameters and more complex functions, but their Ca_i and Ca_SR dependence is still captured in the general forms of the model given above. Further alternatives and updates exist of these baseline models. In Colman et al. (2016, 2017a,2017b) and Colman (2019), a functionally motivated hybrid was developed wherein only one state corresponds to the open condition but the luminal dependence is buffer-mediated rather than direct.

One feature of the models, which is independent of the fundamental model structure but nonetheless important for model behavior, is the Hill coefficient, H, to which Ca_i is raised by for the closed-to-open (left-to-right) transition rates. Many of the models use a simple coefficient of H = 2. However, single-channel and single-nanodomain studies (Sobie et al., 2005; Cannell et al., 2013) indicate that H is species-dependent and varies in the range 2–2.8 which may be implemented in some studies.

Alternative Models of the Ryanodine Receptors

There are alternative formulations to describe RyR kinetics which either: (1) are not of the form of a four-state model or (2) introduce further environmental variables to control gating. In Song et al. (2016), a reduced two-state approximation of the four-state RyR model was introduced in a study focused on elucidating long-lasting Ca²⁺ sparks. A reduced, or minimal, two-state model (Figure 3) was used in studies of single nanodomains (Greenstein and Winslow, 2002; Walker et al., 2014; Sheard et al., 2019). In Greene and Shiferaw (2021) a three-state RyR model (Figure 3) was implemented which does not correspond directly to a reduction of the four-state model. The model included a second closed state after the open state, which was introduced to reproduce “flicker” based on Mukherjee et al. (2012), as well as containing regulation of the RyRs by Calmodulin (CaM). CaM was also included in the deterministic model presented in Wei et al. (2021).

There is also the question of whether allosteric interactions/cooperativity play a role in RyR gating, with Marx et al. (2001) demonstrating that the regulatory subunit FK506-binding protein could functionally couple neighboring RyRs to underlie coordinated gating. These interactions are included in many works (Stern et al., 1999; Sobie et al., 2002; Chen et al., 2009; Greene and Shiferaw, 2021) and have been proposed as one explanation for self-termination of the Ca²⁺ spark, as discussed in Cannell and Kong (2017), although as argued in that review, unlikely to be a major contributor to this phenomenon.

Numerical Solutions to the Ryanodine Receptors Model

In common-pool models of the cardiomyocyte, the solutions to the RyR/LTCC models are typically numerically approximated using deterministic algorithms such as the forward-Euler method. These numerical solutions correspond to tracking only the average state of the system, i.e., the proportion of open RyRs/LTCCs across the whole-cell; information on the state of individual channels or channel clusters is not preserved in such an approximation. One motivation for the development of detailed spatial models of the dyad or cardiomyocyte is to capture the stochastic (random) nature of individual RyRs/LTCCs as well as their local control, due to the relevance of both of these features for both CICR and more complex emergent dynamics such as Ca²⁺-waves. Thus, deterministic solutions are no longer suitable. Instead, stochastic algorithms that explicitly account for randomness and track individual channels are required.

The most straight-forward method is to implement the Monte-Carlo approach: the state of each individual channel is tracked directly, and state-transitions are determined based on random numbers and the probability of transition. For example, for a two-state RyR model corresponding to only closed (C) and open (O) states, the algorithm at each time-step (Δt) might look like:

Loop over all RyRs:

RAND = generate random number between 0 and 1

IF state is equal to C:

IF RAND < K_C–O × Δt: state becomes O

ELSE state remains C

ELSE IF state is equal to O:

IF RAND < K_O–C × Δt: state becomes C

ELSE state remains O

The release flux is then given by the sum of the open channels in each dyad. This approach is ideally suited to cases where each individual channel is required to be tracked, for example in spatial models of the single nanodomain (Mesa et al., 2021). However, in larger models, e.g., of the whole-cell, this then requires the state of ∼50 RyRs + ∼15 LTCCs per dyad × ∼20,000 dyads to be tracked individually, with random numbers generated to determine state transitions for each one, which can be computationally intensive. There are more sophisticated approaches that can capture these same dynamics but at a reduced computational cost, such as a modified Gillespie’s algorithm (Gillespie, 1976; Rathinam et al., 2003; Nivala et al., 2012a; Song et al., 2019), or the Fokker-Planck or Langevin equations (Herzel, 1991; Wang et al., 2015), which can be briefly summarized as the addition of noise to a deterministic solution. Further approaches to capture the stochastic nature of RyR dynamics at a reduced computational cost, suitable for large-scale tissue simulations, are described in the final section of this review: “Simplified, minimal and tissue models.”

Perspectives on Model Differences

The disparity between the gating mechanisms of the different models could have important implications for model dynamics and thus mechanistic conclusions drawn from these simulations. Due to inter-model differences in the setup of the whole-cell (see section “Spatial models of the whole cardiomyocyte”), it is not necessarily trivial to directly evaluate RyR function, as dynamics are intricately linked with other model parameters such as dyadic cleft volume, LTCC formulation/magnitude, local Ca²⁺ buffering, and spatial Ca²⁺ coupling; an RyR model often cannot simply be “dropped in” or “swapped out” in a whole-cell model. One major difference in model function is the typical number (or proportion) of RyRs that open in a given cluster during a triggered Ca²⁺ spark, as directly addressed in the reformulation presented in Sato and Bers (2011). A similar model structure (e.g., the same four states and functional dependence on local Ca²⁺ concentrations) but with different parameters governing the transition rates and maximal Ca²⁺-flux rate may lead to very different outcomes in adaptive function regarding Ca²⁺-spark dynamics and homeostasis.

Another difference is the description of RyR inactivation and refractoriness, corresponding to the multiple approaches discussed above (e.g., no inactivation, luminal gating, buffer-mediated gating). One reason for the disparities between the model structures is the current debate regarding the fundamental physiological relevance of RyR inactivation (Cannell and Kong, 2017), where there is no strong evidence that significant inactivation occurs under physiological conditions. However, there are strong indications of interactions between CaM binding sites and opening properties of the RyR that may be relevant for mathematically equivalent states to inactivation: Wei et al. (2021), for example, demonstrate that inactivation may be mediated by CaM-Ca²⁺ binding to the RyR and that this plays a significant role in CaT alternans. An important implication of the model choice to either include or omit significant RyR inactivation is the degree to which the junctional SR depletes during a triggered Ca²⁺-spark: models which do not include RyR inactivation (e.g., Hake et al., 2012; Cannell et al., 2013) exhibit substantially greater depletion of local SR-Ca²⁺ concentrations than those which do include inactivation (e.g., Restrepo et al., 2008; Colman et al., 2017b; Sutanto et al., 2018). This further highlights the challenge of evaluating RyR models under the same environmental cell conditions: the extent of junctional SR depletion will have large implications on homeostasis when combined with the specific formulations and parameterization of the Ca²⁺ buffers, SERCA and NCX, which primarily control the balance of SR-Ca²⁺-refilling and cellular Ca²⁺ efflux. Recent studies have also highlighted the direct importance of SERCA function for Ca²⁺ homeostasis and the dynamics of Ca²⁺ sparks (Hake et al., 2012; Sato et al., 2016, 2021; Holmes et al., 2021).

A further factor, which has recently been included in Berti et al. (2017), is the role of other ions such as K⁺, Mg²⁺, Cl^–, and the counter-ion fluxes they facilitate during SR-Ca²⁺ release which help to maintain the trans-SR membrane driving force, especially at rapid pacing rates. It is possible that these dynamics could also influence RyR refractoriness by modulating this driving force.

Considerations of further complexity, such as the inclusion/omission of regulation by CaM (or other potential RyR and Ca²⁺ signaling modulators) can be directly motivated by the aim of the specific study; it is generally the perspective that additional complexity should only be included where specifically required, in order to reduce the influence of the propagation of unknown errors. However, it could also be argued that a non-linear, multi-scale complex system such as this presents the possibility for unpredictable emergent phenomena, which may depend on the interaction of factors such as CaM with other variables. In this case, it can be argued that one should aim to include as many (rigorously derived) components of the system as is feasible. Perspectives on this are ultimately philosophical and it would be unwise to disregard either argument.

Whereas the original study of Stern et al. (1999) compared multiple formulations of the RyR, this necessarily did not include the more recent updates described above (Shannon et al., 2004; Restrepo et al., 2008; Sato and Bers, 2011; Song et al., 2016; Greene and Shiferaw, 2021). Thus, a comprehensive benchmarking study which determines the implications of these model differences in relation to multiple dynamic Ca²⁺-handling-mediated phenomena would be hugely beneficial. Previous reviews and studies, such as Cannell et al. (2013), have performed this in specific contexts such as in the evaluation of the mechanism of termination of Ca²⁺ sparks, but a more holistic benchmarking study which considers RyR function in the context of multiple relevant factors simultaneously (such as CICR, Ca²⁺ spark termination, spatial Ca²⁺ coupling, spontaneous Ca²⁺ spark dynamics, responses to changes in pacing rate and cell environment, reproduction of alternans or after-depolarization) has yet to be performed. Such a study, requiring the whole-cell models described later in this review, could be hugely valuable in understanding the features of the different models and revealing fundamental insight into the physiology of cardiac Ca²⁺ handling.

Spatial Models of the Single Dyad/Nanodomain

Only recently have experimental imaging techniques been able to resolve individual RyRs (Protasi et al., 1998; Baddeley et al., 2009; Crossman et al., 2011; Macquaide et al., 2015; Jayasinghe et al., 2018b; Sheard et al., 2019), enabling the structure of dyads (including number and arrangement of RyRs) to be described in detail. However, it is still challenging to correlate local Ca²⁺ concentration to RyR activity and quantify the specific fractional opening of RyRs during a typical spark in imaging data, and optical methods are limited to the close proximity of the cell surface. Thus, computational models of the single, spatially distributed dyad are useful to understand the mechanisms of Ca²⁺ sparks and their dependence on RyR number and arrangement.

In the original paper by Stern et al. (1999), which evaluated multiple RyR Markov-chain structures including the four-state model described in detail in the previous section, dynamics were evaluated using a spatial model of the nanodomain. The 2D model was discretized at a resolution of 10 nm, with RyRs being arranged in a regular lattice with spacing 30 nm (i.e., each RyR of size 30 nm occupies a 3 × 3 grid at dx = 10 nm; Figure 4A). Dynamics were evaluated using different numbers of RyRs in the dyad. Whereas Ca²⁺ could diffuse within the dyadic space, the rapid equilibrium approximation was generally implemented in order to improve computational efficiency, a far more pressing constraint in 1999 than presently.

FIGURE 4

Figure 4. Spatial models of the dyad. (A) Illustration of model structure (left) and different approaches to RyR arrangement (right), ranging from regular lattice through irregular lattice to free-placement. (B) Illustration of a triggered Ca²⁺ spark occurring in two different RyR arrangements, corresponding to control and right-ventricular heart failure (RV-HF) directly reconstructed from imaging data. Scale bar corresponds to 50 nm. Reproduced with data from Sheard et al. (2019).

In Louch et al. (2010) a 3D cylinder model of the dyad was employed to assess the causes of dyssynchronous Ca²⁺ release in heart failure. The model was combined with AP recordings from control and heart failure myocytes. They observed that whereas AP prolongation reduces the driving force for Ca²⁺ entry through the LTCCs, this is balanced by the increase in RyR sensitivity which results from steady-state increased SR-Ca²⁺. Thus, they conclude that other factors (such as T-system disruption, see section “Modelling variability in sub-cellular structure and function”) underlie the loss of spatial synchrony in the CaT.

A detailed 3D model of multiple structures surrounding a dyad, reconstructed from electron tomograms of a mouse ventricular myocyte, was presented by Hake et al. (2012). Whereas not focusing on RyR distribution, the model included a substantial level of detail in the structures of the junctional SR (which was split into three parts of the back, release, and rim locations) and local T-tubules. The model implemented a phenomenological description of RyR gating as well as local buffers (including CASQ) and fluxes (including SERCA and NCX), and revealed the substantial gradients in local Ca²⁺ in the cytosol which emphasized the importance of the location of SERCA and its role in facilitating long-lasting Ca²⁺ sparks.

Walker et al. (2014) developed a 3D model of the dyad, junctional SR and TT based on the previous work of Sobie et al. (2002) and Williams et al. (2011). Whereas Ca²⁺ could diffuse in 3D, the RyRs were arranged on a 2D pancake (Figure 4A). Similar to the Stern et al. (1999) model, the spatial resolution was ∼10 nm with minimum RyR-spacing of ∼30 nm; RyRs were arranged on this regular lattice but now in irregular patterns. This study demonstrated the importance of specific RyR arrangement – not just total number – on Ca²⁺ spark dynamics and therefore highlighted the importance of Ca²⁺ diffusion within the dyad.

A similar model was also presented in Cannell et al. (2013) and Laver et al. (2013), which included a network SR with SERCA and its associated intracellular uptake flux. The RyR model was specifically fit to Ca²⁺ concentrations corresponding to CICR and thus not suitable for quiescent/spontaneous spark analysis. These models supported the mechanism of induction-decay for Ca²⁺-spark termination, described mathematically by Hinch (2004) and referring to the mechanism by which reduced release flux as the SR depletes results in increased closed-times of adjacent RyRs and thus an increase in the probability of spontaneous spark termination, as a direct consequence of the steep relationship between local Ca²⁺ and the closed-to-open transition rate for the RyRs. This work also demonstrated that other more complex mechanisms (such as inactivation of the RyRs, discussed earlier) were not necessary to explain experimental observations, although they do state that they likely still play a role facilitating multiple points of control.

In Sheard et al. (2019) a 2D model was presented based directly on super-resolution imaging data. The model was also discretized at a resolution of 10 nm with minimum RyR-spacing of 30 nm. In this case, RyR positions were not arranged in a lattice, but were freely placed, directly determined by experimental images of multiple dyads in both control and right-ventricular heart failure (RV-HF) conditions (Figure 4). This model also included differentiation between non-phosphorylated and phosphorylated RyRs (modeled as a simple increase in sensitivity to Ca²⁺), and results support the conclusions of both Cannell et al. (2013) and Walker et al. (2014) regarding the importance of both number and specific arrangement of RyRs for triggered Ca²⁺ spark dynamics and fidelity. Mesa et al. (2021) also investigated the functional impact of selected phosphorylated RyRs within a cluster, indicating that they can play a compensatory role in recovering healthy spark dynamics which had been lost through cluster disruption. Iaparov et al. (2021) developed a 2D model which included further possibilities for RyR arrangement and spacing (including both approaches outlined above), number, and dyad extent. They found that Ca²⁺ spark occurrence varied with the spatial arrangement, but did not consistently correlate with total RyR number, the magnitude of Ca²⁺ current or the surface density. This model included allosteric interactions and the authors found RyR coupling strength to be a major factor underlying sparks.

Challenges and Future Directions; Importance of Spatial Ca²⁺ Coupling

As the experimental structural data improves in resolution and quality it will become more feasible and important to develop models based directly on these data, requiring semi-automated image processing pipelines to produce these geometries in both 2D and 3D. This will enable the models to be applied to more specific conditions including disease states, with many pathologies producing complex and heterogeneous subcellular remodeling. However, there are still numerous challenges in regards to model validation of Ca²⁺ spark spatio-temporal functional properties. Because of the difficulty in measuring Ca²⁺ function and underlying sub-cellular structure simultaneously, computational modeling can be useful to help fill in the gaps to relate structure to function, but for the same reasons, challenging to directly validate. It is worth highlighting that the models of single nanodomains do not reproduce well experimentally measured values for the full width half maximum (FWHM) of Ca²⁺, generally resulting in values of ∼1 μm which are below the 1.8–2.2 μm in experiment. Hoang-Trong et al. (2021) for example did simulate a realistic feature of 1.85 μm, but this required using a large RyR cluster combined with two smaller satellite clusters. However, due to spatio-temporal limits on the resolution of functional imaging experiments, sparks smaller than given sizes are not detected experimentally with any given accuracy. This generates an arbitrary experiment-dependent cut-off that affects the spark distribution and thus makes comparisons between simulation and experiment non-trivial.

A potential limitation of these models is the high spatial-resolution (and thus small voxel/element volumes) required to model RyR distribution: whereas this does not pose a problem during CICR, in which local Ca²⁺ concentration is relatively large, it does pose a problem at resting/quiescent Ca²⁺ levels, where the low Ca²⁺concentration of ∼0.1 μM in small volumes of 10^–21–10^–18L (Scriven et al., 2013) corresponds to the presence of countable numbers of Ca²⁺ ions (i.e., ion distribution is discrete and the notion of a well-defined concentration is debatable); similar considerations may also apply in the junctional SR during CICR when its Ca²⁺ load has been depleted (Hake et al., 2012; Cannell et al., 2013). Continuous approximations may therefore no longer be valid, and simulations of spontaneous Ca²⁺-sparks, in particular, may be non-trivial to implement and analyze. Hybrid schemes which implement spatial stochastic methods to capture the trajectories of individual particles, such as presented in skeletal muscle simulations implementing the Mcell framework (Kerr et al., 2008; Holash and MacIntosh, 2019), may offer solutions to this challenge.

These spatial nanodomain models have demonstrated the importance of specific RyR arrangement – and not just the total number – in determining Ca²⁺ spark dynamics. This presents a challenge for the translation to whole-cell modeling, as it is not feasible to simulate ∼20,000 spatially distributed dyads in a whole cell at resolutions of ∼10 nm, even if this high resolution is only adaptively applied in the local dyadic space. Therefore, in order to understand how a heterogeneous system of dyads with different RyR numbers and spatial arrangements coordinates in a whole-cell, coarse-graining methods will need to be developed which capture the features of this dyad heterogeneity at a reduced computational cost.

Spatial Models of the Whole Cardiomyocyte

Spatial models of the cardiomyocyte describe the cell as a 2D area / 3D volume throughout which dyads are distributed and within which Ca²⁺ can diffuse in both the cytosolic and SR spaces. It is useful now to bring in the terminology of a Ca²⁺-release-unit (CRU). Whereas some studies use this to refer to the dyad or individual RyRs, cellular-scale modeling studies commonly refer to all of the intracellular and SR compartments associated with a single dyad as a CRU, i.e., it can be thought of as the entire volume of cell surrounding each dyad. Typically, this will contain the network SR (NSR), bulk cytoplasmic space (cyto), and the restricted volume of the dyad (treated as a common-pool). It is also common to include a distinct junctional SR compartment (JSR), as well as other optional sub-spaces, such as the sub-sarcolemma volume just below the surface or T-system membrane, included in order to preserve the higher local Ca²⁺ concentration close to a dyad in regions where the membrane fluxes are located. A major focus of this section is to discuss the various approaches to modeling sub-cellular structure in regards to spatial Ca²⁺ coupling.

Fundamental Model Equations of Ca²⁺ Transport and Homeostasis

Calcium homeostasis in the compartments is described by:

$\frac{d{[C{a}^{2+}]}_{cyto}}{dt}={\mathrm{\beta }}_{cyto}(\mathbf{\text{D}}{\nabla }^2{[C{a}^{2+}]}_{cyto}+{\mathrm{\varphi }}_{cyto}+\left(v_{ss}\mathrm{/}{v}_{cyto}\right)J_{ss})(24)$

$\frac{d{[C{a}^{2+}]}_{SS}}{dt}={\mathrm{\beta }}_{SS}(\mathbf{\text{D}}{\nabla }^2{[C{a}^{2+}]}_{SS}+{\mathrm{\varphi }}_{SS}-J_{ss}+\left(v_{ds}\mathrm{/}{v}_{ss}\right)J_{ds})(25)$

$\frac{d{[C{a}^{2+}]}_{nSR}}{dt}={\mathrm{\beta }}_{nSR}(\mathbf{\text{D}}{\nabla }^2{[C{a}^{2+}]}_{nSR}+{\mathrm{\varphi }}_{nSR}-\left(v_{jsr}\mathrm{/}{v}_{nsr}\right)J_{jSR})(26)$

$\frac{d{\left[C{a}^{2+}\right]}_{ds}}{dt}=\mathbf{\text{D}}{\nabla }^2{[C{a}^{2+}]}_{ds}+{\mathrm{\varphi }}_{ds}-J_{ds}(27)$

$\frac{d{\left[C{a}^{2+}\right]}_{JSR}}{dt}={\mathrm{\beta }}_{jSR}\left({\mathrm{\varphi }}_{JSR}+J_{jSR}\right)(28)$

Transfer between compartments is given by:

$J_{ss}=\left({\left[C{a}^{2+}\right]}_{SS}-{\left[C{a}^{2+}\right]}_{cyto}\right){\mathrm{\tau }}_{ss}^{-1}(29)$

$J_{ds}=\left({\left[C{a}^{2+}\right]}_{ds}-{\left[C{a}^{2+}\right]}_{SS}\right){\mathrm{\tau }}_{ds}^{-1}(30)$

$J_{jSR}=\left({\left[C{a}^{2+}\right]}_{nSR}-{\left[C{a}^{2+}\right]}_{jSR}\right){\mathrm{\tau }}_{jSR}^{-1}(31)$

And the general form for the reaction terms are:

${\mathrm{\varphi }}_{cyto}=J_{NaCa}+J_{pCa}+J_{Cab}-\left(J_{up}-J_{leak}\right)-J_{trpn}(32)$

${\mathrm{\varphi }}_{nSR}=\left(J_{up}-J_{leak}\right)\left(v_i\mathrm{/}{v}_{nsr}\right)(33)$

${\mathrm{\varphi }}_{ss}=J_{NaCa\mathrm{\_}SS}+J_{pCa\mathrm{\_}SS}+J_{Cab\mathrm{\_}SS}(34)$

${\mathrm{\varphi }}_{ds}=J_{rel}+J_{CaL}(35)$

${\mathrm{\varphi }}_{JSR}=-J_{rel}\left(v_{ds}\mathrm{/}{v}_{jSR}\right)(36)$

Where cyto, SS, nSR, ds, and JSR refer to the Ca²⁺ concentrations in each of the (sometimes optional) compartments, β refers to the instantaneous buffering term, ϕ refers to a general reaction term in each compartment, J_x refers to transfer flux between compartments, ∇² is the spatial Laplacian operator in 2D or 3D, describing coupling between CRUs, D is the diffusion constant, v refers to the volumes of the compartments, and τ to the time constants of diffusion. The concentration in the dyadic space can be described by a quasi-steady-state approximation, motivated by the rapid equilibration of Ca²⁺ in this small volume. By setting:

$\frac{d{\left[C{a}^{2+}\right]}_{ds}}{dt}=0(37)$

An approximation for equation (27) can be obtained as in Hinch (2004):

${\left[C{a}^{2+}\right]}_{ds}={\left[C{a}^{2+}\right]}_{SS}+\frac{{\mathrm{\tau }}_{ds}.(k_{rel}.{\left[C{a}^{2+}\right]}_{jSR}+J_{CaL})}{\left(1+{\mathrm{\tau }}_{ds}.k_{rel}\right)}(38)$

Where k_rel is defined by J_rel = k_rel(Ca_SR-Ca_ds) and therefore corresponds to:

$k_{rel}=n_{RyR\mathrm{\_}O}.g_{RyR}.v_{ds}^{-1}(39)$

Where g_RyR is the conductance of a single RyR channel and n_{RyR_O} is the number of open RyRs in the dyad (corresponding to states O or O+O^*, dependant on the RyR model implemented). This approximation enables less constraint on the time-step for the simulation, allowing faster simulations to be performed. A limitation of this approximation is that the introduction of this type of equation leads to a lack of ionic Ca²⁺ conservation in the models; similar issues arise with the implementation of the rapid buffering approximations [β terms in equations (24) – (28)]. Models that try to analyze homeostatic properties require that the computing algorithm conserves ions at all orders (Conesa et al., 2020).

Ca²⁺-Voltage Coupling; Incorporation With an Action Potential Model

In general, it is only the Ca²⁺ concentrations and Ca²⁺-handling channels which are described spatially, with the membrane voltage (V_m) and other ion-currents/concentrations assumed to be homogeneous throughout the cell. This is justified by the fast diffusion of V_m along the cell membrane, indicating that it is valid to assume all channels “see” the same global voltage, at least, in the time-frames considered by the models; moreover, the Debye length in cardiac cells is approximately 1 nm (Mori et al., 2008). Detailed electro-diffusion models, such as implemented in simulations of neurons (Pods et al., 2013), are therefore not generally used or required for cardiomyocytes. The spatial description of the Ca²⁺ handling system can simply replace the equivalent components of common-pool models and can therefore be integrated with either simplified descriptions of the AP (such as assuming it follows a simple, analytical form) or biophysically detailed models of the primary ion currents and global ionic concentrations.

The Ca²⁺ and voltage systems are coupled through the influence of V_m on the activity of the LTCCs and NCX (and any other voltage-dependent Ca²⁺ channel), and feedback of Ca²⁺ into the voltage is captured if a biophysically detailed model of the AP is included, wherein I_CaL and I_NaCa directly influence V_m. Thus, the interaction between global voltage and local Ca²⁺ dynamics can be described in these models, enabling study of the mechanisms of, for example, Ca²⁺-induced AP duration (APD) alternans and afterdepolarizations (Eisner et al., 2009, 2017; Qu et al., 2014).

Whereas explicit diffusion of V_m along the membrane of a single-cell has generally not been included in the described models, this has been simulated in models presented by Crocini et al. (2014) and Scardigli et al. (2018). These studies simulated the impact of disruption of the T-system in disease on the ability of the AP to propagate along the T-system into the interior of the cell, indicating that AP propagation failure can directly contribute to a loss of intracellular Ca²⁺ synchronization.

Model Discretization

Numerical solutions to the spatio-temporal reaction-diffusion equations above require the cell to be described in a discretized space. Models can be broadly categorized as being one of two approaches (Figure 5A):

(1) “CRU-grid” or “compartmentalized” models, wherein the spatial resolution is ∼ 1 μm × 1 μm × 2 μm and each pixel/voxel corresponds to a single CRU;

(2) “Sub-micron” or “free-diffusion” models, wherein the volumes of the SR and cytoplasm, and SS if present, are discretized within each CRU.

FIGURE 5

Figure 5. Discretization and CRU coupling of whole-cell models of spatial Ca²⁺ handling. (A) Illustration of the CRU-grid and sub-micron approaches to model discretisation. Top panels illustrate the whole-cell whereas bottom panels demonstrate the discretisation of the cell into either CRUs or sub-CRU voxels. (B) Illustration of different model structures for inter-CRU coupling. The primary pathway of CRU coupling in the intracellular space is shown by the colored-block arrows. Note that within the cell, coupling is symmetric in the ±-axes directions; a single direction only is shown for the purpose of clarity. Cell image for illustration provided by Dr. Izzy Jayasinghe, The University of Sheffield.

Details of each approach are described below. In general, the CRU-grid approach is more computationally efficient and consequently suitable for statistical simulations and high-throughput, population-cohort models of spatial structure and heterogeneity. Its simplified structure also facilitates the dissection of the mechanisms by which various components contribute to macroscopic function. The advantages of the sub-micron approach are that one can more precisely control aspects such as heterogeneous inter-dyad spacing and co-localization of different channels, they have a more accurate recapitulation of Ca²⁺ diffusion and cellular geometry, and are better suited to direct incorporation of experimental imaging data.

Ca²⁺-Release-Unit-Grid Models

In these models, each voxel contains all of the compartments contained within a CRU. Spatial coupling is described using the isotropic finite-difference method (FDM), or approximations thereof. Each CRU is coupled to its four or six nearest-neighbors (for 2D or 3D, respectively) along the principal axes:

$\begin{array}{c}\mathbf{\text{D}}{\nabla }^2{[C{a}^{2+}]}_{c=cyto,nSR,(SS)}\hfill \\ \approx \frac{\mathrm{\beta }D}{{}^{e_i}d{x}^2}\sum _{i=1}^{i=3}({}^{e_i+1}[C{a}^{2+}]_c{+}^{e_i-1}{[C{a}^{2+}]}_c-2{}^{e_i}[C{a}^{2+}]_c)\hfill \end{array}(40)$

Where e_i refers to the three dimensions (x,y,z), the subscript c refers to any compartment which is spatially coupled (i.e., the bulk cyto and NSR spaces and any coupled sub-space), D is the diffusion coefficient and dx is the spatial step (i.e., the resolution). The diffusion term (D/dx²) is often approximated with a time-constant of diffusion between spatially coupled compartments (τ_c,ei):

$\begin{array}{c}\mathbf{\text{D}}{\nabla }^2{[C{a}^{2+}]}_c\approx J_{Ca\mathrm{\_}diff\mathrm{\_}c}\hfill \\ =\sum _{i=1}^{i=3}\left(\frac{{}^{e_i+1}[C{a}^{2+}]_c{+}^{e_i-1}{[C{a}^{2+}]}_c-2{}^{e_i}[C{a}^{2+}]_c}{{\mathrm{\tau }}_{c,e_i}}\right)\hfill \end{array}(41)$

Note that whereas these approximations are derived from the isotropic FDM, the models are often discretized at a larger resolution in the longitudinal direction than the transverse, reflecting the larger spacing of dyads along the cell compared to along a TT. Thus, dx² in equation (40) or τ in equation (41) can be larger in the longitudinal (z) direction compared to the transverse. This therefore introduces an anisotropy (preferential propagation along the transverse direction) despite the model equations being derived from isotropic approximations.

Sub-Micron Models

In these models, each CRU is further discretized into voxels or elements. Not every voxel necessarily contains all compartments. For example, nanodomains/jSR will only be present in a small subset of voxels, and in some models, the T-system or even nSR may also only be present in a subset of voxels. Spatial coupling is solved using FDM [equation (40)] on regular structured grids or using the more complex finite element method (FEM) on structured or unstructured meshes. They can be discretized at different choices of resolution, generally between 0.05 and 0.2 μm (Nivala et al., 2012a,2015; Colman et al., 2017b; Marchena and Echebarria, 2018, 2020; Hoang-Trong et al., 2021), although there are other intermediary approaches, such as Sutanto et al. (2018) which uses the CRU-grid approach in the transverse direction but is discretized at half-CRU distance in the longitudinal direction.

Model Compartment Structure: Inter-Ca²⁺-Release-Unit Coupling

In addition to model discretization, there are important differences in underlying model structure, notably regarding the mechanisms of inter-CRU coupling. Different structures of inter-CRU coupling lead to variable strengths of spatial Ca²⁺ coupling, in part due to the variable peak Ca²⁺ concentrations in the compartments selected for coupling: compartments which have a smaller volume and are more directly coupled to the dyad exhibit larger CaTs and thus stronger inter-CRU coupling compared to larger compartments which are less directly coupled to the dyad. In this context, coupling strength refers only to the intracellular space, not the SR, the spatial coupling of which is independently controlled. Broadly, there can be considered four different structures which the models follow (Figure 5B and Table 1):

TABLE 1

Table 1. Summary of properties of whole-cell models of spatio-temporal Ca²⁺ handling.

Type 1

The simplest model structure comprises of the four compartments of the bulk cytoplasm and dyadic space and network and junctional SR. The bulk cytoplasm and network SR spaces are the only ones which are spatially coupled throughout the cell. These models contain the weakest coupling due to the coupling of bulk cytoplasm only, which contains Ca²⁺ concentrations of the same order of magnitude as the whole-cell average. Models using this include Nivala et al. (2012a,b, 2015).

Type 2

Many models also include a sub-sarcolemmal subspace from/into which the other membrane fluxes (NCX and the plasmalemmal Ca²⁺ pump) act. This subspace has a smaller volume and higher Ca²⁺ concentrations at peak than the bulk cytoplasm. This sub-space is generally coupled between CRUs in the transverse direction only, i.e., along the TTs. Longitudinal coupling of this sub-space may be present where ATs are modeled, but does not by default occur for every CRU throughout the cell. Coupling strength is higher than the simplest type 1 models due to the larger CaT in this sub-space. Models of this type include Song et al. (2015, 2016, 2018).

Type 3

Other models implement a sub-space which couples CRUs in both transverse and longitudinal directions independently of the presence or absence of T-system/SL. Introduced in Colman et al. (2017a,b) and Colman (2019), this sub-space contains fewer buffers and represents potential pathways between dyads around the intracellular buffers. Given the reduced buffering (higher CaT peaks) and coupling in all directions, these models have stronger inter-CRU coupling than those which contain a sub-sarcolemmal subspace only.

Type 4

Finally, other models have more direct inter-CRU coupling between dyads, or SR-release spaces (Gaur and Rudy, 2011; Voigt et al., 2014; Sutanto et al., 2018; Vagos et al., 2020). Due to this direct spatial coupling of the compartment into which release occurs, these models have the strongest spatial coupling.

This classification is simplified, but captures the major features of the various approaches. For example, the Heijman-lab models (type 4) were originally designed to represent atrial cells with no T-system: interior compartments did not contain the SL fluxes or associated sub-space. However, models of type 1–3 can be generalized to match this structure by removing these same SL fluxes and associated sub-space (where present) from interior CRUs, e.g., as in Colman et al. (2016), Shiferaw et al. (2017), Song et al. (2018), and Marchena and Echebarria (2020). Similarly, the Heijman-lab models can also be generalized to incorporate a T-system by the inclusion of the SL fluxes and sub-space in interior CRUs, e.g., as was performed in Sutanto et al. (2018). Thus, the above types represent four fundamentally different approaches to inter-CRU coupling, and can be generalized to any cell structure, containing, for example, full, no, or variable T-system density (see next section “Modeling variability in sub-cellular structure and function”).

Implications of Model Structures

One major feature of the differences between these models is the relationship between RyR sensitivity, inter-CRU coupling strength and the size of the CaT in normal pacing, which has significant implications for the robustness of inter-CRU Ca²⁺ propagation. Broadly, models with weaker inter-CRU Ca²⁺ coupling (type 1–2 above) tend to contain either a physiologically sized CaT and operate at the threshold of Ca²⁺ propagation, or contain a substantially larger CaT (>2–8 μM) with more robust Ca²⁺ propagation (Nivala et al., 2015; Song et al., 2015; Marchena and Echebarria, 2020). Alternatively, widespread initiation of CICR (reflecting an extensive T-system bringing LTCC close to RyR throughout the cardiomyocyte) contributes to robust, synchronized Ca²⁺ release in some of these models. Models with stronger inter-CRU coupling (type 3–4 above) tend to contain more robust Ca²⁺ propagation while maintaining physiologically sized CaTs (i.e., below 1 μM), as in Voigt et al. (2014), Colman et al. (2017a), Colman (2019), and Vagos et al. (2020). Such robust propagation is particularly relevant when simulating cells without an extensive T-system that rely more heavily on fire-diffuse-fire mechanisms for Ca²⁺-wave propagation, such as in the atria. The implications of these differences are far-reaching: Ca²⁺ propagation is relevant for graded release of CICR, the dynamics of CaT alternans, success or failure of triggered Ca²⁺ wave propagation into regions without T-system, and the dynamics of spontaneous Ca²⁺ sparks and waves.

These features and limitations, generally discussed openly in the original papers, do not detract from the ambitions of the various studies nor question their analyses, as model choices are motivated by the focus of the specific study. It is important to note that all approaches are capable of reproducing all of the phenomena described above, including properties such as the statistics of SCRE and dynamics of alternans, although various model parameters are likely to be substantially different in order to reproduce these same macroscopic features. Nevertheless, it is still important to carefully consider the motivations of the study for model selection. As an oversimplified example, if the ambition is to study Ca²⁺-voltage interactions during SCRE or alternans, it may be best to prioritize CaT magnitude (as this will determine the degree of Ca²⁺-induced inactivation of I_CaL and the magnitude of I_NCX which feedback into the voltage) and use a more “functional” description of inter-CRU coupling (types 3–4) in order to maintain robust Ca²⁺ propagation. However, if the focus of the study is on the mechanisms and implications of inter-CRU coupling then the simpler and (potentially) more physiologically justified inter-CRU coupling structures (types 1–2) may be prioritized. These considerations may be particularly relevant for the development of atrial cell models and those with variable T-system density, and will be discussed in this context in more detail in the next section.

These disparities and compromises indicate that there are fundamental properties of the CaT and inter-CRU coupling which we do not fully understand, and these gaps in our understanding and model differences reflect the large degree of uncertainty in the experimental data on which the models are based. There is an important possibility that different internal structures need different behavior at the dyadic level since weaker coupling generally requires stronger transients to reproduce wave-like propagation. This might lead to important model-dependent conclusions on new research issues. As models become more sophisticated, getting closer to genuine cell- and species- specificity (indeed as more data become available), and are applied in more complex and clinically oriented studies, it will become imperative to solve these issues and develop models which include fully physiologically justified descriptions of RyRs and spatial Ca²⁺ coupling. Further to this, it is these authors’ opinion that future studies would do well to implement multiple, disparate models in order to navigate the limitation of model-dependent conclusions.

Modeling Variability in Sub-Cellular Structure and Function

Cardiomyocytes demonstrate a large degree of inter-cellular, inter-subject and inter-species heterogeneity in properties such as ion-channel expression and sub-cellular structure, and recent studies have highlighted the importance of including such variabilities in, for example, predictive models of pharmacology (Muszkiewicz et al., 2016; Passini et al., 2017). This section will describe how heterogeneity in sub-cellular structure can be captured using models of spatial Ca²⁺ cycling.

T-System Variability: Models of Atrial and Remodeled Ventricular Myocytes

The structure and density of the T-system is one of the most important factors which determines sub-cellular dynamics. In healthy ventricular myocytes the T-system is generally robust and dense throughout the volume of the cell. However, in atrial myocytes and diseased ventricular myocytes the T-system can be substantially sparser and more variable (Lyon et al., 2009; Richards et al., 2011; Gadeberg et al., 2016; Singh et al., 2017). This reduction in T-system density is generally correlated with alterations to the CaT (primarily, a prolonged time-to-peak, reduced spatial synchronization, and often a small, sometimes large, reduction in magnitude) and possibly linked to an increased vulnerability to pro-arrhythmic dynamics (Trafford et al., 2013; Gadeberg et al., 2016; Shiferaw et al., 2017, 2018). In these conditions, triggered Ca²⁺ sparks will occur only in regions of the cell where LTCCs are closely coupled to RyRs [although even these regions may not exhibit triggered sparks if the AP fails to reach them (Crocini et al., 2014)]. Ca²⁺ may then propagate into regions without the T-system, i.e., where orphaned RyRs are found without coupling to the LTCCs, through spark-induced-spark triggered Ca²⁺ waves, resulting in “u” or “w” shaped linescan images (Figure 6). However, Ca²⁺ may also fail to propagate as a triggered wave, leading to regions of the cell which do not undergo substantial Ca²⁺-release which consequently underlies a substantially smaller whole-cell CaT. The conditions which either enable or inhibit triggered Ca²⁺-wave propagation, which may be species-, cell-, disease-, and environment-dependent, are unclear from experimental studies alone, and these analyses form the focus of many computational studies which implement variable T-system density and structure.

FIGURE 6

Figure 6. Models of variable T-system density in atrial and remodeled ventricular myocytes. (A) Experimental data demonstrating “u” waves in pig atrial myocytes, and the corresponding local CaTs at the surface (black) and interior (gray) of the cell. Data is from Gadeberg et al. (2016). (B) Dynamics in different computational models of atrial cells, illustrating recapitulation of “u” waves, dependence on T-system density, and Ca²⁺ silencing and CaT alternans. Data from Colman et al. (2016); Marchena and Echebarria (2020), and Vagos et al. (2021). (C) Experimental and simulation data from remodeled ventricular myocytes (Singh et al., 2017), demonstrating a loss of T-system density and associated loss of spatial synchrony of the CaT, in relation to the Organizational index (OI; experiment) and the proportion of CRUs associated with LTCCs (q; simulation).

Implementation of variable T-system structure is relatively straight-forward in these spatial models of Ca²⁺ handling, for both the CRU-grid and sub-micron approaches. The inclusion or omission of a TT or AT from a CRU or voxel/element can be trivially implemented by either the inclusion or omission of the membrane Ca²⁺ fluxes (LTCCs, NCX, Ca²⁺ pump) and any associated sub-sarcolemma sub-space. Thus, one needs only create a map which describes which CRUs or voxels/elements contain a TT or AT. Creation of this map could be through random selection, a T-system generating algorithm, or directly based on experimental imaging data.

Modeling Atrial Cardiomyocytes

Atrial cells exhibit variable T-system density in both control and disease conditions (Richards et al., 2011; Park et al., 2020), and thus computational models of atrial cells tend to either not include the T-system or explicitly model its variability. Koivumäki et al. (2011) developed a simplified, 1D model of the atrial myocyte which captured the propagation of Ca²⁺ waves from the surface to the interior. The deterministic model contained four sub-cellular compartments as a coarse-grained discretisation of the transversal direction of the cell. Colman et al. (2016) implemented an atrial cell model with variable T-system density in a “type 3,” 3D CRU-grid model. Variable T-system structure was crudely implemented by removing patches of the T-system of variable and controllable size within the cell, based on randomly generated seeds. This model demonstrated “u” and “w” linescans during normal pacing (Figure 6), a delay in the time-to-peak, and a small reduction in the magnitude of the CaT as T-system density decreased. CaT alternans and SCRE were found to increase as T-system density decreased and the specific dynamics were highly controlled by the structure: alternans involved alternating between successful and failed Ca²⁺ propagation into the non-T-system regions, and spontaneous Ca²⁺ waves preferentially emerged from regions without a TT or AT.

The atrial models developed in the Heijman-lab (Voigt et al., 2014; Vagos et al., 2020) contained no T-system by default. These models, of type 4 structure, demonstrate generally robust Ca²⁺ propagation into the cell interior, the mechanisms and sensitivity of which were investigated in detail in Vagos et al. (2020) and further expanded on in Vagos et al. (2021). A mechanism of CaT alternans similar to that described in Colman et al. (2016) was observed (Figure 6). Sutanto et al. (2018) expanded the approach by incorporating RyR/LTCC expression/distribution based on experimental imaging data (see later sub-section “Pipelines for image-based modeling”), as well as the inclusion of variable T-system in the interior of the cell, demonstrating the important role of both TTs and ATs in facilitating Ca²⁺ propagation.

In a series of studies which integrated computational modeling with experimental functional measurements of atrial cell dynamics Shiferaw et al. (2017, 2018, 2020) comprehensively investigated the mechanisms and sensitivity of triggered Ca²⁺ wave propagation using a model of heterogeneous T-system density, in combination with other factors such as β-adrenergic stimulation. In brief, these studies demonstrated rate-dependent triggered Ca²⁺ waves emerging in SR-Ca²⁺ overload conditions. It is argued that the emergence and dynamics of triggered Ca²⁺ waves exhibit a highly non-linear dependence on SR-Ca²⁺ load and that this feature of atrial cells contributes to dynamic instabilities which may be pro-arrhythmic.

Marchena and Echebarria (2018, 2020) also developed models of the atrial cell using a type 1 sub-micron approach. A comprehensive analysis of the relationship between the fractional area occupied by the T-system and the magnitude and morphology of the whole-cell CaT was provided (Figure 6). These studies indicated that SR-Ca²⁺ load was unaffected by T-system density and that reduced CaTs were a consequence of failure of Ca²⁺-wave to fully propagate from the cell periphery into the interior; however, Ca²⁺-release gain was increased in detubulated cells as a consequence of the activation of at least some orphaned RyRs. No significant differences in spontaneous Ca²⁺ spark dynamics were found, indicating that experimental observations of heterogeneous spark dynamics (Hüser et al., 1996; Kirk et al., 2003; Brette et al., 2005) could be explained by different kinetics or regulation of RyRs in membrane and non-membrane regions, not captured in the models. This is consistent with Brandenburg et al. (2016) and Sutanto et al. (2018) which show, in experiment and modeling, that RyR hyperphosphorylation contributed to these regional differences.

Modeling Ventricular Cardiomyocytes With Disease-Related T-System Remodeling

In a combined experimental-computational study, Wagner et al. (2012) characterized the remodeling of sub-cellular structure (T-system properties and other related proteins such as junctophilin) post myocardial-infarction (MI). They observed a progressive, time-dependent post-MI increase in the cross-sectional area of individual T-tubules, a decreased expression of junctophilin, an orphaning of RyR clusters, and uncoupling of CICR. Simulations were performed using the previously presented model of Williams et al. (2011), which implemented the simplest type 1 structure at CRU-grid resolution. Remodeling was incorporated by an increased spacing between TTs and RyRs in a subset of the model’s compartments. Simulation results demonstrated that RyR orphaning contributed to post-MI associated AP prolongation, especially when combined with remodeling of NCX and SERCA, as well as reduced SR-Ca²⁺ load and increased Ca²⁺-leak, potentially contributing toward arrhythmogenic afterdepolarizations. However, the changes to TT cross-sectional area were not captured in this model and require more detailed models, discussed in the later sub-section: “Toward realistic sub-cellular structure.”

Nivala et al. (2015) implemented a loss of the T-system representing remodeled ventricular myocytes using a sub-micron, type 1 model. Voxels containing LTCCs (i.e., dyads) were randomly selected for removal of the membrane components. This study demonstrated that disruption of the T-system led to perturbed spatial Ca²⁺ handling (delay in the time-to-peak and slightly reduced magnitude of the CaT) which was more pronounced when combined with remodeling of whole-cell parameters associated with heart failure.

Song et al. (2018) expanded on the work of Nivala et al. (2015) using a type 2 CRU-grid model. In this study, an algorithm was developed to generate more realistic variable T-system structures. Results were largely concurrent with the previous studies, with the major differences that: (1) a more substantial reduction in CaT amplitude was observed in detubulated cells; (2) a biphasic relationship between T-system density and arrhythmogenic dynamics was observed, wherein intermediate densities demonstrate the most instabilities.

Singh et al. (2017) combined experimental measurements and computational modeling to explore the relationship between T-system density and the features of the CaT in a rat ventricular model of the progression of heart failure (Figure 6C). Experimental measurements, in agreement with their previous study (Shah et al., 2014), demonstrate that the CaT exhibits a slower upstroke and reduced magnitude associated with substantial loss of the T-system, but Ca²⁺ propagation is not silenced. The computational (type 2) model was in strong agreement with these observations, revealing a non-linear relationship on both distribution of release units and separation between LTCCs and RyRs.

Model Parameterization

Differences in the dynamics and sensitivity of triggered wave propagation observed in the models can be largely related to the underlying model structure (type 1–4 as described in the previous section). The development of models of atrial myocytes presents a major challenge in this context, i.e., in obtaining robust triggered Ca²⁺ wave propagation with physiologically sized CaTs but without an over-propensity for spontaneous activity. In the Colman-lab and Heijman-lab models, stronger inter-CRU coupling (type 3–4) was implemented for this purpose, as it facilitated triggered Ca²⁺ wave propagation with physiologically sized CaTs. Without this strong inter-CRU coupling, it was not possible to reproduce stable Ca²⁺ homeostasis during regular pacing while maintaining robust Ca²⁺ wave propagation and physiological CaTs without also observing highly arrhythmogenic spontaneous activity – either the spontaneous spark rate or Ca²⁺-wave rate was too high. Type 3–4 models somewhat solve this through the inclusion of these stronger coupling mechanisms, but also tend to exhibit low spontaneous spark rates and a higher probability for a spark to develop into a full wave. This indicates the possible physiological relevance of stronger coupling (e.g., pathways between buffers and/or along the SR membrane), but this has not been supported with imaging experiments.

It is generally assumed triggered waves should propagate in normal conditions in healthy atrial myocytes, but the validity of this assumption and whether it also translates to remodeled ventricular myocytes is not clear. It is worth noting that different experimental observations cover the range of behavior exhibited by these models, with some experiments observing generally robust triggered Ca²⁺-wave propagation with “u” or “w” waves and only a small reduction in the magnitude of the CaT in cells with low T-system density (Louch et al., 2010; Trafford et al., 2013; Crocini et al., 2014; Gadeberg et al., 2016; Setterberg et al., 2021), whereas others demonstrate more variability and sensitivity of triggered Ca²⁺ waves in many different pacing conditions, associated with substantially smaller CaTs in detubulated cells in normal pacing conditions (Brette et al., 2005; Shiferaw et al., 2017).

Toward Realistic Sub-Cellular Structure

Modeling Cellular Contraction, Mitochondria and Energetics

Two components of the structure-function relationships governing cardiac cellular electrophysiology which have not yet been discussed in detail are the myofilaments (and related contractile apparatus) and the mitochondria (and associated localized buffering and energetics). Okada et al. (2005) developed a FEM-based model of cellular contraction associated with intracellular Ca²⁺ waves, enabling investigation of the impact of cell shortening on Ca²⁺-wave velocity. The model also revealed the potential for spiral Ca²⁺ waves which could maintain arrhythmicity. Hatano et al. (2011, 2012) expanded this model to include realistic local mitochondrial Ca²⁺ buffering and ATP production, T-tubules, SR structure, and myofilaments. These models revealed slow changes in the average mitochondrial Ca²⁺ during the cardiac cycle and that asynchronous contraction caused by a large detubulated region can lead to impairment of myocyte contractile efficiency. Recently, Xie et al. (2018), Song et al. (2019), and Pandey et al. (2021) developed CRU-grid spatial models which also explicitly accounted for local buffering and dynamics of the mitochondria. These mitochondria were assigned to alternating CRUs in the transverse direction (but every CRU in the longitudinal direction) and the models were applied to study the role mitochondria may play in proarrhythmogenic dynamics including afterdepolarizations.

Heterogeneous Channel Distribution

Beyond the T-system, the distribution and local density of different Ca²⁺-handling transporters is also important for governing Ca²⁺ homeostasis through effects on local Ca²⁺-flux balance. Sub-micron models are ideal for controlling the fine details of channel distribution and co-localization between different channels, but CRU-grid models remain suitable for investigating heterogeneous channel/transporter expression in different regions of the myocyte.

The number of RyRs and LTCCs per dyad, as well as dyad volume, are commonly heterogeneous in the default implementations of many models (assigned by scaling the expression or volume by numbers randomly sampled from a normalized Gaussian or other distribution). These heterogeneous properties have been shown to be important for capturing the features of graded intracellular Ca²⁺-release (Greenstein and Winslow, 2002; Shiferaw et al., 2003; Restrepo et al., 2008). In principle, heterogeneous expression of any transporter or component (e.g., SERCA, NCX, buffer concentration) could be implemented in the same manner by sampling scale-factors from a defined statistical distribution. However, determining heterogeneous structure based on experimental data and with constraints on the spatial variation/correlation provide more powerful and physiologically relevant approaches.

Pipelines for Image-Based Modeling

In Sutanto et al. (2018) a pipeline was developed which enabled the expression of RyRs and LTCCs (and in principle any desired Ca²⁺-handling component) observed in experimental imaging studies to be processed to align with the modified CRU-grid model (Figure 7A). Thus, RyR or LTCC expression in each dyad was determined by the intensity of immunofluorescence image labeling for each of these channels. The method involved processing the real cellular data so that it could be registered on the idealized cellular geometry of the computational model.

FIGURE 7

Figure 7. Image-based modeling approaches. (A) Approaches for image-based modeling in idealized, CRU-grid cell models. The top panel illustrates the processing of LTCC/RyR image data so that it can be aligned and registered on an idealized grid computational model, as presented in Sutanto et al. (2018). The lower panel illustrates the method of calculating the correlation length-scale (λ) in order to generate Gaussian-random field maps which match these parameters in order to describe sub-cellular heterogeneity, as presented in Colman et al. (2020). (B) Approaches for directly incorporating imaging data into realistic, sub-micron cell models. The top-left panel demonstrates processing of ultra-structural data describing T-system structure into a mesh for simulations, and illustrates the large Ca²⁺ gradients that can occur under different flux-distributions, as presented in Cheng et al. (2010). Lower panel illustrates processing of data describing the T-system and SR in a portion of a cell, approaches to tesselate the cell portion into a whole-cell, and the substantial Ca²⁺-gradients observed in this model also, as presented in Colman et al. (2017b). Right panel illustrates the semi-automated pipeline for segmenting multiple structures, including the myofibrils and mitochondria, into a 3D structurally-detailed representation of a sarcomere, as presented in Rajagopal et al. (2015).

In an alternative approach, a method was presented in Colman et al. (2020) which involved the development of image-analysis techniques to extract parameters describing the spatial correlation and distribution of the channels (Figure 7A). This involved calculating the length-scale which describes the distance over which expression is correlated. This parameter can be used to generate Gaussian random fields which produce expression maps with the same spatial correlations, but not limited to the specifics of the imaging data.

The above approaches enable efficient and high-throughput simulations of variable cellular structure to be performed, which are ideally suited to statistical analysis and the extraction of fundamental mechanisms of homeostasis. However, sub-micron models, which are substantially more computationally intensive, provide more possibilities for the direct inclusion of experimental imaging data. In a series of studies (Cheng et al., 2010; Kekenes-Huskey et al., 2012; Hake et al., 2014) an ultra-structure model of a realistic single TT and its surrounding half-sarcomeres was developed (Figure 7B), extending their previous work which implemented an idealistic TT (Lu et al., 2009). The local cellular geometry was reconstructed from light-and electron-microscopy images of a rat ventricular myocyte and the studies investigated the impact of the distribution of Ca²⁺-flux channels and buffers. They observed substantial and rapid spatial gradients in local Ca²⁺ concentration in the sub-sarcolemmal sub-space (Figure 7B) and their analyses demonstrate the importance of accounting for T-system ultra-structure and Ca²⁺ flux distribution. These detailed models present the possibility for features of T-tubule ultrastructural remodeling to be included, such as observed in Wagner et al. (2012) and Crocini et al. (2014), and could also potentially be extended to account for other local ionic concentrations and the counter-ion fluxes relevant for SR-Ca²⁺ release, such as presented in Berti et al. (2017).

In Colman et al. (2017b) a sub-micron, whole-cell model was developed which directly incorporated imaging data describing the structure of the T-system, SR and distribution of the dyads (Figure 7B). A semi-automatic pipeline was developed in which the images were processed and down-sampled to create simulation-ready geometries. Whereas only a portion of the cellular geometry was reconstructed, the model was tessellated to construct a whole-cell, exploiting the periodic structure of cardiac myocytes. This model also observed substantial intracellular Ca²⁺ gradients which emerged only in the realistic geometries, in agreement with Cheng et al. (2010). The imaging data was based on electron-microscopy (Pinali et al., 2013) and thus channel distribution was not captured in the data. However, future studies could use advances in correlative light-electron microscopy to integrate T-system and SR structure with immunofluorescence-based protein levels associated with each membrane, enabling the possibility to directly impose local relative channel expression on these structures.

This model also indicated the importance of realistic (and variable) dyad distribution on the specific dynamics of both CaT alternans and SCRE, providing spatial constraints on the randomness of both. These results are further supported by the recent sub-micron, type 1 model of Hoang-Trong et al. (2021) which included experimentally influenced spatial distributions of multiple Ca²⁺-flux channels, explicit modeling of transverse- and axial-tubules, and Ca²⁺-CaM interactions. This study highlighted the importance of CRU distribution and the presence of rogue (non-junctional or orphaned) RyRs on Ca²⁺-spark propagation and wave dynamics. The model was able to reproduce regenerative Ca²⁺ waves at high Ca²⁺-overload conditions, emerging from the same location on subsequent instances, and the implementation employed several techniques to reduce computational load and memory requirements, enabling efficient implementation of GPU solvers.

Rajagopal et al. (2015, 2018), Ghosh et al. (2018), and Hussain et al. (2018) present a rather different approach and focus to the other whole-cell computational models so-far described, with a higher level of structural detail accounted for Figure 7B. The 3D computational model is generated from images of the myofibrils, mitochondria and RyR clusters; Data from different sources and of different resolutions (e.g., 3D electron microscopy and high-contrast confocal) were fused through spatial statistics techniques (Illian et al., 2008; Theakston et al., 2010). These models focus on many more details of the structure of myocytes and their impact on local Ca²⁺ buffering and regulation, including that of the contractile apparatus and mitochondria, and simulations focus more on the upstroke of the CaT than on long-term homeostasis. The methods for processing high resolution imaging data to generate meshes for simulation are the most advanced in the field, and the potential of these models to understand the super-resolution features of channel distribution, co-localization and Ca²⁺ regulation is currently unparalleled. It would be a significant achievement to develop computational models which account for this level of structural detail that are also sufficiently efficient to simulate long-term dynamics, homeostasis, and Ca²⁺-voltage coupling in a whole-cell.

Simplified, Minimal, and Tissue Models

The complexity of the spatio-temporal models described above hinders ease of analysis and extraction of fundamental mechanisms, and also precludes the efficient cellular simulations required to model hundreds, thousands or millions of cells in cardiac tissue, due to the high computational load of these detailed models. Thus, approaches are required to simplify these complex descriptions into easily analyzable systems and/or efficient computational models while preserving the underlying stochastic dynamics and the emergent phenomena therein. Almost 20-years ago, Shiferaw et al. (2003) presented a number of simplifications to describe spatial Ca²⁺ handling which underlie many of the developments since. Various different approaches have been used including those which explicitly model the cell as a spatial structure but with simplified components of Ca²⁺ release and propagation, and those which develop entirely non-spatial descriptions (Figure 8A). This section focuses primarily on theoretical and numerical approaches which enable large-scale tissue simulations. Analytical and statistical descriptions have also been presented but are not described in detail here; the reader is referred to Rovetti et al. (2007) and Asfaw et al. (2013), for example, where a mean first-passage-time approach was used to demonstrate how Ca²⁺-release depends on local properties within microdomains and to quantify measures of how events synchronize in tissue.

FIGURE 8

Figure 8. Multi-scale modeling approaches. (A) Illustration of different approaches for reducing the computational intensity of models of sub-cellular Ca²⁺ handling. The center, upper panel illustrates a fully detailed 3D cell model, with biophysically detailed descriptions of gating, fluxes and diffusion. Below this panel are illustrated mathematical approaches for model reduction, either resulting in simplified spatial models (left) or non-spatial models (right). The left-most panel indicates how analysis of AP and CaT dynamics on a beat-to-beat basis can lead to the development of a non-spatial, iterated map model. The right-most panel illustrates how analysis of SCRE statistics enables the parameterization of analytical Ca²⁺ release functions which can be imposed in a traditional, non-spatial cell model. (B) Illustration of bi-directional coupling between focal excitations and re-entry. Left: spontaneous excitations interact with each-other to cause conduction block and re-entry, modified from Liu et al. (2015). Right: illustration of the emergence of focal excitations following self-termination after a period of re-entry, from Colman (2019). (C) Structurally detailed tissue models imeplementing SCRE. Left: probability of SCRE as a function of SR-Ca²⁺ in single-cell and 1D, 2D, and 3D with or without the pseudo-1D Purkinje System (PS), with the locations of focal exictations shown in the lower panel, modified from Campos et al. (2015); right: demonstration of a focal excitation emerging within the isthmus of an infarct borderzone, modified from Campos et al. (2018).

Spatial Simplified Models

Williams et al. (2007, 2008) pioneered methods for efficient simulation of stochastic Ca²⁺ dynamics while preserving the importance of local control and features such as graded release. The first study (Williams et al., 2007) developed a probability density approach for modeling local control of CICR and compared results to a Monte-Carlo simulation regarding both validation and computational efficiency. As a simplification based on the “all-or-nothing” response of RyRs within a CRU, the RyRs in each CRU were described as a single “megachannel” i.e., treating them as a single channel which can only occupy a single state at once. Describing this megachannel with a simple two-state model enabled combination with a two-state model of the LTCCs to develop a four-state, minimal model of the state of the CRU. The probability density approach was then applied in the limit that the number of CRUs in the myocyte is sufficiently large (>5,000, which is less than the ∼20,000 observed in cardiomyocytes) by describing the probability that one would find a randomly sampled CRU in a specific state (S), associated with local dyadic Ca²⁺ (Ca_ds) and Ca_SR concentrations:

${\mathrm{\rho }}^i(C{a}_{ds},C{a}_{SR},t)dC{a}_{ds}dC{a}_{SR}=\begin{array}{c}\hfill \mathrm{Pr}\{C{a}_{ds}<\tilde{C}{a}_{ds}(t)<C{a}_{ds}+dC{a}_{ds}\text{and}\hfill \\ \hfill C{a}_{SR}<\tilde{C}{a}_{SR}(t)<C{a}_{SR}+dC{a}_{SR}\text{and}\hfill \\ \hfill \tilde{S}(t)=i\}\hfill \end{array}(42)$

Where i is an index that runs over the four CRU states and tildes indicate random quantities. This must satisfy advection-reaction equations in order to make it equivalent to the Monte-Carlo approach. This approximation was combined with a deterministic solution to the bulk concentrations and a high resolution FDM approximation of spatial dynamics. The authors also derived a univariate approximation to the model based on the marginal density of Ca_SR jointly distributed with the CRU state, although the reader is referred to the original publication for further details beyond the scope of this overview. Importantly, they demonstrated that the approximation agreed strongly with the Monte-Carlo approach, which converges to this result so long as it contains a realistically large number of CRUs. The resulting method is ∼500 times more computationally efficient than the Monte-Carlo approach.

Subsequently, the approach was generalized to describe situations where the dynamics of Ca_ds are much faster than Ca_SR (Williams et al., 2008). The minimal four-state model of the CRU (two states for the RyR megachannel and two for the LTCCs) was expanded to a 12-state model, by extending the two-state RyR model to a six-state model that includes highly cooperative opening and an SR-Ca²⁺ dependence such that depletion of the SR-Ca²⁺ reduces the open probability. The model builds on the univariate approximation presented in the previous study, employing a moment-closure approach truncated at the second-order. The resulting method similarly agrees with the Monte-Carlo approach but with a substantial increase in computational efficiency by a factor of ∼10,000.

Chen et al. (2011) developed a simplified spatial model by imposing a number of reductions to a full description of spatial Ca²⁺ handling. The following two-assumptions drove model development: (1) Whereas Ca²⁺-waves nucleate at some location within the cell and propagate rapidly in the transverse direction, the greater length of the cell compared to its width implies that most Ca²⁺-wave propagation is approximately planar in the longitudinal direction; and (2) in general during a Ca²⁺-wave, all CRUs within a sarcomere are activated at approximately the same time. These assumptions enable a 1D lattice model to be constructed where each node represents a sarcomere (Figure 8A). The state of each sarcomere can be simplified into either being in a non-spark or spark condition [represented by 0 and 1 respectively and conceptually similar to the use of a single RyR megachannel as described in Williams et al. (2007, 2008) above]. The rate of spark recruitment was defined by:

$R(t)=R_{ICaL}(t)+R_{SCRE}(t)(43)$

Where R_ICaL(t) is the recruitment rate due to activation by I_CaL, and R_SCRE(t) is the rate of recruitment via spontaneous sparks and waves. Assuming the release flux (J_rel) associated with each spark is approximately exponential, the spark-rate, R(t), can be associated with J_rel via the following differential equation:

$\frac{d{J}_{rel}(t)}{dt}=g\cdot C{a}_{SR}(t)\cdot R(t)-\frac{J_{rel}(t)\cdot \left(1-{\mathrm{\tau }}_d\left[dC{a}_{SR}\mathrm{/}dt\right]/C{a}_{SR}\right)}{{\mathrm{\tau }}_d}(44)$

Where τ_d is the time-constant of the exponential function assumed to describe the spark flux. The rates of spontaneous spark initiation and a transmission-time and -probability [i.e., the components of R_SCRE(t)] are primarily regulated by Ca_SR and were fit to experimental data to describe Ca²⁺-wave velocity, enabling many features of whole-cell SCRE to be captured at a significantly reduced cost.

In Hernandez-Hernandez et al. (2015) a similar model was constructed in 2D (Figure 8A). The state of each CRU (0 or 1 for non-spark or spark, respectively) was described by the following simple reaction:

$0\underset{\mathrm{\beta }}{\overset{\mathrm{\alpha }=f(C{a}_i)}{\rightleftharpoons }}1(45)$

Inter-CRU coupling was described using a spatially exponential function such that the influence of one CRU on its neighbors decays rapidly as distance increases:

$h_{ij}=r_j{e}^{-{\left|x_i-x_j\right|}^2/l^2}(46)$

Where r_j is the Ca²⁺ released at site j, x_i, and x_j are the locations of the two coupled CRUs and l is the diffusive length-scale. This model was used to evaluate the influence of CRU connectivity on Ca²⁺-wave nucleation and propagation.

In Romero et al. (2019) a 2D model was developed where the dynamics of each CRU were described by a non-linear map which relates Ca²⁺ concentrations from one beat to the next (rather than solving concentrations on a small time-step within each beat). This is a spatial analog of the iterative-map approach described in the next sub-section. Concentrations for beat n+1 are functions of concentrations and fluxes during beat n:

$C{a}_{SR}^x\left(n+1\right)=C{a}_{SR}^x(n)-R^x(n)+U^x(n)(47)$

$C{a}_i^x\left(n+1\right)=C{a}_i^x(n)+R^x(n)-U^x(n)(48)$

Where Ca_SR^x and Ca_i^x are the SR and intracellular Ca²⁺ concentrations at the spatial point x (=i,j in 2D), R^x(n) is the total Ca²⁺ released from the SR at point x during beat n and U^x(n) is the total Ca²⁺ pumped back into the SR at point x during beat n. As a simplification, it is assumed that the total Ca²⁺ is conserved (Ca_SR + Ca_i = C) and can be thusly normalized to 1 (arbitrary units). The stochastic dependence of intracellular Ca²⁺ release (R^x) on voltage can be incorporated by accounting for the probability of release, and intracellular uptake is given as a function of the CaT peak. Ca²⁺ diffusion was then described by accounting for the average Ca²⁺ over the nearest neighbors in the 2D lattice of CRUs.

In Cantalapiedra et al. (2017) a spatial model was developed which included simplified descriptions of the LTCCs and RyRs through exploiting symmetries and other factors to substantially reduce the number of equations/parameters governing these dynamics; spatial coupling was described explicitly as in a CRU-grid fully spatial model. The model was used to study the influence of SR-Ca²⁺ load and RyR refractoriness on the dynamics of CaT alternans.

Non-spatial Simplified Models

Further simplification can be achieved by developing non-spatial (or “0D”) approximations, which have no explicit description of spatial diffusion in the sub-cellular volume and present the opportunity for the largest increases in computational efficiency. Several different techniques have been proposed.

Reduction to (Semi-)deterministic Models

In Chen et al. (2012) a model was developed which describes single-cell SCRE dynamics using a single simple two-state reaction scheme that is able to capture the statistics of SCRE timing and magnitude. This is the same basic scheme as presented by Hernandez-Hernandez et al. (2015) and equation (45), but is now solved in the deterministic limit describing the whole-cell (Figure 8A). This was applied in a 1D model of tissue to analyze the relationship between the statistics of single-cell SCRE (pertaining to timing, magnitude and duration) and the emergence of spontaneous focal (or ectopic) excitation in tissue.

In studies by Shiferaw et al. (2018, 2020) a phenomenological, population-dynamics-like model of spark recruitment was developed that matched behavior of the 3D cell model. Here, rather than modeling CRU dynamics explicitly, the number of active CRUs (or number of sparks) were tracked and dynamically evolves dependent on the number of sparks initiated (Δn⁺) and extinguished (Δn^–) at each time-step (Δt):

$n_i(t+\mathrm{△}t)=n_i(t)+\mathrm{△}{n}_i^{+}-\mathrm{△}{n}_i^{-}(49)$

The number of sparks initiated was determined by the rate at which sparks are recruited (either spontaneous or triggered) and stochasticity was maintained through the use of random number sampling. This used a similar approach to individual CRU recruitment to that of Chen et al. (2012) and Hernandez-Hernandez et al. (2015) but now included distinction between junctional and non-junctional sparks [both described by equation (45) but with different values for the transition rates], enabling the impact of heterogeneous and variable atrial T-systems to be captured in this reduced model. The approximation could reproduce both CaT alternans and SCRE. The model was used to study synchronization of SCRE in atrial tissue, developing focal excitations, conduction-block, and non-stable re-entrant-like excitation patterns.

Iterated-Map Models

Qu et al. (2007) developed an iterated-map model of CRU activity, relating Ca²⁺ and voltage properties on a beat-to-beat basis (Figure 8A), a whole-cell equivalent to the more recent approach presented in Romero et al. (2019). The approach can be briefly summarized by the following equations. Firstly, the APD can depend on both the diastolic interval (DI) and peak of the CaT (Ca^peak). Thus, the APD at the current cycle (n) depends on the previous cycle (n-1) and is given by:

$AP{D}_n=f_1\left(D{I}_{n-1}\right)+f_2\left(C{a}_n^{peak}\right)AP{D}_n(50)$

Where the functions of DI and intracellular Ca²⁺ have been separated: f₁ is the APD restitution function and f₂ accounts for the coupling strength between Ca²⁺ and APD (which can be positive or negative). The peak Ca²⁺ concentration at cycle n can be given by the sum of the diastolic Ca²⁺ from the previous cycle (Ca_n–1) and the total Ca²⁺ released from the SR in the current cycle (r_n):

$C{a}_n^{peak}=C{a}_{n-1}+r_n(51)$

The total Ca²⁺ released from the SR, r_n, is given by:

$r_n=f_3\left(D{I}_{n-1}\right)f_4\left(CaS{R}_{n-1}^{load}\right)(52)$

Where f₃ describes the restitution properties of SR-Ca²⁺ release (e.g., RyR refractoriness) and f₄ describes the dependence of Ca²⁺ release on the SR-Ca²⁺ load. The reader is referred to the original publication for full details on the parameters and functions involved. This model reproduced the non-linear dynamics of Ca²⁺ handling including CaT alternans, and this simplification helped to develop a unified theory of CaT alternans in cardiac cells (Qu et al., 2016). Furthermore, this vast simplification in both space and time produces exceptionally efficient simulations, computationally less intense than standard, non-spatial common-pool models of cardiac cellular electrophysiology.

Models of Imposed, Stochastic Spontaneous Ca²⁺ Release Functions

As an alternative to the above approaches, it is also possible to control spontaneous CaTs in otherwise deterministic cell models by imposing (or clamping) SCRE waveforms (Figure 8A). In Xie et al. (2010), the intracellular release flux (J_rel) associated with SCRE was controlled by imposing a waveform defined by two sigmoidal functions, the parameters of which determined the timing, duration and magnitude of SCRE:

$\begin{array}{c}\hfill J_{rel}^{spont}=G_{spont}{\left(1+e^{-(t-t_0)/{\mathrm{\tau }}_1}\right)}^{-1}\hfill \\ \hfill {\left(1+e^{(t-t_0)/{\mathrm{\tau }}_2}\right)}^{-1}\left(\frac{v_{ds}}{v_{jsr}}C{a}_{SR}-C{a}_{ds}\right)\hfill \end{array}(53)$

Where G_spont is a rate constant (set to 0.0674 ms^–1 in the original study), t₀ was set to 425 ms and τ₁ = 10 ms and τ₂ = 30 ms. The model was used to determine the minimum number of cells undergoing DADs in various tissue conditions in order for this to manifest as a focal excitation. Whereas not performed in the original study, the timing and duration parameters (t_0, τ₁, τ₂) could be randomly sampled from distributions to reproduce stochasticity in independent cellular SCRE. This approach was used in Liu et al. (2015), sampling from Gaussian distributions, in order to study the dynamics of independently timed DADs in tissue, revealing mechanisms of synchronization into focal excitation and DAD-mediated conduction block (e.g., Figure 8B). Liu et al. (2016) and Ko et al. (2017) used a similar but further simplified approach in which the Ca_i associated with SCRE was directly controlled, described by a Gaussian-shaped function with equivalent parameters to control timing, magnitude and duration:

$C{a}_i^{spont}=A{e}^{-{(t-t_0)}^2/2{\mathrm{\sigma }}^2}(54)$

Where A sets the maximum amplitude, t₀ sets the latency (timing) and σ sets the duration. In Colman et al. (2017a) and Colman (2019) an approach was developed which used sigmoidal functions to control the RyR-state directly, reproducing different shapes of long- and short-release waveforms. For short (spike-like) waveforms:

$N_{RyR\mathrm{\_}O}=N_{RyR\mathrm{\_}O}^{peak}{\left(1+e^{-\left(t-t_1\right)/k_1}\right)}^{-1}{\left(1+e^{\left(t-t_2\right)/k_2}\right)}^{-1}(55)$

$t_1=t_i+0.5\left(t_p-t_i\right)(56)$

$t_2=t_p+0.5\left(t_f-t_p\right)(57)$

$k_1=0.1689(t_p-t_i)+0.00255(58)$

$k_2=0.1689\left(t_f-t_p\right)+0.00255(59)$

where t_i is the initiation time (equivalent to latency, t₀, in the above models) of the SCRE, t_f is the end time (duration, λ, thus = t_f-t_i), t_p is the time of the peak of the waveform and N_{RyR_O}^peak is the peak proportion of open RyRs. In this model, N_{RyR_O} [equation (55)] replaces the “O” in the J_rel equation [equation (5)] and so J_rel and the CaT are allowed to dynamically evolve according to their deterministic functions; the magnitude of the CaT associated with the SCRE is therefore dependent on N_{RyR_O}^peak, the SR-Ca²⁺ and J_rel maximal flux rate. This model demonstrated feedback between re-entry and focal excitation in which the rapid activation during re-entry loads the SR-Ca²⁺ to promote focal excitations following re-entry termination (Figure 8B).

In all implementations, the parameters defining the distributions which describe SCRE statistics could be set as a function of environmental variables, such as SR-Ca²⁺, enabling the simplified cell model to respond to pacing with variable SCRE statistics in congruence with the dynamics of spatial Ca²⁺ cell models. All of these models employ an algorithm to determine if Ca_i, J_rel or RyR-state is controlled by the deterministic cell model or undergoes the imposed SCRE clamped waveform, enabling integration with dynamically evolving deterministic cell models in both single-cell and tissue simulations. The computational efficiency of these models is comparable to that of standard common-pool models of cardiac electrophysiology; the largest computational cost is the generation of random numbers (where one must indeed be careful with implementations of parallelization in tissue simulations), but the cost of this inclusion is smaller than the typical differences in computational efficiency between common-pool models which feature a different number of components and governing equations.

Estimating Probabilities of Rare Events

Walker et al. (2017) implemented a study in which 3D cell models were coupled in a 1D tissue strand (or fiber). They investigated the mediators of Ca²⁺ waves and DADs in single cells, and used the 1D fiber model to translate these features to tissue activity. From this, they developed a spatial-average filtering model which aimed to estimate V_m from intracellular release fluxes, enabling the estimation of the probabilities of “extreme” (i.e., rare) events in which multiple cells synchronously undergo large-scale SCRE, i.e., the requirements for focal excitations. This type of approach is powerful and important because the generation of a serious arrhythmia in an individual is often a very rare event that cannot be robustly or consistently captured in simulations of a generally normally functioning heart. In agreement with other studies (Liu et al., 2016; Campos et al., 2017; Colman, 2019), they found reduced I_K1 and inter-cellular coupling to be important for enabling SCRE to overcome electrotonic load and promote focal excitations.

Whole-Heart Models; Integration With Tissue Imaging

A few studies have integrated these reduced models of stochastic sub-cellular Ca²⁺ handling into models of the whole atria or ventricle in order to study the interaction between cellular function and tissue structure in controlling the emergence and dynamics of arrhythmogenic triggers. This is perhaps taken the furthest in a series of papers by Campos et al. (2015, 2017, 2018, 2019). The initial study (Campos et al., 2015) combined the phenomenological model of SCRE as proposed in Chen et al. (2011) with a full bi-ventricular 3D model which included a description of the Purkinje network. This model demonstrated that focal excitations were preferentially located to the Purkinje network due to the reduced electrotonic load in these pseudo-1D-strands, with focal excitations increasingly repressed as dimensionality increased from 1D to 3D (Figure 8C). Campos et al. (2017) subsequently incorporated a description of sodium-channel dysfunction, which promoted focal excitations and conduction block leading to re-entry. The two more recent studies (Campos et al., 2018, 2019) now combined these analyses with structural remodeling associated with infarcts in both idealized 2D sheets (Figure 8C) and realistic geometries in the 3D bi-ventricular model, demonstrating that the macroscopic and microscopic anatomy of the infarct region could promote both focal excitation and re-entry and highlighting the mechanisms by which fibrosis could increase the probability of focal excitations in these conditions.

These studies demonstrate structural features which can co-localize both focal excitation and re-entrant excitation. In Colman (2019) a purely functional mechanism which can co-locate focal and re-entrant excitation was revealed: the in-excited core of re-entrant excitation lead to substantially longer latency times in this region, enabling focal excitation to preferentially emerge from this same location (Figure 8B); interaction of focal excitation with the tail of the previous re-entrant excitation could lead to highly asymmetric focal excitations which themselves may degenerate back into re-entry with a core in the same approximate location. Thus, these studies are revealing both structural and functional mechanisms which spatially relate focal and re-entrant excitation.

Drawing on Physical Analogies for Mechanistic Explanation

Other works have used detailed or simplified models to make analogies that relate cardiac Ca²⁺-handling phenomena to other physical phenomena, which may offer further insight into fundamental underlying mechanisms or provide more predictive power.

Alvarez-Lacalle et al. (2015) used both a 2D sub-cellular spatial model and simplified descriptions in order to analyze the dynamics of CaT alternans. Through scaling analysis of correlations near the transition to alternans it was demonstrated that CaT alternans could be described as an order-disorder phase transition, leading to an analogy to the Ising model of ferromagnetism in statistical mechanics. This was further generalized to describe more features of Ca²⁺ dynamics, including coupling with voltage, in Romero et al. (2019), where the analogy was extended to the more general Potts model. Both discontinuous first-order phase transitions and second-order continuous phase transitions (alternans) were observed to emerge under different conditions in cardiac cellular dynamics, further supporting the idea that statistical mechanics tools may be valuable for understanding cardiac function.

Conesa et al. (2020) presented a novel approach to understanding steady-state activity through analysis of a single beat which is not in homeostatic balance, by reduction to two-variable general equilibrium conditions in analogy to models of macro-economics. Such an approach can help to explain the complex and often counter-intuitive features of Ca²⁺-handling and offers substantial predictive power without the requirement for computer-intensive simulations.

Summary and Conclusion

Recent advances in the robustness, complexity and sophistication of computational models of spatial Ca²⁺-handling from the nanometer-scale to the whole-heart scale are enabling advanced simulations to be performed to reveal fundamental properties of cardiac ECC in both health and disease. Multiple different approaches have been explored to describe nanodomain dynamics, model inter-CRU coupling, implement experimental imaging data, and translate models to the whole-heart scale. There are fundamental differences in the approaches and structure of these models, and in some cases model behavior can differ substantially (e.g., in triggered Ca²⁺ wave propagation). Nevertheless, the general agreement between models is encouraging and the availability of multiple different models provides the opportunity to comprehensively test hypotheses and explore fundamental theories. As computational power increases, experimental imaging data improves, and more powerful coarse-graining techniques are developed, the relevance, scope and power of these models will only continue to increase.

There are a number of challenges and avenues for future development and innovation. Robust validation of the models remains a major difficulty, not-least because of the challenges in obtaining sufficient and congruent experimental data to validate the many interacting model components and emergent functional phenomena. Validation of the governing RyR models themselves, in isolation and in the context of the spatial-nanodomain, is non-trivial - important features of Ca²⁺ sparks, such as the spatial FWHM, do not necessarily match experimental observations. This is not the only difficulty associated with obtaining a realistic and well-validated description of RyR gating: Integration of spark properties in the context of whole-cell homeostasis is also highly challenging due to the interaction with the wider model system and further constraints on model stability and long-term dynamics.

At the whole-cell scale, many approaches have been implemented to describe spatial Ca²⁺ coupling, corresponding to fundamentally different underlying model structures. Each approach has associated features of the CaT and Ca²⁺-handling dynamics, such as the robustness of Ca²⁺-wave propagation; certainly, not all of these structures and parameter combinations can simultaneously accurately describe real myocytes, and there are therefore fundamental questions about the mechanisms of inter-CRU Ca²⁺ diffusion which remain to be resolved. Nevertheless, the success of these models to explain and provide interpretations into experimental results is highly encouraging, especially in the more recent studies which combine experiment and simulation and are indeed beginning to resolve some of these issues.

One of the most exciting and challenging prospects is true multi-scale model integration. The recent spatial models of single nanodomains highlight the importance of specific RyR arrangement and microstructure of the dyad; incorporating these features into a whole-cell model containing tens of thousands of heterogeneous dyads is far from trivial. The detailed models of the local regions of a single T-tubule reveal important features of local Ca²⁺ gradients and channel distribution – combining these models with detailed representations of RyR arrangement in nanodomains, models of counter-ion fluxes, and the electrophysiology of the T-tubule membrane itself presents exciting prospects for powerful and highly accurate models of local control; again, translating these to the whole-cell scale is associated with a number of challenges. Initial success has been achieved in models attempting to preserve the impact of stochastic spatial Ca²⁺ dynamics in reduced, computationally efficient cell models suitable for tissue simulations, revealing the mechanisms of ectopic excitation, its interaction with re-entry, and dependence on tissue structure. Generalizing these models to naturally capture the dynamics of heterogenous populations of cells, themselves depending on heterogeneous sub-cellular structure, remains a major goal.

Finally, there have been substantial advances in approaches for image-based modeling, both for high-throughput, population-cohort simulations, and for direct integration of experimental structures. Closely related to advances in experimental imaging modalities, this is occurring at multiple spatial scales including the nanodomain, localized sub-cellular regions, whole-cells, and whole-heart. The further development and automation of these experimental-simulation frameworks presents exciting prospects for the true mechanistic analysis of structure-function relationships underlying cardiac electrophysiology from the nanometer to the whole-heart.

Author Contributions

MC conceived the review, drafted and edited the manuscript, and prepared the illustrations. All the authors drafted and edited the manuscript.

Funding

This work was supported by a Medical Research Council, United Kingdom, Strategic Skills Fellowship (Grant Number MR/M014967/1) and Career Development Award (Grant Number MR/V010050/1) awarded to MC, by the Netherlands Organization for Scientific Research NWO/ZonMW Vidi 09150171910029 awarded to JH, by MICINN/AEI through research grants SAF-2017-88019-C3-2-R and PID2020-116927RB-C22 awarded to BE and EA-L, and by the National Institutes of Health (Grant Numbers R01-HL149349 and P01-HL141084) awarded to DS.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to thank the following for helpful conversations and perspectives on this topic: Yohannes Shiferaw (California State University, Northridge), Zhen Song and Zhilin Qu (University of California, Los Angeles), Fernando Campos and Martin Bishop (King’s College London), and Izzy Jayasinghe (The University of Sheffield).

References

Alvarez-Lacalle, E., Cantalapiedra, I. R., Peñaranda, A., Cinca, J., Hove-Madsen, L., and Echebarria, B. (2013). Dependency of calcium alternans on ryanodine receptor refractoriness. PLoS One 8:e55042. doi: 10.1371/journal.pone.0055042

PubMed Abstract | CrossRef Full Text | Google Scholar

Alvarez-Lacalle, E., Echebarria, B., Spalding, J., and Shiferaw, Y. (2015). Calcium alternans is due to an order-disorder phase transition in cardiac cells. Phys. Rev. Lett. 114:108101. doi: 10.1103/PhysRevLett.114.108101

PubMed Abstract | CrossRef Full Text | Google Scholar

Asfaw, M., Alvarez-Lacalle, E., and Shiferaw, Y. (2013). The timing statistics of spontaneous calcium release in cardiac myocytes. PLoS One 8:e62967. doi: 10.1371/journal.pone.0062967

PubMed Abstract | CrossRef Full Text | Google Scholar

Baddeley, D., Jayasinghe, I. D., Lam, L., Rossberger, S., Cannell, M. B., and Soeller, C. (2009). Optical single-channel resolution imaging of the ryanodine receptor distribution in rat cardiac myocytes. Proc. Natl. Acad. Sci. U. S. A. 106, 22275–22280. doi: 10.1073/pnas.0908971106

PubMed Abstract | CrossRef Full Text | Google Scholar

Benson, A. P., Bernus, O., Dierckx, H., Gilbert, S. H., Greenwood, J. P., Holden, A. V., et al. (2011). Construction and validation of anisotropic and orthotropic ventricular geometries for quantitative predictive cardiac electrophysiology. Interface Focus 1, 101–116. doi: 10.1098/rsfs.2010.0005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bers, D. M. (2002). Cardiac excitation–contraction coupling. Nature 415, 198–205. doi: 10.1038/415198a

PubMed Abstract | CrossRef Full Text | Google Scholar

Berti, C., Zsolnay, V., Shannon, T. R., Fill, M., and Gillespie, D. (2017). Sarcoplasmic reticulum Ca2+, Mg2+, K+, and Cl- concentrations adjust quickly as heart rate changes. J. Mol. Cell. Cardiol. 103, 31–39. doi: 10.1016/j.yjmcc.2016.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, N., D’Souza, A., Wang, Y., Piggins, H., Dobrzynski, H., Morris, G., et al. (2019). Circadian rhythm of cardiac electrophysiology, arrhythmogenesis, and the underlying mechanisms. Heart Rhythm 16, 298–307. doi: 10.1016/j.hrthm.2018.08.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Brandenburg, S., Kohl, T., Williams, G. S. B., Gusev, K., Wagner, E., Rog-Zielinska, E. A., et al. (2016). Axial tubule junctions control rapid calcium signaling in atria. J. Clin. Invest. 126, 3999–4015. doi: 10.1172/JCI88241

PubMed Abstract | CrossRef Full Text | Google Scholar

Brette, F., Despa, S., Bers, D. M., and Orchard, C. H. (2005). Spatiotemporal characteristics of SR Ca2+ uptake and release in detubulated rat ventricular myocytes. J. Mol. Cell. Cardiol. 39, 804–812. doi: 10.1016/j.yjmcc.2005.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Campos, F. O., Shiferaw, Y., Prassl, A. J., Boyle, P. M., Vigmond, E. J., and Plank, G. (2015). Stochastic spontaneous calcium release events trigger premature ventricular complexes by overcoming electrotonic load. Cardiovasc. Res. 107, 175–183. doi: 10.1093/cvr/cvv149

PubMed Abstract | CrossRef Full Text | Google Scholar

Campos, F. O., Shiferaw, Y., Vigmond, E. J., and Plank, G. (2017). Stochastic spontaneous calcium release events and sodium channelopathies promote ventricular arrhythmias. Chaos 27:093910. doi: 10.1063/1.4999612

CrossRef Full Text | Google Scholar

Campos, F. O., Shiferaw, Y., Weber dos Santos, R., Plank, G., and Bishop, M. J. (2018). Microscopic Isthmuses and Fibrosis Within the Border Zone of Infarcted Hearts Promote Calcium-Mediated Ectopy and Conduction Block. Front. Phys. 6:57. doi: 10.3389/fphy.2018.00057

CrossRef Full Text | Google Scholar

Campos, F. O., Whitaker, J., Neji, R., Roujol, S., O’Neill, M., Plank, G., et al. (2019). Factors promoting conduction slowing as substrates for block and reentry in infarcted hearts. Biophys. J. 117, 2361–2374. doi: 10.1016/j.bpj.2019.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Cannell, M. B., and Kong, C. H. T. (2012). Local control in cardiac E-C coupling. J. Mol. Cell. Cardiol. 52, 298–303. doi: 10.1016/j.yjmcc.2011.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Cannell, M. B., and Kong, C. H. T. (2017). Quenching the spark: termination of CICR in the submicroscopic space of the dyad. J. Gen. Physiol. 149, 837–845. doi: 10.1085/jgp.201711807

PubMed Abstract | CrossRef Full Text | Google Scholar

Cannell, M. B., Kong, C. H. T., Imtiaz, M. S., and Laver, D. R. (2013). Control of Sarcoplasmic Reticulum Ca2+ Release by Stochastic RyR Gating within a 3D Model of the Cardiac Dyad and Importance of Induction Decay for CICR Termination. Biophys. J. 104, 2149–2159. doi: 10.1016/j.bpj.2013.03.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Cantalapiedra, I. R., Alvarez-Lacalle, E., Peñaranda, A., and Echebarria, B. (2017). Minimal model for calcium alternans due to SR release refractoriness. Chaos 27:093928. doi: 10.1063/1.5000709

CrossRef Full Text | Google Scholar

Chen, W., Aistrup, G., Wasserstrom, J. A., and Shiferaw, Y. (2011). A mathematical model of spontaneous calcium release in cardiac myocytes. Am. J. Physiol.-Heart Circ. Physiol. 300, H1794–H1805. doi: 10.1152/ajpheart.01121.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., Asfaw, M., and Shiferaw, Y. (2012). The Statistics of Calcium-Mediated Focal Excitations on a One-Dimensional Cable. Biophys. J. 102, 461–471. doi: 10.1016/j.bpj.2011.12.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., Wasserstrom, J. A., and Shiferaw, Y. (2009). Role of coupled gating between cardiac ryanodine receptors in the genesis of triggered arrhythmias. Am. J. Physiol.-Heart Circ. Physiol. 297, H171–H180. doi: 10.1152/ajpheart.00098.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, H., Lederer, W. J., and Cannell, M. B. (1993). Calcium sparks: elementary events underlying excitation-contraction coupling in heart muscle. Science 262, 740–744. doi: 10.1126/science.8235594

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, Y., Yu, Z., Hoshijima, M., Holst, M. J., McCulloch, A. D., McCammon, J. A., et al. (2010). Numerical Analysis of Ca2+ Signaling in Rat Ventricular Myocytes with Realistic Transverse-Axial Tubular Geometry and Inhibited Sarcoplasmic Reticulum. PLoS Comput. Biol. 6:e1000972. doi: 10.1371/journal.pcbi.1000972

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, J. D., Caldwell, J. L., Horn, M. A., Bode, E. F., Richards, M. A., Hall, M. C. S., et al. (2015). Perturbed atrial calcium handling in an ovine model of heart failure: potential roles for reductions in the L-type calcium current. J. Mol. Cell. Cardiol. 79, 169–179. doi: 10.1016/j.yjmcc.2014.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Colman, M. A. (2019). Arrhythmia mechanisms and spontaneous calcium release: Bi-directional coupling between re-entrant and focal excitation. PLoS Comput. Biol. 15:e1007260. doi: 10.1371/journal.pcbi.1007260

PubMed Abstract | CrossRef Full Text | Google Scholar

Colman, M. A., Aslanidi, O. V., Kharche, S., Boyett, M. R., Garratt, C., Hancox, J. C., et al. (2013). Pro-arrhythmogenic effects of atrial fibrillation-induced electrical remodelling: insights from the three-dimensional virtual human atria. J. Physiol. 591, 4249–4272. doi: 10.1113/jphysiol.2013.254987

PubMed Abstract | CrossRef Full Text | Google Scholar

Colman, M. A., Holmes, M., Whittaker, D. G., Jayasinghe, I., and Benson, A. P. (2020). Multi-scale approaches for the simulation of cardiac electrophysiology: I - Sub-cellular and stochastic calcium dynamics from cell to organ. Methods 185, 49–59. doi: 10.1016/j.ymeth.2020.02.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Colman, M. A., Perez Alday, E. A., Holden, A. V., and Benson, A. P. (2017a). Trigger vs. Substrate: multi-Dimensional Modulation of QT-Prolongation Associated Arrhythmic Dynamics by a hERG Channel Activator. Front. Physiol. 8:757. doi: 10.3389/fphys.2017.00757

PubMed Abstract | CrossRef Full Text | Google Scholar

Colman, M. A., Pinali, C., Trafford, A. W., Zhang, H., and Kitmitto, A. (2017b). A computational model of spatio-temporal cardiac intracellular calcium handling with realistic structure and spatial flux distribution from sarcoplasmic reticulum and t-tubule reconstructions. PLoS Comput. Biol. 13:e1005714. doi: 10.1371/journal.pcbi.1005714

PubMed Abstract | CrossRef Full Text | Google Scholar

Colman, M. A., Sarathy, P. P., MacQuiaide, N., and Workman, A. J. (2016). “A new model of the human atrial myocyte with variable T-tubule organization for the study of atrial fibrillation,” in 2016 Computing in Cardiology Conference (CinC), (Vancouver), 221–224. doi: 10.23919/CIC.2016.7868719

CrossRef Full Text | Google Scholar

Conesa, D., Echebarria, B., Peñaranda, A., Cantalapiedra, I. R., Shiferaw, Y., and Alvarez-Lacalle, E. (2020). Two-variable nullcline analysis of ionic general equilibrium predicts calcium homeostasis in ventricular myocytes. PLoS Comput. Biol. 16:e1007572. doi: 10.1371/journal.pcbi.1007572

PubMed Abstract | CrossRef Full Text | Google Scholar

Crocini, C., Coppini, R., Ferrantini, C., Yan, P., Loew, L. M., Tesi, C., et al. (2014). Defects in T-tubular electrical activity underlie local alterations of calcium release in heart failure. Proc. Natl. Acad. Sci. U. S. A. 111, 15196–15201. doi: 10.1073/pnas.1411557111

PubMed Abstract | CrossRef Full Text | Google Scholar

Crossman, D. J., Ruygrok, P. N., Ruygrok, P. R., Soeller, C., and Cannell, M. B. (2011). Changes in the organization of excitation-contraction coupling structures in failing human heart. PLoS One 6:e17901. doi: 10.1371/journal.pone.0017901

PubMed Abstract | CrossRef Full Text | Google Scholar

Dibb, K. M., Louch, W. E., and Trafford, A. W. (2022). Cardiac Transverse Tubules in Physiology and Heart Failure. Annu. Rev. Physiol. 84. doi: 10.1146/annurev-physiol-061121-040148

PubMed Abstract | CrossRef Full Text | Google Scholar

D’Souza, A., Wang, Y., Anderson, C., Bucchi, A., Baruscotti, M., Olieslagers, S., et al. (2021). A circadian clock in the sinus node mediates day-night rhythms in Hcn4 and heart rate. Heart Rhythm 18, 801–810. doi: 10.1016/j.hrthm.2020.11.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisner, D. A., Caldwell, J. L., Kistamás, K., and Trafford, A. W. (2017). Calcium and Excitation-Contraction Coupling in the Heart. Circ. Res. 121, 181–195. doi: 10.1161/CIRCRESAHA.117.310230

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisner, D. A., Kashimura, T., Venetucci, L. A., and Trafford, A. W. (2009). From the ryanodine receptor to cardiac arrhythmias. Circ. J. 73, 1561–1567. doi: 10.1253/circj.cj-09-0478

PubMed Abstract | CrossRef Full Text | Google Scholar

Gadeberg, H. C., Bond, R. C., Kong, C. H. T., Chanoit, G. P., Ascione, R., Cannell, M. B., et al. (2016). Heterogeneity of T-Tubules in Pig Hearts. PLoS One 11:e0156862. doi: 10.1371/journal.pone.0156862

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaur, N., and Rudy, Y. (2011). Multiscale modeling of calcium cycling in cardiac ventricular myocyte: macroscopic consequences of microscopic dyadic function. Biophys. J. 100, 2904–2912. doi: 10.1016/j.bpj.2011.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh, S., Tran, K., Delbridge, L. M. D., Hickey, A. J. R., Hanssen, E., Crampin, E. J., et al. (2018). Insights on the impact of mitochondrial organisation on bioenergetics in high-resolution computational models of cardiac cell architecture. PLoS Comput. Biol. 14:e1006640. doi: 10.1371/journal.pcbi.1006640

PubMed Abstract | CrossRef Full Text | Google Scholar

Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22, 403–434. doi: 10.1016/0021-9991(76)90041-3

CrossRef Full Text | Google Scholar

Greene, D., and Shiferaw, Y. (2021). Mechanistic link between CaM-RyR2 interactions and the genesis of cardiac arrhythmia. Biophys. J. 120, 1469–1482. doi: 10.1016/j.bpj.2021.02.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Greenstein, J. L., and Winslow, R. L. (2002). An integrative model of the cardiac ventricular myocyte incorporating local control of Ca2+ release. Biophys. J. 83, 2918–2945. doi: 10.1016/S0006-3495(02)75301-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Hake, J., Edwards, A. G., Yu, Z., Kekenes-Huskey, P. M., Michailova, A. P., McCammon, J. A., et al. (2012). Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit. J. Physiol. 590, 4403–4422. doi: 10.1113/jphysiol.2012.227926

PubMed Abstract | CrossRef Full Text | Google Scholar

Hake, J., Kekenes-Huskey, P. M., and McCulloch, A. D. (2014). Computational modeling of subcellular transport and signaling. Curr. Opin. Struct. Biol. 25, 92–97. doi: 10.1016/j.sbi.2014.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatano, A., Okada, J., Hisada, T., and Sugiura, S. (2012). Critical role of cardiac t-tubule system for the maintenance of contractile function revealed by a 3D integrated model of cardiomyocytes. J. Biomech. 45, 815–823. doi: 10.1016/j.jbiomech.2011.11.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatano, A., Okada, J., Washio, T., Hisada, T., and Sugiura, S. (2011). A Three-Dimensional Simulation Model of Cardiomyocyte Integrating Excitation-Contraction Coupling and Metabolism. Biophys. J. 101, 2601–2610. doi: 10.1016/j.bpj.2011.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Heijman, J., Erfanian Abdoust, P., Voigt, N., Nattel, S., and Dobrev, D. (2016). Computational models of atrial cellular electrophysiology and calcium handling, and their role in atrial fibrillation. J. Physiol. 594, 537–553. doi: 10.1113/JP271404

PubMed Abstract | CrossRef Full Text | Google Scholar

Heijman, J., Volders, P. G. A., Westra, R. L., and Rudy, Y. (2011). Local control of β-adrenergic stimulation: effects on ventricular myocyte electrophysiology and Ca2+-transient. J. Mol. Cell. Cardiol. 50, 863–871. doi: 10.1016/j.yjmcc.2011.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernandez-Hernandez, G., Alvarez-Lacalle, E., and Shiferaw, Y. (2015). Role of connectivity and fluctuations in the nucleation of calcium waves in cardiac cells. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 92:052715. doi: 10.1103/PhysRevE.92.052715

PubMed Abstract | CrossRef Full Text | Google Scholar

Herzel, H. (1991). Risken, H., The Fokker-Planck-Equation. Methods of Solution and Applications. Second edition. Berlin etc., Springer-Verlag 1989. XIV, 472 pp., 95 figs., DM 98,—. ISBN3-540-50498-2 (Springer Series in Synergetics 18). ZAMM - J. Appl. Math. Mech. Z. Für Angew. Math. Mech. 71, 67–67. doi: 10.1002/zamm.19910710128

CrossRef Full Text | Google Scholar

Hinch, R. (2004). A mathematical analysis of the generation and termination of calcium sparks. Biophys. J. 86, 1293–1307. doi: 10.1016/S0006-3495(04)74203-4

CrossRef Full Text | Google Scholar

Hoang-Trong, T. M., Ullah, A., Lederer, W. J., and Jafri, M. S. (2021). A Stochastic Spatiotemporal Model of Rat Ventricular Myocyte Calcium Dynamics Demonstrated Necessary Features for Calcium Wave Propagation. Membranes 11:989. doi: 10.3390/membranes11120989

PubMed Abstract | CrossRef Full Text | Google Scholar

Holash, R. J., and MacIntosh, B. R. (2019). A stochastic simulation of skeletal muscle calcium transients in a structurally realistic sarcomere model using MCell. PLoS Comput. Biol. 15:e1006712. doi: 10.1371/journal.pcbi.1006712

PubMed Abstract | CrossRef Full Text | Google Scholar

Holmes, M., Hurley, M., Sheard, T., Benson, A., Jayasinghe, I., and Colman, M. (2021). Increased SERCA2a sub-cellular heterogeneity in right-ventricular heart failure inhibits excitation-contraction coupling and modulates arrhythmogenic dynamics. Philos. Trans. R. Soc. B Biol. Sci. (in press). doi: 10.1098/rstb.2021-0317 [Epub ahead of print].

CrossRef Full Text | Google Scholar

Hüser, J., Lipsius, S. L., and Blatter, L. A. (1996). Calcium gradients during excitation-contraction coupling in cat atrial myocytes. J. Physiol. 494, 641–651. doi: 10.1113/jphysiol.1996.sp021521

PubMed Abstract | CrossRef Full Text | Google Scholar

Hussain, A., Ghosh, S., Kalkhoran, S. B., Hausenloy, D. J., Hanssen, E., and Rajagopal, V. (2018). An automated workflow for segmenting single adult cardiac cells from large-volume serial block-face scanning electron microscopy data. J. Struct. Biol. 202, 275–285. doi: 10.1016/j.jsb.2018.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Iaparov, B. I., Zahradnik, I., Moskvin, A. S., and Zahradníková, A. (2021). In silico simulations reveal that RYR distribution affects the dynamics of calcium release in cardiac myocytes. J. Gen. Physiol. 153:e202012685. doi: 10.1085/jgp.202012685

PubMed Abstract | CrossRef Full Text | Google Scholar

Illian, J., Penttinen, A., Stoyan, H., and Stoyan, D. (2008). Statistical Analysis and Modelling of Spatial Point Patterns. Chichester: John Wiley & Sons Ltd, doi: 10.1002/9780470725160

CrossRef Full Text | Google Scholar

Jayasinghe, I., Clowsley, A. H., Lin, R., Lutz, T., Harrison, C., Green, E., et al. (2018b). True Molecular Scale Visualization of Variable Clustering Properties of Ryanodine Receptors. Cell Rep. 22, 557–567. doi: 10.1016/j.celrep.2017.12.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Jayasinghe, I., Clowsley, A. H., de Langen, O., Sali, S. S., Crossman, D. J., and Soeller, C. (2018a). Shining New Light on the Structural Determinants of Cardiac Couplon Function: insights From Ten Years of Nanoscale Microscopy. Front. Physiol. 9:1472. doi: 10.3389/fphys.2018.01472

PubMed Abstract | CrossRef Full Text | Google Scholar

Kekenes-Huskey, P., Cheng, Y., Hake, J., Sachse, F., Bridge, J., Holst, M., et al. (2012). Modeling Effects of L-Type Ca2+ Current and Na+-Ca2+ Exchanger on Ca2+ Trigger Flux in Rabbit Myocytes with Realistic T-Tubule Geometries. Front. Physiol. 3:351. doi: 10.3389/fphys.2012.00351

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerr, R. A., Bartol, T. M., Kaminsky, B., Dittrich, M., Chang, J.-C. J., Baden, S. B., et al. (2008). Fast Monte Carlo Simulation Methods for Biological Reaction-Diffusion Systems in Solution and on Surfaces. SIAM J. Sci. Comput. 30, 3126–3149. doi: 10.1137/070692017

PubMed Abstract | CrossRef Full Text | Google Scholar

Kirk, M. M., Izu, L. T., Chen-Izu, Y., McCulle, S. L., Wier, W. G., Balke, C. W., et al. (2003). Role of the Transverse-Axial Tubule System in Generating Calcium Sparks and Calcium Transients in Rat Atrial Myocytes. J. Physiol. 547, 441–451. doi: 10.1113/jphysiol.2002.034355

PubMed Abstract | CrossRef Full Text | Google Scholar

Ko, C. Y., Liu, M. B., Song, Z., Qu, Z., and Weiss, J. N. (2017). Multiscale Determinants of Delayed Afterdepolarization Amplitude in Cardiac Tissue. Biophys. J. 112, 1949–1961. doi: 10.1016/j.bpj.2017.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Koivumäki, J. T., Korhonen, T., and Tavi, P. (2011). Impact of Sarcoplasmic Reticulum Calcium Release on Calcium Dynamics and Action Potential Morphology in Human Atrial Myocytes: a Computational Study. PLoS Comput. Biol. 7:e1001067. doi: 10.1371/journal.pcbi.1001067

PubMed Abstract | CrossRef Full Text | Google Scholar

Laver, D. R., Kong, C. H. T., Imtiaz, M. S., and Cannell, M. B. (2013). Termination of calcium-induced calcium release by induction decay: an emergent property of stochastic channel gating and molecular scale architecture. J. Mol. Cell. Cardiol. 54, 98–100. doi: 10.1016/j.yjmcc.2012.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, M. B., de Lange, E., Garfinkel, A., Weiss, J. N., and Qu, Z. (2015). Delayed afterdepolarizations generate both triggers and a vulnerable substrate promoting reentry in cardiac tissue. Heart Rhythm 12, 2115–2124. doi: 10.1016/j.hrthm.2015.06.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, M. B., Ko, C. Y., Song, Z., Garfinkel, A., Weiss, J. N., and Qu, Z. (2016). A Dynamical Threshold for Cardiac Delayed Afterdepolarization-Mediated Triggered Activity. Biophys. J. 111, 2523–2533. doi: 10.1016/j.bpj.2016.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Louch, W. E., Hake, J., Jølle, G. F., Mørk, H. K., Sjaastad, I., Lines, G. T., et al. (2010). Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes. Biophys. J. 99, 1377–1386. doi: 10.1016/j.bpj.2010.06.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, S., Michailova, A., Saucerman, J., Cheng, Y., Yu, Z., Kaiser, T., et al. (2009). Multi-Scale Modeling in Rodent Ventricular Myocytes: contributions of structural and functional heterogeneities to excitation-contraction coupling. IEEE Eng. Med. Biol. Mag. Q. Mag. Eng. Med. Biol. Soc. 28, 46–57. doi: 10.1109/MEMB.2009.931787

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyon, A. R., MacLeod, K. T., Zhang, Y., Garcia, E., Kanda, G. K., Lab, M. J., et al. (2009). Loss of T-tubules and other changes to surface topography in ventricular myocytes from failing human and rat heart. Proc. Natl. Acad. Sci. U. S. A. 106, 6854–6859. doi: 10.1073/pnas.0809777106

PubMed Abstract | CrossRef Full Text | Google Scholar

Macquaide, N., Tuan, H.-T. M., Hotta, J., Sempels, W., Lenaerts, I., Holemans, P., et al. (2015). Ryanodine receptor cluster fragmentation and redistribution in persistent atrial fibrillation enhance calcium release. Cardiovasc. Res. 108, 387–398. doi: 10.1093/cvr/cvv231

PubMed Abstract | CrossRef Full Text | Google Scholar

Maleckar, M. M., Edwards, A. G., Louch, W. E., and Lines, G. T. (2017). Studying dyadic structure–function relationships: a review of current modeling approaches and new insights into Ca2+ (mis)handling. Clin. Med. Insights Cardiol. 11:1179546817698602. doi: 10.1177/1179546817698602

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltsev, A. V., Maltsev, V. A., and Stern, M. D. (2017). Stabilization of diastolic calcium signal via calcium pump regulation of complex local calcium releases and transient decay in a computational model of cardiac pacemaker cell with individual release channels. PLoS Comput. Biol. 13:e1005675. doi: 10.1371/journal.pcbi.1005675

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltsev, V. A., and Lakatta, E. G. (2013). Numerical models based on a minimal set of sarcolemmal electrogenic proteins and an intracellular Ca2+ clock generate robust, flexible, and energy-efficient cardiac pacemaking. J. Mol. Cell. Cardiol. 59, 181–195. doi: 10.1016/j.yjmcc.2013.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchena, M., and Echebarria, B. (2018). Computational Model of Calcium Signaling in Cardiac Atrial Cells at the Submicron Scale. Front. Physiol. 9:1760. doi: 10.3389/fphys.2018.01760

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchena, M., and Echebarria, B. (2020). Influence of the tubular network on the characteristics of calcium transients in cardiac myocytes. PLoS One 15:e0231056. doi: 10.1371/journal.pone.0231056

PubMed Abstract | CrossRef Full Text | Google Scholar

Marx, S. O., Gaburjakova, J., Gaburjakova, M., Henrikson, C., Ondrias, K., and Marks, A. R. (2001). Coupled Gating Between Cardiac Calcium Release Channels (Ryanodine Receptors). Circ. Res. 88, 1151–1158. doi: 10.1161/hh1101.091268

PubMed Abstract | CrossRef Full Text | Google Scholar

Mesa, M. H., van den Brink, J., Louch, W. E., McCabe, K. J., and Rangamani, P. (2021). Nanoscale organization of ryanodine receptor distribution and phosphorylation pattern determines the dynamics of calcium sparks. bioRxiv [Preprint] doi: 10.1101/2021.10.19.465028

CrossRef Full Text | Google Scholar

Mori, Y., Fishman, G. I., and Peskin, C. S. (2008). Ephaptic conduction in a cardiac strand model with 3D electrodiffusion. Proc. Natl. Acad. Sci. U. S. A. 105, 6463–6468. doi: 10.1073/pnas.0801089105

PubMed Abstract | CrossRef Full Text | Google Scholar

Mukherjee, S., Thomas, N. L., and Williams, A. J. (2012). A mechanistic description of gating of the human cardiac ryanodine receptor in a regulated minimal environment. J. Gen. Physiol. 140, 139–158. doi: 10.1085/jgp.201110706

PubMed Abstract | CrossRef Full Text | Google Scholar

Muszkiewicz, A., Britton, O. J., Gemmell, P., Passini, E., Sánchez, C., Zhou, X., et al. (2016). Variability in cardiac electrophysiology: using experimentally-calibrated populations of models to move beyond the single virtual physiological human paradigm. Prog. Biophys. Mol. Biol. 120, 115–127. doi: 10.1016/j.pbiomolbio.2015.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Nivala, M., de Lange, E., Rovetti, R., and Qu, Z. (2012a). Computational Modeling and Numerical Methods for Spatiotemporal Calcium Cycling in Ventricular Myocytes. Front. Physiol. 3:114. doi: 10.3389/fphys.2012.00114

PubMed Abstract | CrossRef Full Text | Google Scholar

Nivala, M., Ko, C. Y., Nivala, M., Weiss, J. N., and Qu, Z. (2012b). Criticality in intracellular calcium signaling in cardiac myocytes. Biophys. J. 102, 2433–2442. doi: 10.1016/j.bpj.2012.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Nivala, M., Song, Z., Weiss, J. N., and Qu, Z. (2015). T-tubule disruption promotes calcium alternans in failing ventricular myocytes: mechanistic insights from computational modeling. J. Mol. Cell. Cardiol. 79, 32–41. doi: 10.1016/j.yjmcc.2014.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Okada, J., Sugiura, S., Nishimura, S., and Hisada, T. (2005). Three-dimensional simulation of calcium waves and contraction in cardiomyocytes using the finite element method. Am. J. Physiol.-Cell Physiol. 288, C510–C522. doi: 10.1152/ajpcell.00261.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandey, V., Xie, L.-H., Qu, Z., and Song, Z. (2021). Mitochondrial Contributions in the Genesis of Delayed Afterdepolarizations in Ventricular Myocytes. Front. Physiol. 12:744023. doi: 10.3389/fphys.2021.744023

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, S. H., Kim, A., An, J., Cho, H. S., and Kang, T. M. (2020). Nanoscale imaging of rat atrial myocytes by scanning ion conductance microscopy reveals heterogeneity of T-tubule openings and ultrastructure of the cell membrane. Korean J. Physiol. Pharmacol. 24, 529–543. doi: 10.4196/kjpp.2020.24.6.529

PubMed Abstract | CrossRef Full Text | Google Scholar

Passini, E., Britton, O. J., Lu, H. R., Rohrbacher, J., Hermans, A. N., Gallacher, D. J., et al. (2017). Human In Silico Drug Trials Demonstrate Higher Accuracy than Animal Models in Predicting Clinical Pro-Arrhythmic Cardiotoxicity. Front. Physiol. 8:668. doi: 10.3389/fphys.2017.00668

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinali, C., Bennett, H., Davenport, J. B., Trafford, A. W., and Kitmitto, A. (2013). Three-dimensional reconstruction of cardiac sarcoplasmic reticulum reveals a continuous network linking transverse-tubules: this organization is perturbed in heart failure. Circ. Res. 113, 1219–1230. doi: 10.1161/CIRCRESAHA.113.301348

PubMed Abstract | CrossRef Full Text | Google Scholar

Pods, J., Schönke, J., and Bastian, P. (2013). Electrodiffusion Models of Neurons and Extracellular Space Using the Poisson-Nernst-Planck Equations—Numerical Simulation of the Intra- and Extracellular Potential for an Axon Model. Biophys. J. 105, 242–254. doi: 10.1016/j.bpj.2013.05.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Protasi, F., Franzini-Armstrong, C., and Allen, P. D. (1998). Role of Ryanodine Receptors in the Assembly of Calcium Release Units in Skeletal Muscle. J. Cell Biol. 140, 831–842. doi: 10.1083/jcb.140.4.831

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, Z., Hu, G., Garfinkel, A., and Weiss, J. N. (2014). Nonlinear and stochastic dynamics in the heart. Phys. Rep. 543, 61–162. doi: 10.1016/j.physrep.2014.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, Z., Liu, M. B., and Nivala, M. (2016). A unified theory of calcium alternans in ventricular myocytes. Sci. Rep. 6:35625. doi: 10.1038/srep35625

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, Z., Shiferaw, Y., and Weiss, J. N. (2007). Nonlinear dynamics of cardiac excitation-contraction coupling: an iterated map study. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 75:011927. doi: 10.1103/PhysRevE.75.011927

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajagopal, V., Bass, G., Ghosh, S., Hunt, H., Walker, C., Hanssen, E., et al. (2018). Creating a Structurally Realistic Finite Element Geometric Model of a Cardiomyocyte to Study the Role of Cellular Architecture in Cardiomyocyte Systems Biology. JoVE J. Vis. Exp. 134:e56817. doi: 10.3791/56817

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajagopal, V., Bass, G., Walker, C. G., Crossman, D. J., Petzer, A., Hickey, A., et al. (2015). Examination of the Effects of Heterogeneous Organization of RyR Clusters, Myofibrils and Mitochondria on Ca2+ Release Patterns in Cardiomyocytes. PLoS Comput. Biol. 11:e1004417. doi: 10.1371/journal.pcbi.1004417

PubMed Abstract | CrossRef Full Text | Google Scholar

Rathinam, M., Petzold, L. R., Cao, Y., and Gillespie, D. T. (2003). Stiffness in stochastic chemically reacting systems: the implicit tau-leaping method. J. Chem. Phys. 119, 12784–12794. doi: 10.1063/1.1627296

CrossRef Full Text | Google Scholar

Restrepo, J. G., Weiss, J. N., and Karma, A. (2008). Calsequestrin-mediated mechanism for cellular calcium transient alternans. Biophys. J. 95, 3767–3789. doi: 10.1529/biophysj.108.130419

PubMed Abstract | CrossRef Full Text | Google Scholar

Richards, M. A., Clarke, J. D., Saravanan, P., Voigt, N., Dobrev, D., Eisner, D. A., et al. (2011). Transverse tubules are a common feature in large mammalian atrial myocytes including human. Am. J. Physiol. - Heart Circ. Physiol. 301, H1996–H2005. doi: 10.1152/ajpheart.00284.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Romero, L., Alvarez-Lacalle, E., and Shiferaw, Y. (2019). Stochastic coupled map model of subcellular calcium cycling in cardiac cells. Chaos Interdiscip. J. Nonlin. Sci. 29:023125. doi: 10.1063/1.5063462

CrossRef Full Text | Google Scholar

Rovetti, R., Das, K. K., Garfinkel, A., and Shiferaw, Y. (2007). Macroscopic consequences of calcium signaling in microdomains: a first-passage-time approach. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 76:051920. doi: 10.1103/PhysRevE.76.051920

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, D., Bartos, D. C., Ginsburg, K. S., and Bers, D. M. (2014). Depolarization of Cardiac Membrane Potential Synchronizes Calcium Sparks and Waves in Tissue. Biophys. J. 107, 1313–1317. doi: 10.1016/j.bpj.2014.07.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, D., and Bers, D. M. (2011). How Does Stochastic Ryanodine Receptor-Mediated Ca Leak Fail to Initiate a Ca Spark? Biophys. J. 101, 2370–2379. doi: 10.1016/j.bpj.2011.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, D., Bers, D. M., and Shiferaw, Y. (2013). Formation of spatially discordant alternans due to fluctuations and diffusion of calcium. PLoS One 8:e85365. doi: 10.1371/journal.pone.0085365

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, D., Despa, S., and Bers, D. M. (2012). Can the Sodium-Calcium Exchanger Initiate or Suppress Calcium Sparks in Cardiac Myocytes? Biophys. J. 102, L31–L33. doi: 10.1016/j.bpj.2012.03.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, D., Shannon, T. R., and Bers, D. M. (2016). Sarcoplasmic Reticulum Structure and Functional Properties that Promote Long-Lasting Calcium Sparks. Biophys. J. 110, 382–390. doi: 10.1016/j.bpj.2015.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, D., Uchinoumi, H., and Bers, D. M. (2021). Increasing SERCA function promotes initiation of calcium sparks and breakup of calcium waves. J. Physiol. 599, 3267–3278. doi: 10.1113/JP281579

PubMed Abstract | CrossRef Full Text | Google Scholar

Scardigli, M., Ferrantini, C., Crocini, C., Pavone, F. S., and Sacconi, L. (2018). Interplay Between Sub-Cellular Alterations of Calcium Release and T-Tubular Defects in Cardiac Diseases. Front. Physiol. 9:1474. doi: 10.3389/fphys.2018.01474

PubMed Abstract | CrossRef Full Text | Google Scholar

Scriven, D. R. L., Asghari, P., and Moore, E. D. W. (2013). Microarchitecture of the dyad. Cardiovasc. Res. 98, 169–176. doi: 10.1093/cvr/cvt025

PubMed Abstract | CrossRef Full Text | Google Scholar

Setterberg, I. E., Le, C., Frisk, M., Li, J., and Louch, W. E. (2021). The Physiology and Pathophysiology of T-Tubules in the Heart. Front. Physiol. 12:718404. doi: 10.3389/fphys.2021.718404

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, S. J., Aistrup, G. L., Gupta, D. K., O’Toole, M. J., Nahhas, A. F., Schuster, D., et al. (2014). Ultrastructural and cellular basis for the development of abnormal myocardial mechanics during the transition from hypertension to heart failure. Am. J. Physiol. Heart Circ. Physiol. 306, H88–H100. doi: 10.1152/ajpheart.00642.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, T. R., Wang, F., Puglisi, J., Weber, C., and Bers, D. M. (2004). A Mathematical Treatment of Integrated Ca Dynamics within the Ventricular Myocyte. Biophys. J. 87, 3351–3371. doi: 10.1529/biophysj.104.047449

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheard, T. M. D., Hurley, M. E., Colyer, J., White, E., Norman, R., Pervolaraki, E., et al. (2019). Three-Dimensional and Chemical Mapping of Intracellular Signaling Nanodomains in Health and Disease with Enhanced Expansion Microscopy. ACS Nano 13, 2143–2157. doi: 10.1021/acsnano.8b08742

PubMed Abstract | CrossRef Full Text | Google Scholar

Shiferaw, Y., Aistrup, G. L., Louch, W. E., and Wasserstrom, J. A. (2020). Remodeling Promotes Proarrhythmic Disruption of Calcium Homeostasis in Failing Atrial Myocytes. Biophys. J. 118, 476–491. doi: 10.1016/j.bpj.2019.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Shiferaw, Y., Aistrup, G. L., and Wasserstrom, J. A. (2017). Mechanism for Triggered Waves in Atrial Myocytes. Biophys. J. 113, 656–670. doi: 10.1016/j.bpj.2017.06.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Shiferaw, Y., Aistrup, G. L., and Wasserstrom, J. A. (2018). Synchronization of Triggered Waves in Atrial Tissue. Biophys. J. 115, 1130–1141. doi: 10.1016/j.bpj.2018.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Shiferaw, Y., Watanabe, M. A., Garfinkel, A., Weiss, J. N., and Karma, A. (2003). Model of intracellular calcium cycling in ventricular myocytes. Biophys. J. 85, 3666–3686. doi: 10.1016/S0006-3495(03)74784-5

CrossRef Full Text | Google Scholar

Singh, J. K., Barsegyan, V., Bassi, N., Marszalec, W., Tai, S., Mothkur, S., et al. (2017). T-tubule remodeling and increased heterogeneity of calcium release during the progression to heart failure in intact rat ventricle. Physiol. Rep. 5:e13540. doi: 10.14814/phy2.13540

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobie, E. A., Dilly, K. W., dos Santos Cruz, J., Lederer, W. J., and Jafri, M. S. (2002). Termination of cardiac Ca(2+) sparks: an investigative mathematical model of calcium-induced calcium release. Biophys. J. 83, 59–78. doi: 10.1016/s0006-3495(02)75149-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobie, E. A., Song, L.-S., and Lederer, W. (2005). Local recovery of Ca2+ release in rat ventricular myocytes. J. Physiol. 565, 441–447. doi: 10.1113/jphysiol.2005.086496

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Z., Karma, A., Weiss, J. N., and Qu, Z. (2016). Long-Lasting Sparks: multi-Metastability and Release Competition in the Calcium Release Unit Network. PLoS Comput. Biol. 12:e1004671. doi: 10.1371/journal.pcbi.1004671

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Z., Ko, C. Y., Nivala, M., Weiss, J. N., and Qu, Z. (2015). Calcium-voltage coupling in the genesis of early and delayed afterdepolarizations in cardiac myocytes. Biophys. J. 108, 1908–1921. doi: 10.1016/j.bpj.2015.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Z., Liu, M. B., and Qu, Z. (2018). Transverse tubular network structures in the genesis of intracellular calcium alternans and triggered activity in cardiac cells. J. Mol. Cell. Cardiol. 114, 288–299. doi: 10.1016/j.yjmcc.2017.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Z., Qu, Z., and Karma, A. (2017). Stochastic initiation and termination of calcium-mediated triggered activity in cardiac myocytes. Proc. Natl. Acad. Sci. U. S. A. 114, E270–E279. doi: 10.1073/pnas.1614051114

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Z., Xie, L.-H., Weiss, J. N., and Qu, Z. (2019). A Spatiotemporal Ventricular Myocyte Model Incorporating Mitochondrial Calcium Cycling. Biophys. J. 117, 2349–2360. doi: 10.1016/j.bpj.2019.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Stern, M. D., Song, L.-S., Cheng, H., Sham, J. S. K., Yang, H. T., Boheler, K. R., et al. (1999). Local Control Models of Cardiac Excitation–Contraction Coupling. J. Gen. Physiol. 113, 469–489.

Google Scholar

Sutanto, H., Lyon, A., Lumens, J., Schotten, U., Dobrev, D., and Heijman, J. (2020). Cardiomyocyte calcium handling in health and disease: insights from in vitro and in silico studies. Prog. Biophys. Mol. Biol. 157, 54–75. doi: 10.1016/j.pbiomolbio.2020.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Sutanto, H., van Sloun, B., Schönleitner, P., van Zandvoort, M. A. M. J., Antoons, G., and Heijman, J. (2018). The Subcellular Distribution of Ryanodine Receptors and L-Type Ca2+ Channels Modulates Ca2+-Transient Properties and Spontaneous Ca2+-Release Events in Atrial Cardiomyocytes. Front. Physiol. 9:1108. doi: 10.3389/fphys.2018.01108

PubMed Abstract | CrossRef Full Text | Google Scholar

Theakston, E., Walker, C., O’Sullivan, M., and Rajagopal, V. (2010). Stochastic modelling of cardiac cell structure. Annu. Int. Conf. IEEE Eng. Med. Biol. Soc. 2010, 3257–3260. doi: 10.1109/IEMBS.2010.5627229

PubMed Abstract | CrossRef Full Text | Google Scholar

Trafford, A. W., Clarke, J. D., Richards, M. A., Eisner, D. A., and Dibb, K. M. (2013). Calcium signalling microdomains and the t-tubular system in atrial mycoytes: potential roles in cardiac disease and arrhythmias. Cardiovasc. Res. 98, 192–203. doi: 10.1093/cvr/cvt018

PubMed Abstract | CrossRef Full Text | Google Scholar

Vagos, M., van Herck, I. G. M., Sundnes, J., Arevalo, H. J., Edwards, A. G., and Koivumäki, J. T. (2018). Computational Modeling of Electrophysiology and Pharmacotherapy of Atrial Fibrillation: recent Advances and Future Challenges. Front. Physiol. 9:1221. doi: 10.3389/fphys.2018.01221

PubMed Abstract | CrossRef Full Text | Google Scholar

Vagos, M. R., Arevalo, H., Heijman, J., Schotten, U., and Sundnes, J. (2020). A Novel Computational Model of the Rabbit Atrial Cardiomyocyte With Spatial Calcium Dynamics. Front. Physiol. 11:556156. doi: 10.3389/fphys.2020.556156

PubMed Abstract | CrossRef Full Text | Google Scholar

Vagos, M. R., Arevalo, H., Heijman, J., Schotten, U., and Sundnes, J. (2021). A Computational Study of the Effects of Tachycardia-Induced Remodeling on Calcium Wave Propagation in Rabbit Atrial Myocytes. Front. Physiol. 12:651428. doi: 10.3389/fphys.2021.651428

PubMed Abstract | CrossRef Full Text | Google Scholar

Voigt, N., Heijman, J., Wang, Q., Chiang, D. Y., Li, N., Karck, M., et al. (2014). Cellular and molecular mechanisms of atrial arrhythmogenesis in patients with paroxysmal atrial fibrillation. Circulation 129, 145–156. doi: 10.1161/CIRCULATIONAHA.113.006641

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, E., Lauterbach, M. A., Kohl, T., Westphal, V., Williams, G. S. B., Steinbrecher, J. H., et al. (2012). Stimulated Emission Depletion Live-Cell Super-Resolution Imaging Shows Proliferative Remodeling of T-Tubule Membrane Structures After Myocardial Infarction. Circ. Res. 111, 402–414. doi: 10.1161/CIRCRESAHA.112.274530

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, M. A., Gurev, V., Rice, J. J., Greenstein, J. L., and Winslow, R. L. (2017). Estimating the probabilities of rare arrhythmic events in multiscale computational models of cardiac cells and tissue. PLoS Comput. Biol. 13:e1005783. doi: 10.1371/journal.pcbi.1005783

PubMed Abstract | CrossRef Full Text | Google Scholar

Walker, M. A., Williams, G. S. B., Kohl, T., Lehnart, S. E., Jafri, M. S., Greenstein, J. L., et al. (2014). Superresolution Modeling of Calcium Release in the Heart. Biophys. J. 107, 3018–3029. doi: 10.1016/j.bpj.2014.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Weinberg, S. H., Hao, Y., Sobie, E. A., and Smith, G. D. (2015). Calcium homeostasis in a local/global whole cell model of permeabilized ventricular myocytes with a Langevin description of stochastic calcium release. Am. J. Physiol.-Heart Circ. Physiol. 308, H510–H523. doi: 10.1152/ajpheart.00296.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, J., Yao, J., Belke, D., Guo, W., Zhong, X., Sun, B., et al. (2021). Ca2+-CaM Dependent Inactivation of RyR2 Underlies Ca2+ Alternans in Intact Heart. Circ. Res. 128, e63–e83. doi: 10.1161/CIRCRESAHA.120.318429

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, G. S. B., Chikando, A. C., Tuan, H.-T. M., Sobie, E. A., Lederer, W. J., and Jafri, M. S. (2011). Dynamics of calcium sparks and calcium leak in the heart. Biophys. J. 101, 1287–1296. doi: 10.1016/j.bpj.2011.07.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, G. S. B., Huertas, M. A., Sobie, E. A., Jafri, M. S., and Smith, G. D. (2007). A Probability Density Approach to Modeling Local Control of Calcium-Induced Calcium Release in Cardiac Myocytes. Biophys. J. 92, 2311–2328. doi: 10.1529/biophysj.106.099861

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, G. S. B., Huertas, M. A., Sobie, E. A., Jafri, M. S., and Smith, G. D. (2008). Moment Closure for Local Control Models of Calcium-Induced Calcium Release in Cardiac Myocytes. Biophys. J. 95, 1689–1703. doi: 10.1529/biophysj.107.125948

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, A., Song, Z., Liu, H., Zhou, A., Shi, G., Wang, Q., et al. (2018). Mitochondrial Ca2+ Influx Contributes to Arrhythmic Risk in Nonischemic Cardiomyopathy. J. Am. Heart Assoc. 7:e007805. doi: 10.1161/JAHA.117.007805

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y., Sato, D., Garfinkel, A., Qu, Z., and Weiss, J. N. (2010). So little source, so much sink: requirements for afterdepolarizations to propagate in tissue. Biophys. J. 99, 1408–1415. doi: 10.1016/j.bpj.2010.06.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaniv, Y., Lakatta, E. G., and Maltsev, V. A. (2015). From two competing oscillators to one coupled-clock pacemaker cell system. Front. Physiol. 6:28. doi: 10.3389/fphys.2015.00028

PubMed Abstract | CrossRef Full Text | Google Scholar

Zahradníková, A., Zahradník, I., Györke, I., and Györke, S. (1999). Rapid Activation of the Cardiac Ryanodine Receptor by Submillisecond Calcium Stimuli. J. Gen. Physiol. 114, 787–798. doi: 10.1085/jgp.114.6.787

PubMed Abstract | CrossRef Full Text | Google Scholar