Arrhythmia mechanisms and spontaneous calcium release: Bi-directional coupling between re-entrant and focal excitation

Michael A. Colman

doi:10.1371/journal.pcbi.1007260

Abstract

Spontaneous sub-cellular calcium release events (SCRE) are conjectured to promote rapid arrhythmias associated with conditions such as heart failure and atrial fibrillation: they can underlie the emergence of spontaneous action potentials in single cells which can lead to arrhythmogenic triggers in tissue. The multi-scale mechanisms of the development of SCRE into arrhythmia triggers, and their dynamic interaction with the tissue substrate, remain elusive; rigorous and simultaneous study of dynamics from the nanometre to the centimetre scale is a major challenge. The aim of this study was to develop a computational approach to overcome this challenge and study potential bi-directional coupling between sub-cellular and tissue-scale arrhythmia phenomena. A framework comprising a hierarchy of computational models was developed, which includes detailed single-cell models describing spatio-temporal calcium dynamics in 3D, efficient non-spatial cell models, and both idealised and realistic tissue models. A phenomenological approach was implemented to reproduce SCRE morphology and variability in the efficient cell models, comprising the definition of analytical Spontaneous Release Functions (SRF) whose parameters may be randomly sampled from appropriate distributions in order to match either the 3D cell models or experimental data. Pro-arrhythmogenic pacing protocols were applied to initiate re-entry and promote calcium overload, leading to the emergence of SCRE. The SRF accurately reproduced the dynamics of SCRE and its dependence on environment variables under multiple different conditions. Sustained re-entrant excitation promoted calcium overload, and led to the emergence of focal excitations after termination. A purely functional mechanism of re-entry and focal activity localisation was demonstrated, related to the unexcited spiral wave core. In conclusion, a novel approach has been developed to dynamically model SCRE at the tissue scale, which facilitates novel, detailed multi-scale mechanistic analysis. It was revealed that complex re-entrant excitation patterns and SCRE may be bi-directionally coupled, promoting novel mechanisms of arrhythmia perpetuation.

Author summary

A loss of the regular rhythm of the beating heart, called arrhythmia, can inhibit its pumping function and even lead to sudden death. Understanding the processes by which normal rhythm is interrupted presents a major yet critical research problem. One challenge is the inherent multi-scale dependence of the electrical activity of the heart: behaviour at the microscopic scales (single proteins) can propagate to the macroscopic (whole-heart). Simultaneously studying phenomena at both of these scales is difficult, if not impossible, to perform experimentally. Developing mathematical models of the heart in order to perform variable and controlled simulations of its electrical activity provides the possibility to undertake such multi-scale analysis. However, new approaches are required to simulate the potential role of “spontaneous calcium release”, which occurs at the sub-cellular scale, in triggering arrhythmia events at the whole-heart scale. This study develops an approach to perform such simulations and applies it to study the long-term interactions of these sub-cellular triggers with the complex electrical activity of the most dangerous arrhythmias–tachycardia and fibrillation. The new mechanistic insight has implications for both diagnosis and treatment of the associated disorders.

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

Editor: Patrick Boyle, University of Washington, UNITED STATES

Received: April 5, 2019; Accepted: July 9, 2019; Published: August 8, 2019

Copyright: © 2019 Michael A. Colman. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the manuscript and its Supporting Information files.

Funding: This work was supported by a Medical Research Council Strategic Skills Fellowship, https://mrc.ukri.org/ (MR/M014967/1). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Cardiovascular disease is one of the major healthcare problems faced by the developed world, with increasing prevalence associated with aging populations [1–3]. Improved understanding of the mechanisms underlying cardiac arrhythmias, a major component of cardiovascular diseases’ impact on morbidity and mortality, is vital to the effort to improve both quality and duration of life. Rapid arrhythmias, such as tachycardia and fibrillation, are associated with disordered and incomplete contraction, reducing cardiac output and potentially leading to sudden cardiac death. The underlying rapid and irregular electrical activation of cardiac tissue may be mediated by abnormal spontaneous pacing (focal ectopic activity; “arrhythmia triggers”), self-perpetuating re-entrant excitation (“arrhythmia substrate”), or a complex interplay between both mechanisms (trigger-substrate interactions) [4,5]. Management of these arrhythmias is typically challenging, often requiring invasive procedures such as implanted defibrillators or catheter ablation; even these interventions have limited success rates [6,7]. Understanding of the mechanisms underlying the genesis, perpetuation and recurrence of rapid arrhythmias will ultimately lead to the development of improved treatment strategies.

Whereas multiple experimental and simulation studies have investigated the tissue substrate for the emergence and perpetuation of re-entrant excitation [8,9], the multi-scale mechanisms of arrhythmia triggers, and their dynamic interaction with the tissue substrate, remain elusive [10]. Malfunction of the intracellular calcium (Ca²⁺) handling system has been implicated in the development of rapid arrhythmias, linking sub-cellular spontaneous Ca²⁺ release events (SCRE) to pro-arrhythmic triggers in single cell [11–15]. However, translation of these cellular data to assess the mechanisms and importance of SCRE in tissue-scale arrhythmia is a significant challenge [10], namely because of the dependence of SCRE on stochastic fluctuations at the sub-cellular, nanometre scale which must propagate to the whole-organ, centimetre scale. These complex multi-scale mechanisms are discussed in detail in previous reviews [11,12,16], and are briefly summarised below.

A feedback mechanism emerging from the structure-function relationships underlying cardiac cellular excitation-contraction coupling presents a potential pro-arrhythmogenic pathway for the propagation of random state-transitions at the microscopic scale (sub-cellular) to the macroscopic (whole-organ): (i) spontaneous opening of the ryanodine receptors (RyRs), which control release of Ca²⁺ from the Sarcoplasmic Reticulum (SR) into the intracellular space, can trigger spontaneous Ca²⁺ sparks in restricted nanodomains called dyads; (ii) the sub-cellular spatial distribution of dyads presents a substrate for the propagation of Ca²⁺ sparks as a whole-cell event; (iii) this spontaneous whole-cell Ca²⁺ transient can lead to cellular delayed-after-depolarisations (DADs) and triggered action potentials (TA) through activation of the sodium-calcium exchanger (NCX); (iv) TA may propagate in tissue as ectopic focal excitation, potentially leading to the initiation of transient or sustained arrhythmia. The dynamics of Ca²⁺ homeostasis, also strongly dependent on the refilling of the SR through the SR-Ca²⁺ pump (SERCA), is critical for the emergence of SCRE: cardiac cells typically exhibit a minimum SR-Ca²⁺ load threshold above which whole-cell SCRE occur. The cytosolic Ca²⁺ concentration and kinetics of the RyRs (which may be more prone to spontaneous release in disease [14,17]) also strongly affect SCRE dynamics and vulnerability, which can be reflected in adaptations of this SR-Ca²⁺ load dependence (wherein lower thresholds correspond to larger vulnerability to SCRE).

Computational modelling provides a viable approach for detailed simultaneous multi-scale evaluation of cardiac arrhythmia mechanisms. Nevertheless, simulating SCRE in tissue-scale models is non-trivial due to the contrasting computational requirements of models at these different scales: single-cell models capable of reproducing SCRE are computationally intensive and unsuitable for simulation of the thousands or millions of coupled cells comprising tissue models appropriate for studying arrhythmia mechanisms. New approaches must therefore be developed to realise this goal. Only a few previous studies have attempted such modelling, for example investigating: (i) the minimum tissue substrate for the emergence of focal excitations resulting from non-stochastic DADs [18]; (ii) the emergence of focal excitation from stochastic SCRE [19–23] and its potential interaction with extracellular matrix remodelling [24]; (iii) SCRE as a mechanism for both triggered activity and conduction block [25]; and (iv) the potential complex considerations for pharmacological action on both triggers and substrate [22].

In a previous study, a phenomenological approach was introduced to model the synchronisation of cellular triggers in tissue [22]. The present study aims to refine this approach, and extend it to allow dynamic modelling of trigger-substrate interactions. The novel approach, comprising a hierarchy of computational models, is applied to: (i) directly translate behaviour observed in detailed cell models to the tissue-scale and (ii) study the dynamic interactions of SCRE and re-entrant excitation. The approach is also generalised to allow fully controllable investigations and integrate limited experimental data.

Methods

This section first outlines the computational framework and the hierarchy of cell and tissue models of which it comprises. Then, analysis of stochastic SCRE observed using the spatial single cell models is discussed in relation to the development of analytical waveform functions which approximate SCRE in the non-spatial cell- and multi-dimensional tissue-models, termed Spontaneous Release Functions (SRF). Multiple implementations are presented in which the distributions from which waveform parameters are sampled are determined: (i) by directly controlled user inputs (“Direct Control” model); (ii) dynamically to reproduce the behaviour of the spatial cell model under multiple conditions (“Dynamic Fit” model); or (iii) dynamically based on user-controlled inputs (“General Dynamic” model). Finally, simulation protocols for studying different arrhythmia dynamics are described.

The computational framework

The developed framework consists of a hierarchy of computational models (Fig 1):

The microscopic or 3D spatio-temporal cell model is the baseline model for the framework. It accounts for spatially distributed dyads, Ca²⁺ diffusion within the cell, and stochastic state transitions in the L-type Ca²⁺ channel (LTCC) and RyR models (Fig 1Ai and 1B). This model was used to study SCRE at the cellular scale.
The deterministic, 0D or non-spatial cell model is derived from the microscopic model. It does not contain a distributed Ca²⁺ handling structure or account for stochastic state transitions (Fig 1Aii and 1B). This manuscript focuses on the derivation and introduction of analytical Spontaneous Release Functions (SRF) into 0D models to reproduce behaviour from the 3D spatio-temporal cell model.
The Tissue model refers to coupled 0D cell models in either idealised 2D sheets or 3D reconstructions of cardiac anatomy (Fig 1C).

Download:

Fig 1. Components of the multi-scale computational framework.

A–Schematic of single cell Ca²⁺ handling and ion current models. (i)– 3D, microscopic Ca²⁺ handling model, illustrating the 3D grid of calcium release units (CRUs; upper panel) and the compartments and Ca²⁺ fluxes within a single CRU (lower panel). Labelled are the dyadic cleft space (DS), sub-space (SS), bulk cytosolic space (CYTO), network and junctional SR spaces (NSR, JSR), and a T-tubule (TT); fluxes through the LTCCs (J_CaL), RyRs (J_rel), NCX (J_NaCa) and SERCA (J_up) are illustrated according to the key; double-headed black arrows indicate transfer between compartments; double-headed red arrows indicate diffusion between neighbouring CRUs. (ii)– 0D, non-spatial cell model, illustrating the same fluxes as in (i) but without inter-CRU diffusion. Global ion currents are illustrated along the membrane (which apply to both models). B–Whole-cell voltage (i) and calcium transient (ii) of the different ion-current models used in the present study, showing the hybrid-minimal model (left; minimal), O’Hara et al., human ventricular model [26] (middle; ORd) and Colman et al., human atrial model [27] (right; Col 2013). C–Tissue models, showing schematic of a 2D sheet model (i), and the 3D anatomical reconstructions: (ii)—human ventricular wedge [28]; (iii) whole canine ventricle [29]; and (iv) whole human atria [27,30–32].

https://doi.org/10.1371/journal.pcbi.1007260.g001

The microscopic, 3D spatio-temporal Ca²⁺ handling model.

The microscopic model, a simplification of a previously presented and structurally detailed model [33], consists of an idealised 3D intracellular Ca²⁺ system coupled to a point-source V_m and ion-current model (Fig 1Ai and 1B). The model is capable of reproducing normal and abnormal Ca²⁺ dynamics [33] and captures the general properties of experimentally observed SCRE [11,14,34,35] in-line with the most recent modelling studies [19,25,36]. The intracellular model incorporates a 3D grid of CRUs, each containing five compartments (three intracellular compartments comprising the bulk-cytoplasm, sub-space, and restricted dyadic cleft, as well as the network and junctional SR spaces; Fig 1Ai). The bulk cytosolic, sub-space and network SR compartments are coupled to neighbouring CRUs. Each dyad contains ~10 LTCCs and ~100 RyRs (varies between models and can be heterogeneous within a cell), described by stochastic differential equations [33]. Full details and model equations are provided in the S1 Text (Model Description). Most relevant for this study is the reaction for the dyadic cleft, Φ_ds, which contains intracellular Ca²⁺ release, J_rel. For a single dyad n, J_rel is described by: (1)

Where g_RyR is the maximal flux factor for Ca²⁺ release, ⁿv_ds is the volume of the dyadic cleft space, and ⁿN_{RyR_O} is the summed number of open RyR channels in dyad n, determined through the stochastic Monte-Carlo method. N_{RyR_O} is therefore the primary variable which describes the dynamics of SCRE, and thus the variable used to integrate spatial and non-spatial cell models.

The deterministic, 0D model.

The deterministic model structure is identical to a single CRU of the microscopic model and of the same form as the majority of contemporary cardiac cell models based on the Hogkin-Huxley approach (Fig 1Aii). Both the RyR and LTCC models are described by ordinary differential equations and solved through the deterministic forward-Euler method. N_{RyR_O} now represents the whole-cell average of the equivalent dyad variable in the 3D cell model. Modifications to the RyR kinetics model were required due to the poor recapitulation of whole cell Ca²⁺-induced-Ca²⁺-release (CICR) using deterministic approaches–see S1 Text (Model Description) for details.

Action potential and tissue models.

For the purpose of demonstrating the general potential of the developed framework for integration with contemporary cell models, multiple ionic models describing the non-Ca²⁺ dependent membrane currents were integrated with the Ca²⁺ handling system (Fig 1B): the model was integrated with simplified versions of the O’Hara et al., 2011 human ventricular AP model [26] and Colman et al., 2013 human atrial AP model [27,37]. Furthermore, a hybrid-minimal model (comprising of a minimal setup suitable for coupling with physiological Ca²⁺ currents) was developed which describes human AP morphology in the three transmural cell types of the ventricles as well as the atrial myocardium. This cell-type heterogeneity pertained only to differences in the ionic model and not the intracellular Ca²⁺ handling model. The minimal model is loosely based on a combination of the two biophysically detailed models, but with fewer components not corresponding directly to specific physiological ion currents, and designed with the motivation for full controllability of AP morphology.

Idealised 2D sheet models consist of a 2D array of coupled cells in isotropic medium (Fig 1Ci), including both homogeneous sheets of either the ventricular epicardial layer or the right atrial wall, and a model of the transmural heterogeneity in the ventricular wall (simple ratio of 1:1:1 ENDO:M:EPI cells). 3D models consist of a wedge reconstruction of the human ventricular wall (Fig 1Cii, [28]), a reconstruction of the whole canine ventricle [29] (Fig 1Ciii), and a reconstruction of the whole human atria [27,31] (Fig 1Civ). Further details of the tissue models are presented in S1 Text (Model Description); all tissue models were electrically homogeneous except for the 2D ventricular transmural sheet and human ventricular wedge models, which implemented ENDO, M and EPI cells.

Pro-arrhythmic conditions

To induce prominent full-cell release events through different cellular conditions, representative (but non-specific) models were included for isoprenaline (ISO, sympathetic response which enhances CICR) and two types of pro-SCRE general disease remodelling mimicking features observed in conditions such as AF and HF (e.g., [14,38]): (i) SERCA was up-regulated and NCX was down-regulated (R_SERCA/NCX); (ii) the SR-Ca²⁺ threshold for release was lowered through increased inter-CRU coupling (R_CRU-CRU). R_SERCA/NCX also involved remodelling of the ion currents, including a reduction of I_K1, to provide a different AP environment coupled to the Ca²⁺ handling model. Details are provided in S1 Text (Model Description). The purpose of these models, which are not biophysically detailed or representative of specific regulation or diseased conditions, was to induce pro-arrhythmic behaviour which is not observed under control conditions, and test and demonstrate the ability of the 0D approximations to accurately capture variable cellular conditions and translate these to the tissue-scale.

Partial Ca²⁺-clamp protocol

A partial Ca²⁺-clamp ladder protocol (Fig 2) was implemented to both derive and validate the SRF in the 0D model: At each step (total duration 2s), the intracellular- and SR-Ca²⁺ concentrations were initially clamped to specified values, with the SR-Ca²⁺ clamped concentration incrementing for each successive step (Fig 2Aii). When spontaneous intracellular release occurs, the Ca²⁺ concentrations were allowed to dynamically evolve (clamp constraint removed). The membrane potential was allowed to evolve during this protocol, and the conductances of I_Na and I_CaL were set to zero to prevent excitation and resulting interruption of the SCRE by CICR. This protocol illustrates the variety of SCRE and their underlying spatio-temporal dynamics which must be captured in the 0D model (Fig 2B; S1 Video).

Download:

Fig 2. Illustration of Ca²⁺ clamp protocol.

A–Ca²⁺ clamp protocol illustrated for 9 steps of SR-Ca²⁺, showing traces for: (i) proportion open RyR; (ii) intracellular- (purple) and SR- (blue) Ca²⁺ concentration. B–Snapshots of the spatio-temporal Ca²⁺ dynamics at different SR-Ca²⁺ concentrations, showing: (i) non-propagating sparks; (ii) slow Ca²⁺ wave; (iii) multiple and rapid Ca²⁺ waves; the time range for the snapshots is shown in the square brackets. The data shown are clipped to the first of the two seconds associated with each clamp step in order to clearly visualise the waveforms.

https://doi.org/10.1371/journal.pcbi.1007260.g002

Derivation of the Spontaneous Release Functions (SRF)

The phenomenological approach involved the development of SRF which describe whole-cell RyR dynamics associated with SCRE observed in the 3D cell model, based on an extension of previous work [21,22]. These functions are waveforms which approximate the range of morphologies for the time series of the whole-cell N_{RyR_O}/N_RyR observed during SCRE (Fig 3A).

Download:

Fig 3. Derivation of the Spontaneous Release Functions.

A–Traces of open RyR (N_{RyR_O}/N_RyR) associated with 250 simulations of SCRE at different SR-Ca²⁺ concentrations (low = 1125 mM, purple; and high = 1200 mM, blue; 100 traces for each condition shown) in the 3D cell model. Each trace represents an individual simulation, and two simulations are highlighted for each of the two SR-Ca²⁺ values (note the two blue traces overlap, but are separate simulations); A(inset)–spontaneous Ca²⁺ transients which correspond to the four highlighted RyR traces. B–Histogram describing the initiation time for SCRE and C–the corresponding cumulative frequency plots, associated with 250 simulations of each condition. Inset–example of fitting the cumulative frequency with two sigmoidal functions (F₁(t_i) and F₂(t_i), orange and green), separated at a specific point (t_{i_Sep}, F(_ti)|_{ti = ti_Sep}, red triangular marker). D–Examples of two types of waveform, corresponding to those highlighted in A (upper panel) and the SRF which approximate them (lower panel). Labelled are the parameters which fully describe the waveforms: the initiation time (t_i), peak time (t_p), final time (t_f), initiation time of spike during plateau (t_i^spike), peak open RyR (N_{RyR_O}^peak), plateau open RyR (N_{RyR_O}^plateau). E–Histogram illustrating the distribution of RyR waveform duration for the two conditions, with the median duration (MD) and width of distribution either side of the median (DW₁, DW₂) labelled for Distribution 1 (upper panel); relationship between the distribution widths and median (lower panel, points–data; lines–fit by Eqs (18 and 19)). F–Correlation between peak of open RyR and the duration, shown for the two conditions featured in the Figure (purple, blue) and all other simulations (yellow). Low amplitude SCRE occurring near the threshold SR-Ca²⁺ are shown in orange. Fit by Eq (12) is shown by the red line.

https://doi.org/10.1371/journal.pcbi.1007260.g003

The concept of these SRF was introduced in a previous work [22] in the context of a statically defined model (corresponding to the Direct Control approach, below). This required the parameters describing the SRF to be defined by the user for specific simulations and therefore is not suitable for investigating natural cellular responses to different conditions or coupled dynamic simulations. Here, the functions themselves are refined, a method to derive the function parameters dynamically from model variables is introduced, leading to behaviour directly in-line with the 3D cell models, and the approaches are generalised to be fully controllable and suitable for parameterisation to experimental data. For completeness and context, the SRF are described here in full.

The Spontaneous Release Functions (SRF).

Results from the Ca²⁺ clamp ladder protocol across all SR-Ca²⁺ values were used to derive the analytical formulations of the SRF describing the variability in spontaneous release N_{RyR_O} waveforms. Comparison of SCRE from 250 simulations at low and high SR-Ca²⁺ illustrates the range of waveforms observed (Fig 3A) and the corresponding probability density function and cumulative frequency of initiation times, t_i (Fig 3B and 3C).

The N_{RyR_O} waveforms can be grouped into two primary types: spike-like associated with short, large-amplitude release, and plateau-like associated with long, small-amplitude release (Fig 3D). For the spike-like morphology, the waveform can be well approximated with the simple function: (2) (3) (4) (5) (6) where t_i is the initiation time 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 of open proportion RyR (Fig 3D). The constants in Eqs (5 and 6) were obtained from best fits to the waveforms observed. The function for the plateau-like waveform (corresponding to durations longer than 300 ms) is derived from the same parameters: (7)

Where N_{RyR_O}^plateau is the amplitude of the plateau (Fig 3D). This equation assumes the same form for the spike occurring within the plateau, with its upstroke time being 50 ms and its decay time 35ms; t_i^spike (Fig 3D) therefore corresponds to t_p-50 (and its half maximal activation time t_p-25). The constants were similarly obtained from the best fits to the waveforms observed.

The waveform is therefore completely described by four or five parameters (for spike-like and plateau-like waveforms, respectively): (1) initiation time, t_i; (2) duration (λ = t_f-t_i); (3) peak time, t_p; and (4–5) amplitude (N_{RyR_O}^peak; N_{RyR_O}^plateau). In order to maintain physiological waveforms and randomly sample the parameter values from appropriate distributions, the nature of stochastic variation of these four parameters must be examined, as is discussed below.

Parameter distributions and inverse functions.

(1)—t_i: The probability density functions for the initiation time associated with each SR-Ca²⁺ value exhibit skewed distributions (Fig 3B). The cumulative frequency (Fig 3C) was well approximated by the use of two simple sigmoidal functions (Fig 3C-inset), maintaining the desire for restraint in the number of parameters and allowing simple and intuitive controllability: (8)

The distribution for t_i is therefore determined by four parameters: the initiation time corresponding to the point where the functions are separated (t_i,Sep); the cumulative frequency at this point (CF_ti,Sep = F(t_i)|_{ti = ti,Sep}), and the gradient parameter of each function (k_F1, k_F2—corresponding to the width of the distribution either side of t_i,Sep; Fig 4C-inset.)

Download:

Fig 4. SR-Ca²⁺ dependence of SRF distribution parameters.

A–The Dynamic Fit SRF parameters: summary data (points) and the fit from the relevant functions (lines) for the three Ca²⁺ handling conditions (purple–control; blue–CRU-CRU coupling enhancement; orange–SERCA upregulation and NCX downregulation model) against SR-Ca²⁺ for: i) probability of whole-cell SCRE; ii) initiation time corresponding to the separation point, t_i; iii) the normalised cumulative frequency at this point, CF_ti,Sep = F(_ti)|_{ti = ti,Sep}; iv) the k parameter for F₁(t_i) and (v) for F₂(t_i); (vi) the median duration, MD. B–The General Dynamic SRF parameters: i) Illustration of the curve for probability of SCRE and its relation to the user defined parameters (CaSR_thresdhold, CaSR_{P_range}) and derived parameters (CaSR_min). ii) Illustration of the function form describing t_i,Sep and median duration (MD) (purple line) and its relation to CaSR_min, CaSR_max, t_i,Sep^min, t_i,Sep^max, MD^min, MD^max, and how the range of the distributions varies with SR-Ca²⁺ (shaded regions) at two different non-linearity factors (H_width = 0.75, blue; = 2.5, orange).

https://doi.org/10.1371/journal.pcbi.1007260.g004

(2)–λ: The distributions for the duration are also non-normal, and well approximated by two sigmoidal functions describing the cumulative frequency for half of the data either side of the median duration (MD; Fig 3E): (9)

Where the widths (DW₁, DW₂, in ms) are a function of the MD (Fig 3E), given by: (10) (11)

Default parameters, fit to the waveforms observed in the control model, are given in S1 Text (Model Description). The duration distribution for each model condition (control and remodelling models) is therefore completely described by one variable, the median, MD, and the above model-specific parameters. Note that the widths (DW₁, DW₂) could also be specified directly for complete control over the variability in duration.

(3)—t_p: The timing of the peak varies approximately evenly within the duration of the wavefrom, occurring between 25 ms after the initiation (t_i) and 52 ms before the final time (t_f).

(4)—N_{RyR_O}^peak; N_{RyR_O}^plateau: The expectation value of the amplitude correlates strongly with duration, λ (Fig 3F): (12) (13)

With uniform variation around these expectation values. Note that it would not be appropriate to define the amplitude independently from the duration, due to the correlation between these two parameters corresponding to the total amount of Ca²⁺ released.

With this setup, therefore, all parameters of the waveform are derived from two primary waveform properties: the initiation time, t_i, and the duration, λ, which also determine the peak time and amplitude. The distributions describing the variability of these properties are entirely described by 5–7 parameters (t_i = f(t_{i_sep}, CF_{ti_sep}, k_F1, k_F2); λ = f(MD, DW₁, DW₂) where DW₁, DW₂ = f(MD) or manually specified). Producing a single SRF waveform therefore requires only that t_i and λ are sampled from the relevant distributions, by passing a random number between 0 and 1 into the corresponding inverse functions:

Initiation time, t_i; inverse function of Eq (8): (14)

Duration, λ; inverse function of Eq (9): (15)

And the derived variables, t_p and N_{RyR_O}^peak, N_{RyR_O}^plateau, can be randomly sampled from uniform distributions: (16) (17) (18)

The different implementations of the SRF therefore pertain to different methods to define the distributions from which the SRF parameters will be sampled, as well as the probability of release, P(SCR).

Sampling the parameters: The direct control SRF model.

For the simplest implementation of SCRE in 0D models, the user “Direct Control” model, SRF variability and morphology can therefore be described simply by explicitly defining the 5–7 parameters which describe the t_i and λ distributions as desired (for example, to fit a single dataset, or for controlled investigation), and then sampling from the resulting inverse functions on a beat-by-beat basis. The next section describes an approach to derive these parameters dynamically based on relevant Ca²⁺ handling environment variables.

Sampling the parameters: The dynamic fit SRF model.

The Dynamic Fit SRF model was derived through correlation of the parameters defining the t_i and λ distributions with the primary environment variable controlling SCRE, the SR-Ca²⁺ concentration (Fig 4A; S1 Text (Model Description)). Relation to this single variable was chosen for practicality and simplicity of the resulting equations. Other factors which affect SCRE dynamics (e.g. RyR dysfunction, remodelling of SERCA/NCX, changes in diastolic Ca²⁺) can then be indirectly accounted for by their impact on the SR-Ca²⁺ dependence; using these simple functions, the coefficients which fit the data from the two remodelling conditions (which both directly affected SCRE dynamics) could also be obtained (Fig 4A; S1 Text (Model Description)). Note that both remodelling models reduce the SR-Ca²⁺ threshold for SCRE compared to control. For R_CRU-CRU, this is a direct result of the reduced time constant of sub-space coupling; for R_SERCA/NCX, it is primarily the reduced Ca²⁺ efflux due to the lower expression of NCX which causes this shift, by more easily allowing spontaneous Ca²⁺ sparks to propagate to neighbouring CRUs.

Sampling the parameters: The general dynamic SRF model.

Whereas the Dynamic Fit SRF model described above ensures congruence between the 3D and 0D models, a General Dynamic SRF model was also derived which provides full control over the distributions determining the dynamic behaviour. The approach therefore involved defining the functions which correlate SR-Ca²⁺ with the SRF input parameters (P(SCR), t_{i_sep}, CF_ti,Sep, k_F1, k_F1, MD) based on intuitive controllable parameters (Fig 4B): (i)—The threshold for SCRE (CaSR_threshold); (ii)—The SR-Ca²⁺ range over which P(SCR) varies from 0 to 1 (CaSR_{P_range}); (iii)—The maximal SR-Ca²⁺ above which SCRE distributions converge (CaSR_max > CaSR_threshold + CaSR_{P_range}); (iv)—The minimum and maximum t_i,Sep and MD (t_i,Sep^min, t_i,Sep^max, MD^min, MD^max); (v)—The t_i and λ distribution widths at these extremes (t_i,width^min, t_i,width^max, λ_width^min, λ_width^max); And (vi)—the non-linearity of width variance (H_width); CF_ti,Sep was set to 0.4. Equations are given in S1 Text (Model Description).

Implementation with the 0D cell model.

These SRF were implemented within the 0D cell models using a simple algorithm (Fig 5): the input parameters (P(SCR), t_i,Sep, CF_ti,Sep, k_F1, k _F2, MD; see previous section: Parameter distributions) are defined at a certain time (see below) and then the waveform parameters (t_i, and λ, which in turn define t_p and N_{RyR_O}^peak; see previous section: Parameter distributions) are randomly sampled from the associated inverse functions (Eqs 14–18 above). When the SRF have been initiated (i.e., t_i < t < t_i + λ and N_{RyR_O}^SRF is > 0), the N_{RyR_O} in the cell model is set to N_{RyR_O}^SRF and the model thus evolves as if the equivalent SCRE was occurring in the 3D cell model.

Download:

Fig 5. Implementation algorithm.

Schematic of the algorithms used to integrate the SRF with the non-spatial cell models for the three implementations. Note that the Δ[Ca²⁺]_SR clause is only calculated if the SRF parameters have been set but SCRE has not yet been initiated, and RyR recovery is only calculated if the model has undergone an AP.

https://doi.org/10.1371/journal.pcbi.1007260.g005

For the simplest implementation of the SRF, the Direct Control model, this calculation is determined at the time of cellular excitation, setting t_i and λ based on the distributions defined by user input parameters (P(SCR), t_i,Sep, CF_{ti_sep}, k_F1, k_F2, MD) and five random numbers input into the corresponding inverse functions. The model will set SRF parameters based on these single distributions with every cellular excitation (note that the parameters may give an SCRE timing later than the next stimulated excitation, in which case it will not occur and be reset on the next excitation).

For the Dynamic Fit and General Dynamic models, the calculation is performed multiple times, dynamically determined during the simulation: Initially, after the cell model has undergone an AP and the RyR availability has recovered above a set threshold (inactivation state-occupancy reaches below 0.2), the SR-Ca²⁺ is input to first define the probability of release, P(SCR), from S1 Text (Model Description)–equations (121; 128), and then the remaining input parameters are defined from the SR-Ca²⁺ according to the appropriate functions [S1 Text (Model Description)–equations (122–126; 127–134)]. Once these parameters have been calculated, if the SR-Ca²⁺ concentration changes more than by a predefined value (0.01 mM) before the SCRE has been initiated, then the parameters are recalculated based on this new SR-Ca²⁺. SRF parameters are not calculated during excitation (i.e., when the RyR availability is low) but recalculated upon RyR recovery, to impose that the SRF do not interrupt deterministic CICR.

Simulation protocols

Rapid pacing protocol.

To study the emergence of triggered activity in single cell and tissue, a rapid pacing protocol was applied to promote SR-Ca²⁺ loading and the emergence of whole-cell SCRE: The models were paced to stable state under varying basic cycle lengths (BCL = 200–600 ms); the state variables were saved at the stable values and used as initial conditions to efficiently run any number of simulations of a short pacing period followed by a quiescent (non-applied pacing) period, within which the statistics of SCRE can be analysed. This protocol was applied to the 3D and 0D single cell models (implementing the Dynamic Fit SRF approach) to compare behaviour between these models and translate cell-type and model-condition differences to the tissue-scale.

Assessing the relationship between SR-Ca²⁺ and focal excitation.

The General Dynamic SRF implementation was used to assess the relationship between SR-Ca²⁺ and focal excitation in tissue, allowing simulation of both homogeneous and heterogeneous media with respect to SCRE properties (electrical heterogeneity, pertaining to the different cell-types, was not included in this analysis). A baseline SRF model was set with an SR-Ca²⁺ threshold for SCRE of 1.10 mM. Four alternative models were introduced with this SR-Ca²⁺ threshold varied by ± 0.05 and 0.1 mM. All other parameters are either identical or set relative to the threshold (S1 Text (Model Description)). Note that at a specific SR-Ca²⁺ concentration, implementations with a lower SR-Ca²⁺ threshold exhibit larger amplitude and more rapidly timed SCRE than those with a higher SR-Ca²⁺ threshold. A small amount of heterogeneity was first introduced by randomly assigning 60% of the cells to the baseline, and 10% to each of the four other distributions; a larger amount of heterogeneity was introduced by randomly assigning 20% of the cells to each of the five distributions. Note that this heterogeneity is designed to represent normal inter-cellular variability, rather than heterogeneous patches associated with disease remodelling, and is therefore implemented to preserve the mean for direct comparison to the homogeneous case.

Studying unidirectional conduction block.

Two potential mechanisms of SCRE-induced unidirectional excitation conduction block, which presents the possibility to degenerate into re-entry, have been shown in the studies of Liu et al., 2015 [25] and Campos et al., 2017 [20]: a well-established mechanism, in which focal excitation (in this instance as a result of SCRE) interacts with heterogeneity in refractory period and leads to asymmetric propagation patterns. And a novel mechanism by which DADs in combination with sodium channelopathy led to heterogeneous inactivation of the fast sodium current and conduction block. To demonstrate the capability of the SRF implementation to reproduce these dynamics, and to provide independent verification of the results in those previous studies, two protocols were used: (i) DAD mediated conduction block was simulated in a homogeneous 2D sheet of the epicardial layer of the ventricle in combination with impaired sodium current (conductance of I_Na scaled by 30%; 5 mV shift of the inactivation steady-state). These conditions were combined with SRF parameters which were below the threshold for focal excitation, but in which significant DADs were observed; (ii) Conduction block due to focal excitation and repolarisation heterogeneity was simulated in a 2D sheet model of the transmural cell-types in the ventricular wall (simple 1:1:1 ratio of EPI, M and ENDO cells), with the pre-pacing stimulus applied uniformly to the edge of the ENDO region; a small region was set to reduced I_K1 to promote localized excitation. The Direct Control implementation was used to impose early-timed and controllable SCRE.

Re-entry and SCRE interactions.

To study the long-term interactions of re-entry and SCRE, it was required to initiate a transient period of re-entrant excitation, and then for this to terminate to allow the emergence of SCRE. Re-entrant excitation was initiated using two different methods: the phase-distribution method [39,40], and cross-field stimulation. In general, a reduction in the conduction velocity (D reduced by up to half) and/or shortening of the AP (increase in I_Kr and/or I_Ks) was also employed to shorten the wavelength and promote sustained excitation. Different modifications were required by different cell models and conditions in order to lead to this complex excitation; it should be noted, however, that the aim of these studies was to determine how re-entry may promote SCRE, rather than to study the dynamics of re-entry itself. In many cases, re-entry self-terminated within 20s, representing natural self-termination of a transient re-entrant event. In other cases, a block of I_Kr or an increase in D (back to control values) was necessary to terminate re-entry which would otherwise sustain–the required termination mechanism also depended highly on the underlying cell model and conditions, and was not the focus of the simulations. These simulations were combined with the General Dynamic SRF model to investigate the possible influence of re-entry under different vulnerabilities to SCRE: simulations associated with multiple re-entrant conditions were analysed under multiple SRF parameter combinations. A rapid regular pacing protocol was also applied to match SR-Ca²⁺ loading as observed in different re-entrant cases, to study the potential differences between behaviour emerging from the two different loading mechanisms under comparable environmental conditions.

Results

In this section is first described validation of the Dynamic Fit SRF approach for the three Ca²⁺ system conditions before application of the framework, in order to: (i) demonstrate the mechanism of stochastically mediated SCRE-induced ectopic beats; (ii) examine the complexities underlying the SR-Ca²⁺ dependence of focal excitation; (iii) demonstrate two independent mechanisms of conduction block associated with SCRE; and (iv) study the potential multi-scale, pro-arrhythmic interactions between SCRE and re-entrant excitation.

Validation of the dynamic fit spontaneous release functions

The 0D model implementing the Dynamic Fit SRF model was first validated by comparison of whole-cell SCRE under Ca²⁺ clamp conditions with a second set of simulations (i.e., not those on which the model was derived) of the 3D cell model, for the control and remodelling conditions (S2 Text (Validation and Results)). These simulations highlight the strong agreement for waveform morphology and its variation, the distributions and their summary properties, and the differences between the control and remodelled conditions.

Secondly, the 0D and 3D models were compared across the range of ionic models (see Methods: Action Potential and Tissue Models) and remodelling/ISO conditions (see Methods: Pro-arrhythmic conditions). Examples of dynamics emerging in conditions close to the SR-Ca²⁺ threshold (Fig 6Aa) and high above it (Fig 6Ab) show good agreement between the 3D and 0D models (compare Fig 6Ai-ii with 6Aiv-v), importantly capturing the key features and differences between the conditions. The case for low SR-Ca²⁺ (Fig 6Aa) was intentionally selected as one of those in which the match was the poorest, in order to fully illustrate the quality of the approximation in an upfront manner; even here, the match is reasonable, and the main features of the behaviour are preserved (DADs with some TA; wide initiation time distribution). The distributions of initiation time and the probability of triggered APs across all conditions tested which resulted in notable SCRE also show good agreement (Fig 6B), confirming the ability of the 0D models to dynamically reproduce the SCRE of the 3D cell models and capture cellular and condition dependent behaviour. Note that the differences between the cell-types and conditions is not of primary interest, but, rather, it is the match between 3D and 0D models, including the reproduction of these differences, which is of interest, with a key feature being the probability of TA or DADs.

Download:

Fig 6. Validation of the SRF under dynamic pacing conditions.

Results of 250 simulations for 20 different pacing conditions in which notable SCRE occurred. A– 100 examples (with one highlighted) of SCRE occurring in the 3D model (purple, i-iii) and 0D model (blue, iv-v) for two different conditions which resulted in SR-Ca²⁺ close to threshold (a, corresponding to Vent EPI, BCL = 400 ms, R_CRU-CRU) and above it (b, corresponding to Vent ENDO, BCL = 400 ms, ISO + R_SERCA/NCX). The linescans in (iii) correspond to the highlighted trace in (i-ii). Final paced beat and subsequent quiescent period is shown. B–Histograms of SCRE initiation time (left of each panel) and incidence of DADs and TA (bars, right of each panel) for the 20 different conditions (panel titles correspond to cell model, pre-pacing BCL and pro-SCRE conditions); the x-axis label for the histogram plots refers to the total range over which the plot is shown, rather than absolute values. Col. RA refers to the simplified Colman et al. 2013 [27] human atrial model; ORd refers to the simplified O’Hara et al. [26] human ventricular cell model; all other labels refer to the cell-type used in the hybrid minimal model presented in this study.

https://doi.org/10.1371/journal.pcbi.1007260.g006

Emergence of focal excitation and translation of model differences

The potential for the models to simulate SCRE at the tissue and organ scale was illustrated through pacing tissue models under an equivalent rapid pacing–quiescent protocol to the single cells (see Methods: Simulation Protocols). Under the right conditions (i.e., significant SR-Ca²⁺ loading and thus large-scale release events) a triggered action potential emerged from a single focus and propagated throughout the tissue, observed in multiple tissue models (S2 Text (Validation and Results)). Multi-focal activations were also observed.

The different SCRE dynamics observed across the range of 3D single cells models under pro-SCRE conditions (Fig 6B) was accentuated in tissue, wherein focal excitations only emerged under conditions which resulted in significant TA (Fig 7). Note also the important impact of reduced I_K1, present in the atrial cell models and R_{SERCA_NCX} remodelling conditions, on allowing the emergence of tissue focal excitation.

Download:

Fig 7. The emergence of SCRE at the tissue scale in different models.

The BCL range over which activity corresponding to single cell DADs (purple) and TA (blue) and tissue focal excitation (orange) is shown for each of the different cell models/regions (A-F) under different conditions (x-axis labels–remodelling, and control and remodelling + ISO). The BCL range for which at least one SCRE/TA/focal-excitation occurs is indicated by the extent of the lines. Note that no significant SCRE was observed for control conditions for any model, and so this condition has not been included in the figure.

https://doi.org/10.1371/journal.pcbi.1007260.g007

Mechanisms of synchronisation of focal excitation

Evaluation of the focal activation emerging in the 2D sheet illustrates the mechanism by which these small scale cellular events result in full tissue excitation and the dual role of electrotonic coupling (Fig 8; S2 and S3 Videos): SCRE occurring in independent cells initially causes only a small depolarisation of the cell membrane, repressed by electrotonic coupling (Fig 8C, purple trace); however, this depolarisation is also electrotonically spread to surrounding cells (Fig 8B); due to this depolarised resting potential, SCRE occurring later can act to further depolarise the surrounding tissue (Fig 8C, blue trace); once local tissue is sufficiently depolarised, the probability of DADs manifesting as TA significantly increases and appropriately-timed firing cells can much more easily initiate a focal excitation (Fig 8C, orange trace).

Download:

Fig 8. The role of electrotonic coupling in overcoming source-sink mismatch.

A–Activation maps for intracellular Ca²⁺ (i) and membrane potential (ii) associated with a spontaneous focal excitation. Note that the time of the initiation of focal excitation (t = 0 ms) corresponds to the halfway point of the colour map. B–Temporal snapshots of Ca²⁺ (i) and V_m (ii) associated with the onset of focal excitation. C–V_m (i) and Ca²⁺ (ii) traces from individual cells from two regions with the tissue (labelled in Aii) illustrating the independent (Ca²⁺) and coupled (V_m) cellular behaviour.

https://doi.org/10.1371/journal.pcbi.1007260.g008

Inter-cellular variability and the SR-Ca²⁺ dependence of ectopic activity

In homogeneous tissue, using a baseline General Dynamic SRF model parameter set (see Methods: Simulation Protocols), the probability curve for ectopic activity was substantially steeper than that for the emergence of TA from DADs in single cell (Fig 9Ai). Reduction of the density of I_K1 shifted both the single cell and tissue TA probability curves to the left (lower SR-Ca²⁺), towards the curve describing DADs (Fig 9Aii); however, the steep relationship in tissue remains.

Download:

Fig 9. The effect of SCRE heterogeneity on the SR-Ca²⁺-TA relationship.

A–Dependence of single cell DADs (purple), single cell TA (blue), and ectopic focal activity in 2D tissue (green; square markers) on the SR-Ca²⁺ concentration in control (i) and reduced I_K1 conditions (ii); homogeneous SCRE dynamics. B–Dependence of ectopic focal activity in tissue on SR-Ca²⁺ in the heterogeneous SRF conditions with small and large variability (blue and orange) and with small and large variability with higher SR-Ca²⁺ threshold cells reassigned to baseline (blue and red; triangular markers) compared to the homogeneous condition (green); single cell DADs in the baseline condition are shown for reference (purple).

https://doi.org/10.1371/journal.pcbi.1007260.g009

With control I_K1 density, both heterogeneity models (involving either 10% or 20% of cells randomly allocated to each of the additional four General Dynamic SRF models with different SR-Ca²⁺ thresholds) shifted the SR-Ca²⁺ relationship to the right (higher SR-Ca²⁺ concentrations), despite the presence of cells more susceptible to SCRE. With reduced I_K1 density, the small heterogeneity model still shifted the relationship to the right, but the large heterogeneity model now shifted it to the left (Fig 9B). In general, the introduction of heterogeneity reduced the steepness of the probability curve at low probabilities (Fig 9B). Heterogeneity models with the less vulnerable cells removed (higher SR-Ca²⁺ threshold cells reassigned to baseline) produced the same overall features: three cases still shifted to the right despite the presence of more pro-SCRE cells and absence of less vulnerable cells.

Spontaneous Ca²⁺ release as a mechanism for conduction block

SCRE leading to non-TA inducing DADs was demonstrated to produce unidirectional conduction block during applied pacing in a homogeneous 2D model of the human ventricular wall (Fig 10A; S4 and S5 Videos) under simulated sodium channelopathy conditions (see Methods: Simulation protocols), wherein the SCRE induced DADs inactivate the sodium channel and result in non-uniform propagation following a stimulus uniformly applied to one edge (Fig 10A).

Download:

Fig 10. Mechanisms of SCRE mediated conduction block.

A–Demonstration of DAD-mediated conduction block in 2D (upper panels) and 3D (lower panels). In both cases, two stimuli were applied to one side of the tissue (left edge of the 2D sheet; ENDO surface in 3D) at a coupling interval of 500 ms, with SCRE induced DADs interrupting the second applied stimulus. Spatial snapshots cover the time just before and during this second stimulus. The locations of the cells from which the AP traces are taken are indicated by the triangular markers in the 2D sheets; for the 3D case, the traces correspond to a region which did (blue) and did not (purple) exhibit conduction block. Solid white lines represent sites of conduction block. B–Demonstration of spontaneous focal excitation leading to different behaviour in electrically homogeneous or heterogeneous tissue. The purple trace corresponds to the homogeneous condition, in which the focal excitation propagates uniformly; the blue trace corresponds to the heterogeneous condition in which focal excitation propagates non-uniformly following conduction block. The triangular marker indicates the site from which the AP traces were extracted, and the region of reduced I_K1 is highlighted by the dashed-white rectangle. The stimulus is applied once to the left edge (ENDO region) of the tissue at t = 0 ms; the second excitation is spontaneously induced.

https://doi.org/10.1371/journal.pcbi.1007260.g010

The potential for SCRE mediated focal excitation to result in unidirectional conduction block due to regional APD heterogeneity was demonstrated in the 2D transmural model of the human ventricular wall. Focal excitations originating in a small temporal window could exhibit a conduction block with the still refractory M-cell region; this was not observed in the homogeneous tissue (Fig 10B; S6 Video). Note that the site of conduction block, unlike in simulations with an applied focal stimulus, is not necessarily clearly at the boundary between the EPI and M cells: in the illustrated case, the focus emerges from the far side of the boundary.

Re-entry promotes spontaneous Ca²⁺ release and focal excitation

Sustained re-entry combined with pro-arrhythmic conditions resulted in significant SR-Ca²⁺ loading, which promoted the onset of SCRE mediated focal excitations following termination. In the illustrated 2D simulations using the minimal RA AP model, re-entry sustained for ~13 s before self-terminating, loading the SR-Ca²⁺ to ~ 1.25 mM (Fig 11A and 11B; S7–S10 Videos). By implementing a General Dynamic SRF model with the SR-Ca²⁺ threshold set below, but close to, the maximal SR-Ca²⁺ observed (CaSR_threshold = 1.125 mM), a single, delayed focal excitation was observed following termination (Fig 11A and 11Bii). With a lower threshold (CaSR_threshold = 1.0 mM) multiple and rapid focal excitations were observed which perpetuated the arrhythmic dynamics for the duration of the simulation (Fig 11A and 11Biii). In general, where the tissue was sufficiently vulnerable for focal excitation to emerge, its origin was localised to the region of the scroll wave core on its final excitation pathway (Fig 11B; S8 and S9 Videos); focal excitations emerging away from the scroll wave core were also observed closer to the threshold (Fig 11C).

Download:

Fig 11. Coupling between re-entry and SCRE.

A–SR-Ca²⁺ concentration (i) and V_m (ii) associated with sustained re-entry followed by self-termination (at around 13 s), simulations without SCRE (purple) and with the General Dynamic SRF model with two different thresholds (orange, 1.125 mM; blue, 1.0 mM). B–Temporal snapshots of voltage in the 2D sheet associated with the traces shown in A, showing self-termination (i) and the emergence of delayed (ii, corresponding to the orange traces in A) and rapid (iii, corresponding to the blue trace in A) focal excitations. C–Examples of non-localised focal excitations emerging in the 2D sheet (i) and 3D whole atria models (ii). Baseline General Dynamic SRF parameters, corresponding to panel A(orange)/Bii: CaSR_threshold = 1.25 mM; CaSR_max = 1.525 mM; CaSR_{P_range} = 0.05 mM; t_i,Sep^min = 300 ms; t_i,Sep^max = 870 ms; t_i,width^min = 200 ms; t_i,width^max = 1000 ms; MD^min = 150ms; MD^max = 600 ms; λ_width^min = 70 ms; λ_width^max = 300 ms; H_width = 2.5. Parameter differences for panel A(blue)/Biii: CaSR_threshold = 1.00 mM; CaSR_max = 1.2 mM; t_i,Sep^min = 30 ms; MD^min = 50ms; λ_width^min = 20 ms. For panel Ci: MD^min = 160ms; λ_width^min = 75 ms; λ_width^max = 300 ms. For panel Cii: CaSR_threshold = 0.9 mM; CaSR_max = 1.3 mM; t_i,Sep^min = 30 ms; t_i,width^min = 20 ms; MD^min = 80ms; MD^max = 800 ms; λ_width^min = 20 ms.

https://doi.org/10.1371/journal.pcbi.1007260.g011

Localisation of focal excitation to the scroll wave core is determined by the dynamics of re-entrant excitation: the scroll wave core remains unexcited associated with each re-entrant cycle (Fig 12A), which leads to an island of high SR-Ca²⁺ in this region (Fig 12Ac). Due to this large SR-Ca²⁺ and longest recovery time, this region thus presents the earliest SCRE which may manifest as TA. This relationship is clearly illustrated by comparing the spatial distribution of recovery time (time since last AP which induced CICR) at the moment before focal excitation occurs with the associated focal activation map: the site of activation can correspond exactly (Fig 12Bi, iv) or approximately (Fig 12Bii-iii, v) to the location of longest recovery time, in both single- (Fig 12Bi-iv) and double- scroll wave (Fig 12Bv) simulations; the earlier the focal excitation occurs relative to the final re-entrant scroll wave, the stronger the correlation between location of longest recovery time and activation source (Fig 12Bvi).

Download:

Fig 12. Functional localisation of re-entry and focal activity.

A–Temporal snapshots (i-viii) of voltage (a), intracellular Ca²⁺ (b), and SR-Ca²⁺ (c) associated with the final, self-terminating re-entrant cycle and first focal excitation. Arrows in (a) indicate the conduction of the re-entrant and then focal excitations. Highlighted region (circle) illustrates the island of large SR-Ca²⁺ associated with the unexcited scroll wave core (i-iv) and its correlation with the focus of ectopic activation (vi-vii). B–Examples of recovery time maps (a) and focal activation maps (b) for 5 independent simulations (i-v; selected from simulations covering the full range of coupled AP models and tissue parameters which led to transient re-entry which sustained sufficiently to load the SR-Ca²⁺) associated with the self-termination of re-entry followed by ectopic excitation. The contour surrounding the region of longest recovery time (corresponding to the unexcited core illustrated in A) is highlighted in red in the recovery time map and green in the activation maps. (vi)—summary of the correlation between distance between from centre of the focal source to the closest edge of the region of longest recovery and the time of the focal excitation, t(focal), relative to the latest activation of the non-focal excitation. C–Mechanism switching between re-entrant and focal excitation, showing the AP from a randomly selected cell (a) and temporal snapshots associated with the transition from re-entry to focal activity (b) and focal activity to re-entry (c). Snapshots corresponds to the temporal range illustrated by the grey and white bars with solid (re-entry to focal) and dashed (focal to re-entry) borders. Parameters which led to the mechanism switching simulation: CaSR_threshold = 1.0 mM; CaSR_max = 1.2 mM; CaSR_{P_range} = 0.05 mM; t_i,Sep^min = 30 ms; t_i,Sep^max = 870 ms; t_i,width^min = 20 ms; t_i,width^max = 200 ms; MD^min = 50ms; MD^max = 800 ms; λ_width^min = 20 ms; λ_width^max = 300 ms; H_width = 2.5.

https://doi.org/10.1371/journal.pcbi.1007260.g012

Under pro-SCRE conditions comprising tight synchronisation (small width of t_i distribution at short-coupled intervals) and short-duration, large-amplitude releases, focal activity can occur at a rate comparable to re-entry; such pacing can transiently drive the excitation, potentially eventually terminating arrhythmia as the SR-Ca²⁺ depletes. Furthermore, the wavefront could degenerate back into a sustained re-entrant excitation (Fig 12C; S11 Video). The latter is promoted by the asymmetric conduction patterns emerging from rapid focal excitation interacting with the tail of the previous scroll wave or asymmetric focal excitation. During this mechanism switching, the underlying driving mechanism (focal or re-entrant) was not clear in AP traces from randomly selected cells in the 2D tissue, and could only be determined by spatio-temporally high-resolution analysis of the excitation patterns (Fig 12C).

SR-Ca²⁺ loading as a result of re-entry led to different emergent behaviour than that induced by regular rapid pacing with matched average SR-Ca²⁺ peak concentrations (Fig 13). The impact of the islands of high SR-Ca²⁺ with longer waiting times is clear, as focal excitations in re-entry simulations could emerge almost immediately following termination whereas those following regular rapid pacing exhibited a minimum delay (approximately equal to the sum of SR refilling time and earliest t_i; Fig 13Ai). However, when SR-Ca²⁺ is close to the threshold for release, focal excitations were observed following rapid pacing but not following re-entry. This was due to the heterogeneous SR-Ca²⁺ loading (and subsequent timing of SCRE) associated with re-entry, whereas regular rapid pacing led to almost homogeneous SR-Ca²⁺ loading (offset by activation time; Fig 13Aii). Focal excitations following regular rapid pacing were always plane wave, not occurring early enough (in electrically homogeneous media) to interact with the tail of the previous excitation; those following re-entry, however, exhibited asymmetric excitation patterns when occurring sufficiently early to interact with the tail of the previous excitation (Fig 13B). Note that this is due to both the earlier excitation time and the asymmetry of the previous re-entrant excitation compared to the plane wave of regular pacing.

Download:

Fig 13. Comparison of focal excitations following re-entry vs regular pacing.

Ai–Focal excitation time, t(focal), at different SR-Ca²⁺ thresholds for release under re-entry (purple, circle markers) and matched regular pacing (blue, triangular markers) conditions. t(focal) is calculated relative to the latest activation time of the final paced or re-entrant excitation. (ii)–Temporal snapshots of spatial SR-Ca²⁺ concentration in the 2D sheets, comparing the heterogeneity at equivalent time points. Bi–Categorisation of the outcomes of focal excitation under the two conditions into either non-focal (blue), symmetric focal (orange) or asymmetric focal (red). (ii)–Examples of symmetrical vs asymmetrical conduction patterns emerging following each type of excitation. The times labelled in Bii are relative to the onset of focal excitation.

https://doi.org/10.1371/journal.pcbi.1007260.g013

Discussion

Summary

An increased incidence of spontaneous calcium release events (SCRE) is frequently observed in isolated cardiomyocytes from the diseased myocardium [11–15], and their pro-arrhythmic coupling to the membrane potential through activating inward NCX current has led to the hypothesis that these dysfunctional Ca²⁺ handling phenomena play a role in the initiation and dynamics of complex arrhythmia conduction patterns. However, investigating these multi-scale mechanisms presents a significant challenge, both for experimental and simulation approaches, and thus the precise mechanisms and potential importance of these events manifesting as tissue-scale arrhythmia have yet to be fully described.

In this study, a multi-scale computational approach was developed to simulate the dynamics of stochastic SCRE in organ-scale models of cardiac excitation. The computational framework comprises a hierarchy of models (Fig 1) encompassing the microscopic- (3D cell model), mesoscopic- (0D cell model) and macroscopic-scales (tissue models). Spontaneous Release Functions (SRF; Fig 3) were used to reproduce the morphology of SCRE in the 0D cell and tissue models (Figs 6 and 7), directly translating the single-cell modelling to the tissue-scale.

The approaches were then applied to study arrhythmia mechanisms. Firstly, the dual role of electrotonic coupling in the emergence of SCRE as a full tissue focal excitation was illustrated (Fig 8). Secondly, the role of cellular variability in SCRE dynamics on the SR-Ca2+—focal excitation relationship was investigated (Fig 9). Thirdly, two mechanisms of SCRE mediated conduction block were demonstrated (Fig 10), illustrating mechanisms by which re-entry may be initiated. Finally, the long-term interactions with re-entry were investigated (Figs 11 and 12), demonstrating that sustained re-entrant excitation can load the SR-Ca²⁺ and promote SCRE mediated focal excitation, and revealing a purely functional mechanism of localisation. This mechanism also resulted in re-entrant pre-pacing leading to earlier and more asymmetric focal excitations than observed in matched regular pacing scenarios (Fig 13).

The spontaneous release functions

The present study presents a novel approach to efficiently simulate stochastic sub-cellular Ca²⁺ release dynamics through the implementation of simple and controllable analytical wave functions—the SRF. The use of these simple functions allowed direct control over SCRE dynamics as well as models fit to the behaviour of 3D single cells, facilitating both general mechanistic analysis and translation of single-cell investigations to the organ-scale.

The phenomenological approach considered whole-cell behaviour only, and was not based on capturing underlying detail. Sub-cellular heterogeneity of factors such as t-tubule density can influence the vulnerability to and dynamics of whole-cell SCRE [36,41]. This may be particularly relevant for atrial cells (which, for large mammals, exhibit variable size and t-tubule density, and for small mammals in general exhibit a lack of t-tubules [42,43]) and disease conditions (in which t-tubule density is often observed to decrease [44–46]). This top-down approach facilitates direct incorporation of the dynamics of these heterogeneous conditions, which may be essential for accurately modelling focal excitation in the atria and disease states. In general, adaptations to the SR-Ca²⁺ threshold for SCRE, frequently observed to shift towards lower concentrations in pathophysiological conditions such as heart failure [47,48] (possibly as a result of sub-cellular structural remodelling [36]), can be easily accounted for using this approach.

Conditions for the emergence of SCRE mediated focal activity

Analysis of the impact of cellular variability in the dynamics of SCRE on the probability of the emergence of an ectopic beat highlighted the complex considerations surrounding the amplitude and timing of SCRE, and its modulation by I_K1: introducing a proportion of cells more vulnerable to SCRE in combination with reduced I_K1 led to a negative shift the SR-Ca²⁺ relationship (i.e., towards lower SR-Ca²⁺ values) as well as the emergence of lower probability events, contrasting to the almost step-function response in homogeneous tissue; the same SCRE variability without reduced I_K1 right-shifted the probability curve, despite the presence of larger amplitude SCRE compared to the homogeneous condition. This important role of reduced I_K1, as previously highlighted [19,20,25,49], is due to a combination of it resulting in (i) resting potentials closer to the I_Na activation threshold and (ii) smaller repolarising current opposing the DAD, and may be of particular relevance to arrhythmia associated with heart failure remodelling, in which I_K1 is generally reduced.

Preliminary analysis was also performed to assess the criticality of parameter distributions for the emergence of focal excitation (S2 Text (Validation and Results)–section 3, Figure S4). These simulations revealed clear constraints determining thresholds between distributions which did and did not result in focal activity, and the modulation of these constraints by I_K1; however, they also revealed the challenges of predicting emergent behaviour within the threshold region.

Pro-arrhythmic feedback between re-entry and SCRE

The mechanisms by which a suitably timed focal excitation can result in conduction block and the onset of self-perpetuating re-entrant excitation have been extensively studied previously (e.g. [27]). The present study also demonstrates a novel feedback mechanism by which re-entrant excitation promotes SCRE: Combining simulations of SCRE with those of sustained re-entrant excitation demonstrated that the resulting SR-Ca²⁺ loading can promote the emergence of focal excitation following termination of re-entry, perpetuating the arrhythmic conduction patterns.

The simulations revealed a purely functional mechanism of localised ectopic and re-entrant excitation, without the requirement for a specifically vulnerable region, determined by the additional relaxation time associated with the unexcited core of a scroll wave. This localization combined with rapid focal activity resulted in excitation patterns which were highly asymmetric and almost indistinguishable from the re-entry, due to conduction block with the tail of the previous re-entrant or focal excitation. Furthermore, these asymmetric conduction patterns could result in a complete re-entrant circuit and the re-initiation of sustained re-entry. This mechanism switching may not clearly present in mapping or ECG measurements, but may have significant implications for pharmacological intervention and provide one explanation for variable and limited success.

Generalisation of the approaches

In the first instance, the Dynamic Fit SRF models are naturally generalised within the range of modifications to the AP model which do not include any differences in the underlying Ca²⁺ handling system, i.e., remodelling and regulation of the sodium and potassium currents and I_CaL: these currents do not directly affect the probability of spontaneous Ca²⁺ sparks or their propagation as whole-cell events in the model, and so the influence of their modulation on SCRE is only through their effect on SR-Ca²⁺ loading and electrotonic load. This can be further generalised if small modifications were made to the Ca²⁺ handling system in relation to J_up, J_leak and to a lesser extent I_NCX, which have only a small effect on Ca²⁺ spark propagation (in the model) and do not significantly alter the SR-dependence. This was demonstrated in the present study by the ISO model (which had no effect on RyR or NCX, but did affect J_up and LTCC open channel availability): the SRF did not require an ISO-dependent parameter in order for the 0D model to reproduce the 3D model behaviour under ISO conditions (Fig 6).

Further modifications to the Ca²⁺ handling system, which do significantly alter the SR-dependence, required the SRF models to be rederived, as was demonstrated by the two remodelling models. Whereas time-consuming, such an approach could be used to reproduce different dynamics emerging from any individual or combination of regulation and/or remodelling conditions.

The process and approaches presented in this study can be further generalised to be incorporated into other cell models. Theoretically, any “standard” non-spatial cell model which contains a rigorous model of CICR and includes RyR open state dynamics can directly include the analytical SRF. As an example, the Grandi et al., 2011 human atrial cell model [50] was selected as this includes one of the most physiological descriptions of RyR dynamics in a non-spatial cell model. Incorporation of the SRF required only a few additional lines of code (pertaining to the implementation algorithm) further to the SRF themselves, with the magnitude rescaled to reflect the maximum open-state occupancy observed in that model (S3 Text (Generalisation of approaches)). This demonstrates the potential suitability for direct integration with available contemporary, non-spatial AP models, without the requirement to replace the native intracellular Ca²⁺ handling system or develop a spatial model equivalent.

The General Dynamic approach also facilitates parameterisation to experimental data, even with limited information. To demonstrate this functionality, the data provided by Workman et al., 2012 [34] were selected (S3 Text (Generalisation of approaches)), as these data pertain only to measurement of membrane potential (and not direct measurement of Ca²⁺ dynamics) and are therefore not an ideal dataset for parameterisation of SCRE. Even under these limited conditions, it was possible to broadly reproduce the emergence of DADs and TA observed experimentally (S3 Text (Generalisation of approaches)).

Comparison to other methods for tissue-scale simulation of SCRE

Work from only a few research groups has attempted to simulate stochastic SCRE at the organ-scale [19–25]. These independent studies used alternative phenomenological approaches to overcome the inherent challenges of this multi-scale simulation as presented in this manuscript.

Further to providing an independent approach which forms a complementary tool to the previous models, which is of particular importance for theoretical investigation of systems with many unknowns and highly non-linear behaviour, the present approach differs from those previously namely in: (i) the motivation to reproduce SCRE in tissue models in-line with that observed in specific (and variable) 3D cell models, for direct translation of single cell modelling studies to the tissue-scale; (ii) the ability to directly control waveform parameters and relate observed behaviour to these parameters; (iii) an approach which readily allows direct incorporation of both limited and detailed experimental data; and (iv) the presentation of an open-source computational framework for congruent investigation of SCRE at single cell- and tissue-scales (S1 Code).

Comparison to previous mechanistic studies

Only a few computational studies have attempted to dissect the multi-scale mechanisms involved in SCRE mediated arrhythmia. Initially, the minimum tissue substrate for the emergence of focal excitations resulting from non-stochastic EADs and DADs was investigated [18], followed by demonstration of independent cellular events emerging as a focal excitation [19–25] and SCRE as a mechanism for both triggered activity and conduction block [25]. Potential interaction with extracellular matrix remodelling was demonstrated [24], as well as the potential non-linear considerations for pharmacological action on both triggers and substrate [22].

Important features observed in these previous studies are in agreement with that of the present study, i.e. in relation to the mechanism of synchronisation overcoming electrotonic load, the steep SR-Ca²⁺ relationship of ectopic activity in homogeneous tissue, the importance of I_K1 in governing the vulnerability to ectopic activity, and the potential role of non-TA inducing DADs to cause conduction abnormalities [19,20,25,49]; this study therefore provides independent validation of these features. The present study also provides novel analyses and mechanistic insight, pertaining to: (i) the analysis of SCRE vulnerability variability on the SR-Ca²⁺-TA relationship, and the quantification of parameter distribution thresholds to predict focal excitation; (ii) investigation of the potentially pro-arrhythmic bi-directional coupling between SCRE and re-entry; and (iii) demonstration of a mechanism of localisation of these phenomena which can also lead to focal-re-entrant mechanism switching.

Limitations

Due to the substantial components of this paper, full limitations associated with the models and approaches are discussed in detail in the S4 Text (Limitations) and the reader is referred there for full description of the applicability of the models in present form; here, key model limitations and those regarding the novel analysis are discussed.

The method to derive the SRF in order to reproduce behaviour of the single cell model (i.e., not the general implementations) required a large volume of computationally intensive simulations to be performed (~ 5 000 hours of computation time per condition), although these do not need to be repeated once the parameters have been derived. Alternative approaches were also considered which have potential advantages. For example, deriving a similar iterative-map approach to that presented in [51], parameterized to the spatial cell model dynamics, would require significantly less intensive simulations and perhaps provide a more robust underlying dynamic system. A purely mathematical derivation would completely circumvent this computationally intensive requirement, and was considered in the early stages of this research. Such a derivation remains an attractive prospect for a truly rigorous and portable efficient model of SCRE. However, the analytical waveform approach presented also has advantages: based on whole-cell behaviour, parameter sets could be derived to describe underlying spatial models in limitless different conditions including sub-cellular heterogeneity and its variability, which may be significantly more challenging to reproduce with an underlying dynamical system. The General Dynamic implementation also negates the requirement for simulations on which to derive the model, and permits both comprehensive mechanistic analysis and direct parameterization to experimental data.

Simulations of SCRE emerging as tissue-scale arrhythmia in general required model parameters and conditions towards the extreme of physiologically observed behaviour (i.e., corresponding to highly diseased myocardium), and therefore do not provide a complete picture of the role of SCRE across the spectrum of patients presenting arrhythmia. Further analysis and new methodological approaches will be required to practically simulate lower probability events and fully assess the role of SCRE in cardiac arrhythmia.

In the wider context, these simulations represent only initial detailed in silico analysis of the impact of SCRE at the organ scale: experimental validation, integration with biophysically detailed species and disease specific models, and parameterisation to experimentally measured SCRE statistics, are essential to translate the approaches and mechanistic insight of the present study to arrhythmia in patients.

Application of the framework

It is intended that these approaches will be further developed and incorporated with sophisticated biophysically detailed cell models and experimentally validated simulation of cellular SCRE in multiple cardiac conditions, in order to suggest new experiments and contribute to detailed analysis of the role of SCRE in cardiac arrhythmia. Demonstration of the generalisation potential of the approaches is hoped to encourage those interested researchers in the community to integrate the presented framework with their cell models and simulation studies; open-source C++ code of the entire framework and detailed documentation is therefore provided in S1 Code and available on the GitHub repository (https://github.com/michaelcolman/MSCSF).

Conclusions

The multi-scale cardiac modelling approaches described in this manuscript and accompanying model code present the possibility to model the impact of stochastic, sub-cellular calcium dynamics on organ scale arrhythmic excitation patterns with congruent detailed cellular and tissue simulations or parameterised to specific experimental datasets. Such approaches revealed multi-scale coupling between SCRE and re-entrant excitation and a purely functional mechanism for their localisation. The mechanistic insight gained from the application of these approaches may help to improve understanding and management of cardiac arrhythmia.

Supporting information

S1 Text. Model description.

Document containing full model equations and parameters for all of the cell and tissue models and SRF presented in this study.

https://doi.org/10.1371/journal.pcbi.1007260.s001

(PDF)

S2 Text. Validation and Results.

Document containing additional figures and text for the validation of the methods and additional results.

https://doi.org/10.1371/journal.pcbi.1007260.s002

(PDF)

S3 Text. Generalisation of approaches.

Document containing additional figures and text for the generalisation of approaches.

https://doi.org/10.1371/journal.pcbi.1007260.s003

(PDF)

S4 Text. Limitations.

Document containing detailed description of the limitations associated with this study.

https://doi.org/10.1371/journal.pcbi.1007260.s004

(PDF)

S1 Code. C/C++ code and documentation containing all implementations presented in the study.

https://doi.org/10.1371/journal.pcbi.1007260.s005

(TGZ)

S1 Video. Illustration of Ca²⁺ clamp and SCRE dynamics.

Video corresponds to snapshots shown in Fig 2. Left panels show proportion of open RyR (upper) and the intracellular (purple/left axis) and SR (blue/right axis) Ca²⁺ concentrations (lower). Right panel shows the intracellular Ca²⁺ concentration in the 3D volume of the idealised cell model. Video covers SR-Ca²⁺ below threshold, just above, and significantly above.

https://doi.org/10.1371/journal.pcbi.1007260.s006

(MP4)

S2 Video. Emergence of focal excitation in 2D tissue.

Intracellular Ca²⁺ concentration (left) and membrane potential (right) illustrating the emergence of a focal excitation. Corresponds to the condition shown in Fig 8.

https://doi.org/10.1371/journal.pcbi.1007260.s007

(MP4)

S3 Video. Emergence of focal excitation in 3D tissue.

Intracellular Ca²⁺ concentration (left) and membrane potential (right) illustrating the emergence of a focal excitation in the 3D atrial model.

https://doi.org/10.1371/journal.pcbi.1007260.s008

(MP4)

S4 Video. DAD mediated conduction block in 2D.

2D homogeneous sheet, showing the emergence of DADs following regular pacing before an S2 stimulus is applied to the upper edge, resulting in conduction block. Corresponds to Fig 10A—upper.

https://doi.org/10.1371/journal.pcbi.1007260.s009

(MP4)

S5 Video. DAD mediated conduction block in 3D.

3D homogeneous ventricular wedge, showing the final S1 excitation and an S2 excitation leading to conduction block; both stimuli applied to the ENDO wall. Corresponds to Fig 10A–lower.

https://doi.org/10.1371/journal.pcbi.1007260.s010

(MP4)

S6 Video. Focal mediated conduction block.

Simulation in homogeneous (upper) and heterogeneous (lower) 2D sheets, with focal excitation emerging at a similar time, leading to conduction block only in the heterogeneous case. Corresponds to Fig 10B.

https://doi.org/10.1371/journal.pcbi.1007260.s011

(MP4)

S7 Video. Focal excitation following re-entry pt 1.

Left panels show voltage (upper) and SR-Ca²⁺ (lower) from a randomly selected cell in 2D tissue, right panels show the spatial distribution of voltage (left), intracellular Ca²⁺ (middle) and SR-Ca²⁺ (right). Data show the no SRF case (purple lines; upper 2D panels) and with SRF included (red lines, lower 2D panels), set with the threshold being lower than but close to the SR-Ca²⁺ peak reached during re-entry. This simulation illustrates a single, not-functionally localised focal excitation, corresponding to Fig 11Ci.

https://doi.org/10.1371/journal.pcbi.1007260.s012

(MP4)

S8 Video. Focal excitation following re-entry pt 2.

Left panels show voltage (upper) and SR-Ca²⁺ (lower) from a randomly selected cell in 2D tissue, right panels show the spatial distribution of voltage (left), intracellular Ca²⁺ (middle) and SR-Ca²⁺ (right). Data show the no SRF case (purple lines; upper 2D panels), the SRF case from S7 Video (middle 2D panels), and an SRF case with the threshold further below the SR-Ca²⁺ than S7 Video (red lines, lower 2D panels). This new simulation illustrates a single, functionally localised focal excitation, corresponding to Fig 11Bii.

https://doi.org/10.1371/journal.pcbi.1007260.s013

(MP4)

S9 Video. Focal excitation following re-entry pt 3.

Left panels show voltage (upper) and SR-Ca²⁺ (lower) from a randomly selected cell in 2D tissue, right panels show the spatial distribution of voltage (left), intracellular Ca²⁺ (middle) and SR-Ca²⁺ (right). Data show the no SRF case (purple lines; upper 2D panels), the SRF cases from S7 and S8 Videos (middle 2D panels), and an SRF case with the threshold set further below the SR-Ca²⁺ than S7 and S8 Videos (red lines, lower 2D panels). This new simulation illustrates multiple, functionally localised focal excitations, corresponding to Fig 11Biii.

https://doi.org/10.1371/journal.pcbi.1007260.s014

(MP4)

S10 Video. Focal excitation following re-entry in 3D.

A simulation in the 3D human atrial model in which re-entry self terminates and focal excitations emerge. Left panel is voltage, middle panel is intracellular Ca²⁺ and right panel is SR-Ca²⁺_. Simulation corresponds to that shown in Fig 11Cii.

https://doi.org/10.1371/journal.pcbi.1007260.s015

(MP4)

S11 Video. Mechanism switching.

AP traces (upper panels) and snapshots of voltage and intracellular and SR Ca²⁺ (lower panels). A number of re-entrant cycles are shown, followed by rapid focal activity which degenerates back into re-entry. The driving mechanism is indicated by the colour of the AP traces–purple for re-entry and orange for focal.

https://doi.org/10.1371/journal.pcbi.1007260.s016

(MP4)

Acknowledgments

This work was undertaken on ARC2, part of the High Performance Computing facilities at the University of Leeds, UK, and would not have been possible otherwise. I would also like to thank Dr Al Benson, University of Leeds, for providing comments and feedback on the manuscript.

References

1. Krijthe BP, Kunst A, Benjamin EJ, Lip GYH, Franco OH, Hofman A, et al. Projections on the number of individuals with atrial fibrillation in the European Union, from 2000 to 2060. Eur Heart J. 2013 Sep 14;34(35):2746–51. pmid:23900699
- View Article
- PubMed/NCBI
- Google Scholar
2. Bhatnagar P, Wickramasinghe K, Williams J, Rayner M, Townsend N. The epidemiology of cardiovascular disease in the UK 2014. Heart. 2015 May 6;heartjnl-2015-307516.
- View Article
- Google Scholar
3. Ceornodolea AD, Bal R, Severens JL. Epidemiology and Management of Atrial Fibrillation and Stroke: Review of Data from Four European Countries [Internet]. Stroke Research and Treatment. 2017 [cited 2018 Mar 31]. Available from: https://www.hindawi.com/journals/srt/2017/8593207/
- View Article
- Google Scholar
4. Qu Z, Weiss JN. Mechanisms of Ventricular Arrhythmias: From Molecular Fluctuations to Electrical Turbulence. Annu Rev Physiol. 2015;77(1):29–55.
- View Article
- Google Scholar
5. Nattel S. New ideas about atrial fibrillation 50 years on. Nature. 2002 Jan 10;415(6868):219–26. pmid:11805846
- View Article
- PubMed/NCBI
- Google Scholar
6. Calkins H, Kuck KH, Cappato R, Brugada J, Camm AJ, Chen S-A, et al. 2012 HRS/EHRA/ECAS Expert Consensus Statement on Catheter and Surgical Ablation of Atrial Fibrillation: recommendations for patient selection, procedural techniques, patient management and follow-up, definitions, endpoints, and research trial design. Eur Eur Pacing Arrhythm Card Electrophysiol J Work Groups Card Pacing Arrhythm Card Cell Electrophysiol Eur Soc Cardiol. 2012 Apr;14(4):528–606.
- View Article
- Google Scholar
7. Ganesan AN, Shipp NJ, Brooks AG, Kuklik P, Lau DH, Lim HS, et al. Long-term outcomes of catheter ablation of atrial fibrillation: a systematic review and meta-analysis. J Am Heart Assoc. 2013 Mar 18;2(2):e004549. pmid:23537812
- View Article
- PubMed/NCBI
- Google Scholar
8. Lopez-Perez A, Sebastian R, Ferrero JM. Three-dimensional cardiac computational modelling: methods, features and applications. Biomed Eng OnLine [Internet]. 2015 Apr 17 [cited 2019 Jun 13];14. Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4424572/
- View Article
- Google Scholar
9. Vagos M, van Herck IGM, Sundnes J, Arevalo HJ, Edwards AG, Koivumäki JT. Computational Modeling of Electrophysiology and Pharmacotherapy of Atrial Fibrillation: Recent Advances and Future Challenges. Front Physiol [Internet]. 2018 [cited 2019 Jun 13];9. Available from: https://www.frontiersin.org/articles/10.3389/fphys.2018.01221/full
- View Article
- Google Scholar
10. Qu Z, Garfinkel A, Weiss JN, Nivala M. Multi-scale modeling in biology: how to bridge the gaps between scales? Prog Biophys Mol Biol. 2011 Oct;107(1):21–31. pmid:21704063
- View Article
- PubMed/NCBI
- Google Scholar
11. Eisner DA, Kashimura T, Venetucci LA, Trafford AW. From the ryanodine receptor to cardiac arrhythmias. Circ J Off J Jpn Circ Soc. 2009 Sep;73(9):1561–7.
- View Article
- Google Scholar
12. Nattel S, Dobrev D. The multidimensional role of calcium in atrial fibrillation pathophysiology: mechanistic insights and therapeutic opportunities. Eur Heart J. 2012 Aug;33(15):1870–7. pmid:22507975
- View Article
- PubMed/NCBI
- Google Scholar
13. Eisner D, Bode E, Venetucci L, Trafford A. Calcium flux balance in the heart. J Mol Cell Cardiol. 2013 May;58:110–7. pmid:23220128
- View Article
- PubMed/NCBI
- Google Scholar
14. Voigt N, Heijman J, Wang Q, Chiang DY, Li N, Karck M, et al. Cellular and molecular mechanisms of atrial arrhythmogenesis in patients with paroxysmal atrial fibrillation. Circulation. 2014 Jan 14;129(2):145–56. pmid:24249718
- View Article
- PubMed/NCBI
- Google Scholar
15. Landstrom AP, Dobrev D, Wehrens XHT. Calcium Signaling and Cardiac Arrhythmias. Circ Res. 2017 Jun 9;120(12):1969–93. pmid:28596175
- View Article
- PubMed/NCBI
- Google Scholar
16. Qu Z, Weiss JN. Mechanisms of ventricular arrhythmias: from molecular fluctuations to electrical turbulence. Annu Rev Physiol. 2015;77:29–55. pmid:25340965
- View Article
- PubMed/NCBI
- Google Scholar
17. Belevych AE, Terentyev D, Terentyeva R, Nishijima Y, Sridhar A, Hamlin RL, et al. The relationship between arrhythmogenesis and impaired contractility in heart failure: role of altered ryanodine receptor function. Cardiovasc Res. 2011 Jun 1;90(3):493–502. pmid:21273243
- View Article
- PubMed/NCBI
- Google Scholar
18. Xie Y, Sato D, Garfinkel A, Qu Z, Weiss JN. So little source, so much sink: requirements for afterdepolarizations to propagate in tissue. Biophys J. 2010 Sep 8;99(5):1408–15. pmid:20816052
- View Article
- PubMed/NCBI
- Google Scholar
19. Campos FO, Shiferaw Y, Prassl AJ, Boyle PM, Vigmond EJ, Plank G. Stochastic spontaneous calcium release events trigger premature ventricular complexes by overcoming electrotonic load. Cardiovasc Res. 2015 Jul 1;107(1):175–83. pmid:25969391
- View Article
- PubMed/NCBI
- Google Scholar
20. Campos FO, Shiferaw Y, Vigmond EJ, Plank G. Stochastic spontaneous calcium release events and sodium channelopathies promote ventricular arrhythmias. Chaos Woodbury N. 2017 Sep;27(9):093910.
- View Article
- Google Scholar
21. Colman MA, Parra-Rojas C, Perez-Alday EA. From microscopic calcium sparks to the ECG: Model reduction approaches for multi-scale cardiac simulation. In: 2015 Computing in Cardiology Conference (CinC). 2015. p. 325–8.
- View Article
- Google Scholar
22. Colman MA, Alday P, A E, Holden AV, Benson AP. Trigger vs. Substrate: Multi-Dimensional Modulation of QT-Prolongation Associated Arrhythmic Dynamics by a hERG Channel Activator. Front Physiol [Internet]. 2017 [cited 2017 Oct 5];8. Available from: https://www.frontiersin.org/articles/10.3389/fphys.2017.00757/full
- View Article
- Google Scholar
23. Walker MA, Gurev V, Rice JJ, Greenstein JL, Winslow RL. Estimating the probabilities of rare arrhythmic events in multiscale computational models of cardiac cells and tissue. PLoS Comput Biol. 2017 Nov;13(11):e1005783. pmid:29145393
- View Article
- PubMed/NCBI
- Google Scholar
24. Campos FO, Shiferaw Y, Weber dos Santos R, Plank G, Bishop MJ. Microscopic Isthmuses and Fibrosis Within the Border Zone of Infarcted Hearts Promote Calcium-Mediated Ectopy and Conduction Block. Front Phys [Internet]. 2018 [cited 2018 Sep 27];6. Available from: https://www.frontiersin.org/articles/10.3389/fphy.2018.00057/full
- View Article
- Google Scholar
25. Liu MB, de Lange E, Garfinkel A, Weiss JN, Qu Z. Delayed afterdepolarizations generate both triggers and a vulnerable substrate promoting reentry in cardiac tissue. Heart Rhythm Off J Heart Rhythm Soc. 2015 Oct;12(10):2115–24.
- View Article
- Google Scholar
26. O’Hara T, Virág L, Varró A, Rudy Y. Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLoS Comput Biol. 2011 May;7(5):e1002061. pmid:21637795
- View Article
- PubMed/NCBI
- Google Scholar
27. Colman MA, Aslanidi OV, Kharche S, Boyett MR, Garratt C, Hancox JC, et al. Pro-arrhythmogenic effects of atrial fibrillation-induced electrical remodelling: insights from the three-dimensional virtual human atria. J Physiol. 2013 Sep 1;591(17):4249–72. pmid:23732649
- View Article
- PubMed/NCBI
- Google Scholar
28. Benson AP, Halley G, Li P, Tong WC, Holden AV. Virtual cell and tissue dynamics of ectopic activation of the ventricles. Chaos Interdiscip J Nonlinear Sci. 2007 Mar 1;17(1):015105.
- View Article
- Google Scholar
29. Benson AP, Aslanidi OV, Zhang H, Holden AV. The canine virtual ventricular wall: a platform for dissecting pharmacological effects on propagation and arrhythmogenesis. Prog Biophys Mol Biol. 2008 Apr;96(1–3):187–208. pmid:17915298
- View Article
- PubMed/NCBI
- Google Scholar
30. Dössel O, Krueger MW, Weber FM, Schilling C, Schulze WHW, Seemann G. A framework for personalization of computational models of the human atria. Conf Proc Annu Int Conf IEEE Eng Med Biol Soc IEEE Eng Med Biol Soc Annu Conf. 2011;2011:4324–8.
- View Article
- Google Scholar
31. Krueger MW, Seemann G, Rhode K, Keller DUJ, Schilling C, Arujuna A, et al. Personalization of atrial anatomy and electrophysiology as a basis for clinical modeling of radio-frequency ablation of atrial fibrillation. IEEE Trans Med Imaging. 2013 Jan;32(1):73–84. pmid:22665507
- View Article
- PubMed/NCBI
- Google Scholar
32. Morgan R, Colman MA, Chubb H, Seemann G, Aslanidi OV. Slow Conduction in the Border Zones of Patchy Fibrosis Stabilizes the Drivers for Atrial Fibrillation: Insights from Multi-Scale Human Atrial Modeling. Front Physiol. 2016;7:474. pmid:27826248
- View Article
- PubMed/NCBI
- Google Scholar
33. Colman MA, Pinali C, Trafford AW, Zhang H, Kitmitto A. 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. 2017 Aug 31;13(8):e1005714. pmid:28859079
- View Article
- PubMed/NCBI
- Google Scholar
34. Workman AJ, Marshall GE, Rankin AC, Smith GL, Dempster J. Transient outward K+ current reduction prolongs action potentials and promotes afterdepolarisations: a dynamic-clamp study in human and rabbit cardiac atrial myocytes. J Physiol. 2012 Sep 1;590(17):4289–305. pmid:22733660
- View Article
- PubMed/NCBI
- Google Scholar
35. Wasserstrom JA, Shiferaw Y, Chen W, Ramakrishna S, Patel H, Kelly JE, et al. Variability in timing of spontaneous calcium release in the intact rat heart is determined by the time course of sarcoplasmic reticulum calcium load. Circ Res. 2010 Oct 29;107(9):1117–26. pmid:20829511
- View Article
- PubMed/NCBI
- Google Scholar
36. Song Z, Liu MB, Qu Z. Transverse tubular network structures in the genesis of intracellular calcium alternans and triggered activity in cardiac cells. J Mol Cell Cardiol. 2018 Jan;114:288–99. pmid:29217432
- View Article
- PubMed/NCBI
- Google Scholar
37. Courtemanche M, Ramirez RJ, Nattel S. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am J Physiol-Heart Circ Physiol. 1998 Jul 1;275(1):H301–21.
- View Article
- Google Scholar
38. Venetucci LA, Trafford AW, Eisner DA. Increasing ryanodine receptor open probability alone does not produce arrhythmogenic calcium waves: threshold sarcoplasmic reticulum calcium content is required. Circ Res. 2007 Jan 5;100(1):105–11. pmid:17110597
- View Article
- PubMed/NCBI
- Google Scholar
39. Biktashev VN, Holden AV. Reentrant waves and their elimination in a model of mammalian ventricular tissue. Chaos Interdiscip J Nonlinear Sci. 1998 Mar 1;8(1):48–56.
- View Article
- Google Scholar
40. Colman MA, Ni H, Liang B, Schmitt N, Zhang H. In silico assessment of genetic variation in KCNA5 reveals multiple mechanisms of human atrial arrhythmogenesis. PLOS Comput Biol. 2017 Jun 16;13(6):e1005587. pmid:28622331
- View Article
- PubMed/NCBI
- Google Scholar
41. Colman MA, Sarathy PP, MacQuiaide N, Workman AJ. 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). 2016. p. 221–4.
- View Article
- Google Scholar
42. Richards MA, Clarke JD, Saravanan P, Voigt N, Dobrev D, Eisner DA, et al. Transverse tubules are a common feature in large mammalian atrial myocytes including human. Am J Physiol—Heart Circ Physiol. 2011 Nov;301(5):H1996–2005. pmid:21841013
- View Article
- PubMed/NCBI
- Google Scholar
43. Arora R, Aistrup GL, Supple S, Frank C, Singh J, Tai S, et al. Regional Distribution of T-tubule Density in Left and Right Atria in Dogs. Heart Rhythm. 2017 Feb;14(2):273–81. pmid:27670628
- View Article
- PubMed/NCBI
- Google Scholar
44. Crossman DJ, Ruygrok PN, Ruygrok PR, Soeller C, Cannell MB. Changes in the organization of excitation-contraction coupling structures in failing human heart. PloS One. 2011 Mar 9;6(3):e17901. pmid:21408028
- View Article
- PubMed/NCBI
- Google Scholar
45. Kemi OJ, Hoydal MA, Macquaide N, Haram PM, Koch LG, Britton SL, et al. The effect of exercise training on transverse tubules in normal, remodeled, and reverse remodeled hearts. J Cell Physiol. 2011 Sep;226(9):2235–43. pmid:21660947
- View Article
- PubMed/NCBI
- Google Scholar
46. Poláková E, Sobie EA. Alterations in T-tubule and dyad structure in heart disease: challenges and opportunities for computational analyses. Cardiovasc Res. 2013 May 1;98(2):233–9. pmid:23396602
- View Article
- PubMed/NCBI
- Google Scholar
47. Brillantes AM, Allen P, Takahashi T, Izumo S, Marks AR. Differences in cardiac calcium release channel (ryanodine receptor) expression in myocardium from patients with end-stage heart failure caused by ischemic versus dilated cardiomyopathy. Circ Res. 1992 Jul;71(1):18–26. pmid:1318794
- View Article
- PubMed/NCBI
- Google Scholar
48. Marks AR. Cardiac intracellular calcium release channels: role in heart failure. Circ Res. 2000 Jul 7;87(1):8–11. pmid:10884365
- View Article
- PubMed/NCBI
- Google Scholar
49. Greene D, Shiferaw Y. Approximate analytical solutions for excitation and propagation in cardiac tissue. Phys Rev E. 2015 Apr 30;91(4):042719.
- View Article
- Google Scholar
50. Grandi E, Pandit SV, Voigt N, Workman AJ, Dobrev D, Jalife J, et al. Human Atrial Action Potential and Ca2+ Model: Sinus Rhythm and Chronic Atrial Fibrillation. Circ Res. 2011 Oct 14;109(9):1055–66. pmid:21921263
- View Article
- PubMed/NCBI
- Google Scholar
51. Alvarez-Lacalle E, Echebarria B, Spalding J, Shiferaw Y. Calcium alternans is due to an order-disorder phase transition in cardiac cells. Phys Rev Lett. 2015 Mar 13;114(10):108101. pmid:25815968
- View Article
- PubMed/NCBI
- Google Scholar

[ref1] 1. Krijthe BP, Kunst A, Benjamin EJ, Lip GYH, Franco OH, Hofman A, et al. Projections on the number of individuals with atrial fibrillation in the European Union, from 2000 to 2060. Eur Heart J. 2013 Sep 14;34(35):2746–51. pmid:23900699
View Article
PubMed/NCBI
Google Scholar

[2] View Article

[3] PubMed/NCBI

[4] Google Scholar

[ref2] 2. Bhatnagar P, Wickramasinghe K, Williams J, Rayner M, Townsend N. The epidemiology of cardiovascular disease in the UK 2014. Heart. 2015 May 6;heartjnl-2015-307516.
View Article
Google Scholar

[6] View Article

[7] Google Scholar

[ref3] 3. Ceornodolea AD, Bal R, Severens JL. Epidemiology and Management of Atrial Fibrillation and Stroke: Review of Data from Four European Countries [Internet]. Stroke Research and Treatment. 2017 [cited 2018 Mar 31]. Available from: https://www.hindawi.com/journals/srt/2017/8593207/
View Article
Google Scholar

[9] View Article

[10] Google Scholar

[ref4] 4. Qu Z, Weiss JN. Mechanisms of Ventricular Arrhythmias: From Molecular Fluctuations to Electrical Turbulence. Annu Rev Physiol. 2015;77(1):29–55.
View Article
Google Scholar

[12] View Article

[13] Google Scholar

[ref5] 5. Nattel S. New ideas about atrial fibrillation 50 years on. Nature. 2002 Jan 10;415(6868):219–26. pmid:11805846
View Article
PubMed/NCBI
Google Scholar

[15] View Article

[16] PubMed/NCBI

[17] Google Scholar

[ref6] 6. Calkins H, Kuck KH, Cappato R, Brugada J, Camm AJ, Chen S-A, et al. 2012 HRS/EHRA/ECAS Expert Consensus Statement on Catheter and Surgical Ablation of Atrial Fibrillation: recommendations for patient selection, procedural techniques, patient management and follow-up, definitions, endpoints, and research trial design. Eur Eur Pacing Arrhythm Card Electrophysiol J Work Groups Card Pacing Arrhythm Card Cell Electrophysiol Eur Soc Cardiol. 2012 Apr;14(4):528–606.
View Article
Google Scholar

[19] View Article

[20] Google Scholar

[ref7] 7. Ganesan AN, Shipp NJ, Brooks AG, Kuklik P, Lau DH, Lim HS, et al. Long-term outcomes of catheter ablation of atrial fibrillation: a systematic review and meta-analysis. J Am Heart Assoc. 2013 Mar 18;2(2):e004549. pmid:23537812
View Article
PubMed/NCBI
Google Scholar

[22] View Article

[23] PubMed/NCBI

[24] Google Scholar

[ref8] 8. Lopez-Perez A, Sebastian R, Ferrero JM. Three-dimensional cardiac computational modelling: methods, features and applications. Biomed Eng OnLine [Internet]. 2015 Apr 17 [cited 2019 Jun 13];14. Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4424572/
View Article
Google Scholar

[26] View Article

[27] Google Scholar

[ref9] 9. Vagos M, van Herck IGM, Sundnes J, Arevalo HJ, Edwards AG, Koivumäki JT. Computational Modeling of Electrophysiology and Pharmacotherapy of Atrial Fibrillation: Recent Advances and Future Challenges. Front Physiol [Internet]. 2018 [cited 2019 Jun 13];9. Available from: https://www.frontiersin.org/articles/10.3389/fphys.2018.01221/full
View Article
Google Scholar

[29] View Article

[30] Google Scholar

[ref10] 10. Qu Z, Garfinkel A, Weiss JN, Nivala M. Multi-scale modeling in biology: how to bridge the gaps between scales? Prog Biophys Mol Biol. 2011 Oct;107(1):21–31. pmid:21704063
View Article
PubMed/NCBI
Google Scholar

[32] View Article

[33] PubMed/NCBI

[34] Google Scholar

[ref11] 11. Eisner DA, Kashimura T, Venetucci LA, Trafford AW. From the ryanodine receptor to cardiac arrhythmias. Circ J Off J Jpn Circ Soc. 2009 Sep;73(9):1561–7.
View Article
Google Scholar

[36] View Article

[37] Google Scholar

[ref12] 12. Nattel S, Dobrev D. The multidimensional role of calcium in atrial fibrillation pathophysiology: mechanistic insights and therapeutic opportunities. Eur Heart J. 2012 Aug;33(15):1870–7. pmid:22507975
View Article
PubMed/NCBI
Google Scholar

[39] View Article

[40] PubMed/NCBI

[41] Google Scholar

[ref13] 13. Eisner D, Bode E, Venetucci L, Trafford A. Calcium flux balance in the heart. J Mol Cell Cardiol. 2013 May;58:110–7. pmid:23220128
View Article
PubMed/NCBI
Google Scholar

[43] View Article

[44] PubMed/NCBI

[45] Google Scholar

[ref14] 14. Voigt N, Heijman J, Wang Q, Chiang DY, Li N, Karck M, et al. Cellular and molecular mechanisms of atrial arrhythmogenesis in patients with paroxysmal atrial fibrillation. Circulation. 2014 Jan 14;129(2):145–56. pmid:24249718
View Article
PubMed/NCBI
Google Scholar

[47] View Article

[48] PubMed/NCBI

[49] Google Scholar

[ref15] 15. Landstrom AP, Dobrev D, Wehrens XHT. Calcium Signaling and Cardiac Arrhythmias. Circ Res. 2017 Jun 9;120(12):1969–93. pmid:28596175
View Article
PubMed/NCBI
Google Scholar

[51] View Article

[52] PubMed/NCBI

[53] Google Scholar

[ref16] 16. Qu Z, Weiss JN. Mechanisms of ventricular arrhythmias: from molecular fluctuations to electrical turbulence. Annu Rev Physiol. 2015;77:29–55. pmid:25340965
View Article
PubMed/NCBI
Google Scholar

[55] View Article

[56] PubMed/NCBI

[57] Google Scholar

[ref17] 17. Belevych AE, Terentyev D, Terentyeva R, Nishijima Y, Sridhar A, Hamlin RL, et al. The relationship between arrhythmogenesis and impaired contractility in heart failure: role of altered ryanodine receptor function. Cardiovasc Res. 2011 Jun 1;90(3):493–502. pmid:21273243
View Article
PubMed/NCBI
Google Scholar

[59] View Article

[60] PubMed/NCBI

[61] Google Scholar

[ref18] 18. Xie Y, Sato D, Garfinkel A, Qu Z, Weiss JN. So little source, so much sink: requirements for afterdepolarizations to propagate in tissue. Biophys J. 2010 Sep 8;99(5):1408–15. pmid:20816052
View Article
PubMed/NCBI
Google Scholar

[63] View Article

[64] PubMed/NCBI

[65] Google Scholar

[ref19] 19. Campos FO, Shiferaw Y, Prassl AJ, Boyle PM, Vigmond EJ, Plank G. Stochastic spontaneous calcium release events trigger premature ventricular complexes by overcoming electrotonic load. Cardiovasc Res. 2015 Jul 1;107(1):175–83. pmid:25969391
View Article
PubMed/NCBI
Google Scholar

[67] View Article

[68] PubMed/NCBI

[69] Google Scholar

[ref20] 20. Campos FO, Shiferaw Y, Vigmond EJ, Plank G. Stochastic spontaneous calcium release events and sodium channelopathies promote ventricular arrhythmias. Chaos Woodbury N. 2017 Sep;27(9):093910.
View Article
Google Scholar

[71] View Article

[72] Google Scholar

[ref21] 21. Colman MA, Parra-Rojas C, Perez-Alday EA. From microscopic calcium sparks to the ECG: Model reduction approaches for multi-scale cardiac simulation. In: 2015 Computing in Cardiology Conference (CinC). 2015. p. 325–8.
View Article
Google Scholar

[74] View Article

[75] Google Scholar

[ref22] 22. Colman MA, Alday P, A E, Holden AV, Benson AP. Trigger vs. Substrate: Multi-Dimensional Modulation of QT-Prolongation Associated Arrhythmic Dynamics by a hERG Channel Activator. Front Physiol [Internet]. 2017 [cited 2017 Oct 5];8. Available from: https://www.frontiersin.org/articles/10.3389/fphys.2017.00757/full
View Article
Google Scholar

[77] View Article

[78] Google Scholar

[ref23] 23. Walker MA, Gurev V, Rice JJ, Greenstein JL, Winslow RL. Estimating the probabilities of rare arrhythmic events in multiscale computational models of cardiac cells and tissue. PLoS Comput Biol. 2017 Nov;13(11):e1005783. pmid:29145393
View Article
PubMed/NCBI
Google Scholar

[80] View Article

[81] PubMed/NCBI

[82] Google Scholar

[ref24] 24. Campos FO, Shiferaw Y, Weber dos Santos R, Plank G, Bishop MJ. Microscopic Isthmuses and Fibrosis Within the Border Zone of Infarcted Hearts Promote Calcium-Mediated Ectopy and Conduction Block. Front Phys [Internet]. 2018 [cited 2018 Sep 27];6. Available from: https://www.frontiersin.org/articles/10.3389/fphy.2018.00057/full
View Article
Google Scholar

[84] View Article

[85] Google Scholar

[ref25] 25. Liu MB, de Lange E, Garfinkel A, Weiss JN, Qu Z. Delayed afterdepolarizations generate both triggers and a vulnerable substrate promoting reentry in cardiac tissue. Heart Rhythm Off J Heart Rhythm Soc. 2015 Oct;12(10):2115–24.
View Article
Google Scholar

[87] View Article

[88] Google Scholar

[ref26] 26. O’Hara T, Virág L, Varró A, Rudy Y. Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLoS Comput Biol. 2011 May;7(5):e1002061. pmid:21637795
View Article
PubMed/NCBI
Google Scholar

[90] View Article

[91] PubMed/NCBI

[92] Google Scholar

[ref27] 27. Colman MA, Aslanidi OV, Kharche S, Boyett MR, Garratt C, Hancox JC, et al. Pro-arrhythmogenic effects of atrial fibrillation-induced electrical remodelling: insights from the three-dimensional virtual human atria. J Physiol. 2013 Sep 1;591(17):4249–72. pmid:23732649
View Article
PubMed/NCBI
Google Scholar

[94] View Article

[95] PubMed/NCBI

[96] Google Scholar

[ref28] 28. Benson AP, Halley G, Li P, Tong WC, Holden AV. Virtual cell and tissue dynamics of ectopic activation of the ventricles. Chaos Interdiscip J Nonlinear Sci. 2007 Mar 1;17(1):015105.
View Article
Google Scholar

[98] View Article

[99] Google Scholar

[ref29] 29. Benson AP, Aslanidi OV, Zhang H, Holden AV. The canine virtual ventricular wall: a platform for dissecting pharmacological effects on propagation and arrhythmogenesis. Prog Biophys Mol Biol. 2008 Apr;96(1–3):187–208. pmid:17915298
View Article
PubMed/NCBI
Google Scholar

[101] View Article

[102] PubMed/NCBI

[103] Google Scholar

[ref30] 30. Dössel O, Krueger MW, Weber FM, Schilling C, Schulze WHW, Seemann G. A framework for personalization of computational models of the human atria. Conf Proc Annu Int Conf IEEE Eng Med Biol Soc IEEE Eng Med Biol Soc Annu Conf. 2011;2011:4324–8.
View Article
Google Scholar

[105] View Article

[106] Google Scholar

[ref31] 31. Krueger MW, Seemann G, Rhode K, Keller DUJ, Schilling C, Arujuna A, et al. Personalization of atrial anatomy and electrophysiology as a basis for clinical modeling of radio-frequency ablation of atrial fibrillation. IEEE Trans Med Imaging. 2013 Jan;32(1):73–84. pmid:22665507
View Article
PubMed/NCBI
Google Scholar

[108] View Article

[109] PubMed/NCBI

[110] Google Scholar

[ref32] 32. Morgan R, Colman MA, Chubb H, Seemann G, Aslanidi OV. Slow Conduction in the Border Zones of Patchy Fibrosis Stabilizes the Drivers for Atrial Fibrillation: Insights from Multi-Scale Human Atrial Modeling. Front Physiol. 2016;7:474. pmid:27826248
View Article
PubMed/NCBI
Google Scholar

[112] View Article

[113] PubMed/NCBI

[114] Google Scholar

[ref33] 33. Colman MA, Pinali C, Trafford AW, Zhang H, Kitmitto A. 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. 2017 Aug 31;13(8):e1005714. pmid:28859079
View Article
PubMed/NCBI
Google Scholar

[116] View Article

[117] PubMed/NCBI

[118] Google Scholar

[ref34] 34. Workman AJ, Marshall GE, Rankin AC, Smith GL, Dempster J. Transient outward K+ current reduction prolongs action potentials and promotes afterdepolarisations: a dynamic-clamp study in human and rabbit cardiac atrial myocytes. J Physiol. 2012 Sep 1;590(17):4289–305. pmid:22733660
View Article
PubMed/NCBI
Google Scholar

[120] View Article

[121] PubMed/NCBI

[122] Google Scholar

[ref35] 35. Wasserstrom JA, Shiferaw Y, Chen W, Ramakrishna S, Patel H, Kelly JE, et al. Variability in timing of spontaneous calcium release in the intact rat heart is determined by the time course of sarcoplasmic reticulum calcium load. Circ Res. 2010 Oct 29;107(9):1117–26. pmid:20829511
View Article
PubMed/NCBI
Google Scholar

[124] View Article

[125] PubMed/NCBI

[126] Google Scholar

[ref36] 36. Song Z, Liu MB, Qu Z. Transverse tubular network structures in the genesis of intracellular calcium alternans and triggered activity in cardiac cells. J Mol Cell Cardiol. 2018 Jan;114:288–99. pmid:29217432
View Article
PubMed/NCBI
Google Scholar

[128] View Article

[129] PubMed/NCBI

[130] Google Scholar

[ref37] 37. Courtemanche M, Ramirez RJ, Nattel S. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am J Physiol-Heart Circ Physiol. 1998 Jul 1;275(1):H301–21.
View Article
Google Scholar

[132] View Article

[133] Google Scholar

[ref38] 38. Venetucci LA, Trafford AW, Eisner DA. Increasing ryanodine receptor open probability alone does not produce arrhythmogenic calcium waves: threshold sarcoplasmic reticulum calcium content is required. Circ Res. 2007 Jan 5;100(1):105–11. pmid:17110597
View Article
PubMed/NCBI
Google Scholar

[135] View Article

[136] PubMed/NCBI

[137] Google Scholar

[ref39] 39. Biktashev VN, Holden AV. Reentrant waves and their elimination in a model of mammalian ventricular tissue. Chaos Interdiscip J Nonlinear Sci. 1998 Mar 1;8(1):48–56.
View Article
Google Scholar

[139] View Article

[140] Google Scholar

[ref40] 40. Colman MA, Ni H, Liang B, Schmitt N, Zhang H. In silico assessment of genetic variation in KCNA5 reveals multiple mechanisms of human atrial arrhythmogenesis. PLOS Comput Biol. 2017 Jun 16;13(6):e1005587. pmid:28622331
View Article
PubMed/NCBI
Google Scholar

[142] View Article

[143] PubMed/NCBI

[144] Google Scholar

[ref41] 41. Colman MA, Sarathy PP, MacQuiaide N, Workman AJ. 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). 2016. p. 221–4.
View Article
Google Scholar

[146] View Article

[147] Google Scholar

[ref42] 42. Richards MA, Clarke JD, Saravanan P, Voigt N, Dobrev D, Eisner DA, et al. Transverse tubules are a common feature in large mammalian atrial myocytes including human. Am J Physiol—Heart Circ Physiol. 2011 Nov;301(5):H1996–2005. pmid:21841013
View Article
PubMed/NCBI
Google Scholar

[149] View Article

[150] PubMed/NCBI

[151] Google Scholar

[ref43] 43. Arora R, Aistrup GL, Supple S, Frank C, Singh J, Tai S, et al. Regional Distribution of T-tubule Density in Left and Right Atria in Dogs. Heart Rhythm. 2017 Feb;14(2):273–81. pmid:27670628
View Article
PubMed/NCBI
Google Scholar

[153] View Article

[154] PubMed/NCBI

[155] Google Scholar

[ref44] 44. Crossman DJ, Ruygrok PN, Ruygrok PR, Soeller C, Cannell MB. Changes in the organization of excitation-contraction coupling structures in failing human heart. PloS One. 2011 Mar 9;6(3):e17901. pmid:21408028
View Article
PubMed/NCBI
Google Scholar

[157] View Article

[158] PubMed/NCBI

[159] Google Scholar

[ref45] 45. Kemi OJ, Hoydal MA, Macquaide N, Haram PM, Koch LG, Britton SL, et al. The effect of exercise training on transverse tubules in normal, remodeled, and reverse remodeled hearts. J Cell Physiol. 2011 Sep;226(9):2235–43. pmid:21660947
View Article
PubMed/NCBI
Google Scholar

[161] View Article

[162] PubMed/NCBI

[163] Google Scholar

[ref46] 46. Poláková E, Sobie EA. Alterations in T-tubule and dyad structure in heart disease: challenges and opportunities for computational analyses. Cardiovasc Res. 2013 May 1;98(2):233–9. pmid:23396602
View Article
PubMed/NCBI
Google Scholar

[165] View Article

[166] PubMed/NCBI

[167] Google Scholar

[ref47] 47. Brillantes AM, Allen P, Takahashi T, Izumo S, Marks AR. Differences in cardiac calcium release channel (ryanodine receptor) expression in myocardium from patients with end-stage heart failure caused by ischemic versus dilated cardiomyopathy. Circ Res. 1992 Jul;71(1):18–26. pmid:1318794
View Article
PubMed/NCBI
Google Scholar

[169] View Article

[170] PubMed/NCBI

[171] Google Scholar

[ref48] 48. Marks AR. Cardiac intracellular calcium release channels: role in heart failure. Circ Res. 2000 Jul 7;87(1):8–11. pmid:10884365
View Article
PubMed/NCBI
Google Scholar

[173] View Article

[174] PubMed/NCBI

[175] Google Scholar

[ref49] 49. Greene D, Shiferaw Y. Approximate analytical solutions for excitation and propagation in cardiac tissue. Phys Rev E. 2015 Apr 30;91(4):042719.
View Article
Google Scholar

[177] View Article

[178] Google Scholar

[ref50] 50. Grandi E, Pandit SV, Voigt N, Workman AJ, Dobrev D, Jalife J, et al. Human Atrial Action Potential and Ca2+ Model: Sinus Rhythm and Chronic Atrial Fibrillation. Circ Res. 2011 Oct 14;109(9):1055–66. pmid:21921263
View Article
PubMed/NCBI
Google Scholar

[180] View Article

[181] PubMed/NCBI

[182] Google Scholar

[ref51] 51. Alvarez-Lacalle E, Echebarria B, Spalding J, Shiferaw Y. Calcium alternans is due to an order-disorder phase transition in cardiac cells. Phys Rev Lett. 2015 Mar 13;114(10):108101. pmid:25815968
View Article
PubMed/NCBI
Google Scholar

[184] View Article

[185] PubMed/NCBI

[186] Google Scholar

Figures

Abstract

Author summary

Introduction

Methods

The computational framework

The microscopic, 3D spatio-temporal Ca2+ handling model.

The deterministic, 0D model.

Action potential and tissue models.

Pro-arrhythmic conditions

Partial Ca2+-clamp protocol

Derivation of the Spontaneous Release Functions (SRF)

The Spontaneous Release Functions (SRF).

Parameter distributions and inverse functions.

Sampling the parameters: The direct control SRF model.

Sampling the parameters: The dynamic fit SRF model.

Sampling the parameters: The general dynamic SRF model.

Implementation with the 0D cell model.

Simulation protocols

Rapid pacing protocol.

Assessing the relationship between SR-Ca2+ and focal excitation.

Studying unidirectional conduction block.

Re-entry and SCRE interactions.

Results

Validation of the dynamic fit spontaneous release functions

Emergence of focal excitation and translation of model differences

Mechanisms of synchronisation of focal excitation

Inter-cellular variability and the SR-Ca2+ dependence of ectopic activity

Spontaneous Ca2+ release as a mechanism for conduction block

Re-entry promotes spontaneous Ca2+ release and focal excitation

Discussion

Summary

The spontaneous release functions

Conditions for the emergence of SCRE mediated focal activity

Pro-arrhythmic feedback between re-entry and SCRE

Generalisation of the approaches

Comparison to other methods for tissue-scale simulation of SCRE

Comparison to previous mechanistic studies

Limitations

Application of the framework

Conclusions

Supporting information

S1 Text. Model description.

S2 Text. Validation and Results.

S3 Text. Generalisation of approaches.

S4 Text. Limitations.

S1 Code. C/C++ code and documentation containing all implementations presented in the study.

S1 Video. Illustration of Ca2+ clamp and SCRE dynamics.

S2 Video. Emergence of focal excitation in 2D tissue.

S3 Video. Emergence of focal excitation in 3D tissue.

S4 Video. DAD mediated conduction block in 2D.

S5 Video. DAD mediated conduction block in 3D.

S6 Video. Focal mediated conduction block.

S7 Video. Focal excitation following re-entry pt 1.

S8 Video. Focal excitation following re-entry pt 2.

S9 Video. Focal excitation following re-entry pt 3.

S10 Video. Focal excitation following re-entry in 3D.

S11 Video. Mechanism switching.

Acknowledgments

References

The microscopic, 3D spatio-temporal Ca²⁺ handling model.

Partial Ca²⁺-clamp protocol

Assessing the relationship between SR-Ca²⁺ and focal excitation.

Inter-cellular variability and the SR-Ca²⁺ dependence of ectopic activity

Spontaneous Ca²⁺ release as a mechanism for conduction block

Re-entry promotes spontaneous Ca²⁺ release and focal excitation

S1 Video. Illustration of Ca²⁺ clamp and SCRE dynamics.