Calcium Sparks in Cardiac Cells in Silico

We simulate elementary calcium release events (sparks) in a single calcium release unit in ventricular myocyte. Previously developed and tested electron-conformational model of the stochastic dynamics of RyR-channels is integrated to the calcium dynamics model in the cardiac cell. This approach allows to observe RyRs opening/closing in details on the macromolecular level during the calcium dynamics course. We simulate calcium diffusion in the dyadic space and “domino-like” RyR’s activation during the so-called “calcium induced-calcium release process”. Ca2+ sparks initiation, spread and termination are investigated in the computer experiments. Sparks’ initiation and termination rate dependence on the Ca2+ diffusion velocity is observed. We show that sarcoplasmic reticulum lumen local depletion and RyR’s stochastic attrition could be the reasons of Ca2+ spark


Introduction
In the cardiac cells action potential generation and contraction is a result of a complex process of triggered calcium release from sarcoplasmic reticulum.Elementary events of Ca 2+ release from the sarcoplasmic reticulum (SR) through ryanodine receptors (RyR-channels or RyRs) are called calcium sparks [1] which are initiated by Ca 2+ influx to the cell through L-type channels rising Ca 2+ concentration enough to activate RyR-channels.This phenomenon is called Calcium Induced Calcium Release (CICR) [2] and it plays a major role in the calcium signaling in cardiac and muscle cells.A detailed study of calcium sparks is essential for understanding the role of calcium signaling in cardiac cells [3][4][5].
Without the extracellular signal spontaneous calcium sparks can be a reason of the action potential generation automaticity in cardiac pacemaker cells, however in ventricular cardiac cells they may trigger calcium waves and can be a reason of extrasystolic contraction and act as a root of arrhythmias [6].
There is an urgent need of an adequate mathematical model which can be used to simulate calcium sparks and can describe RyR's activation process in details.Here we introduce computational model of spark dynamics which includes Electron-Conformation Model (ECM) [7][8][9][10] (a detailed theory of RyR-channels stochastic dynamics) and a theory of calcium and calcium buffers (calmodulin and calsequestrin) dynamics in the calcium release units in cardiac cells.ECM is a simple and biophysically reasonable theory of RyR's stochastic dynamics describes RyR behavior in terms of open and closed electronic state.In addition to the fast electronic degree of freedom, RyR channels are characterized by the slow classical conformational coordinate, which specifies the RyR channel conductance.The RyR gating implies a conformational Langevin dynamics, Ca 2+ -induced electronic transitions, quantum tunneling and thermal transitions [8,9].
The problem of the calcium diffusion process in the compartments of the cardiac cells traditionally is solved by finite elements method due to the complex geometry of the calcium release system [2,11].For the quick computations we use finite differences method to solve the 2D problem of Ca 2+ diffusion in the dyadic space of RU which is described by the standard diffusion equation.The main advantage of our approach is a robustness of the computer simulations.
In the current paper we present our computer experiments results of calcium sparks initiationspread-termination process.
The reason of Ca 2+ spark termination is still a subject of discussion.No single mechanism is able to explain the process of CICR extinguishing.There are different points of view on this problem.First: local Ca 2+ depletion.Calcium release reduces Ca 2+ lumen concentration thereby drops RyR's activation probability from the lumen side [1,2].Second: RyR's stochastic attrition.Random closing of RyRchannels is able to stop spark spreading process [2].One of the main aims of our research is to prove these ideas.

Model geometry
Calcium release system in cardiac cells (pic.1a)consists of sarcoplasmic reticulum network and SR terminal cisternae (lumens) which are located near t-tubule (small tubules which run transversely through a muscle fiber and through which electrical impulses are transmitted from the sarcoplasm to the fiber's interior).
Pic. 1. Illustration of the model geometry.a) Cross-section of a cylindrical model with a t-tubule in the center.b) A single Ca 2+ -release unit, containing SR-lumen (jSR, junctional sarcoplasmic reticulum), a cluster (9x9) RyR-channels and a dyadic space.RyR's square cluster is located in jSR membrane.T-tubule membrane contains potential activated L-type calcium channels.Released from jSR calcium diffuses across the dyadic space in xy-directions.Calcium buffers: calmodulin (CM) is located in jSR and calsequestrin is in the dyadic space.jSR is connected with SR-network and is refilled by a constant flux jrefill.

Model equations
Ca 2+ concentrations in the release unit (RU) compartments and Ca 2+ buffers dynamics are described by the following system: where CajSR is a lumen concentration, CaSS is an average dyadic space concentration, fCQ is a calsequestrin concentration, fCM is a calmodulin concentration, CQtot: total calsequestrin concentration.CMtot: total calmodulin concentration.kfCM , kfCQ -calmodulin and calsequestrin association constants, .kbCM , kbCQ -calmodulin and calsequestrin dissociation constants, jrefill -constant Ca 2+ flux to the SR lumen, jrel -total release flux from the lumen through all opened RyR-channels which equals , where Nopen is a number of open RyRs.
Calcium dynamics across the dyadic space is described in 2D case by the diffusion equation: where CaSS_mesh is a local dyadic space Ca 2+ concentration, d is a diffusion coefficient, frel+buff is a nonlinear term, describing Ca 2+ release to the dyadic space and calcium buffers influence.RyR's transition probability in the Electron-Conformational model is divided into fast electronic transition probability and fast tunneling probability.Electronic transition probability depends on local Ca 2+ concentration near each RyR and assumed as "threshold-like": where CaSS_crit is the threshold concentration of Ca 2+ in the dyadic space, which defines the beginning of the electronic activation of the channels.Tunneling probability depends on the lumen concentration Ptun= Ptun (CajSR) and the tunneling from the closed to the open state occurs for CajSR>KCa, as well as for CajSR<KCa transitions from the open to the closed state take place, where KCa is a critical concentration of calcium ions in the SR lumen at which minimums of conformational potential are balanced.For details see references [7][8][9].

Сomputational Methods
For the numerical solution of (2) a rather simple implicit Euler numerical scheme of the finite differences method was used.To solve the matrix system of equations taking into account the question of the problem scalability, the PETSc library is proposed (https://www.mcs.anl.gov/petsc/) for C++.This package provides for use a set of distributed data structures like parallel matrices and vectors, as well as a set of the most efficient parallel solvers of the linear equations systems and preconditioners.
Numerical simulation is formulated by considering integers m1, m2 such that hx = Dx/m1, hy = Dy/m2 are characteristics of space grid.Orthogonal uniform computational mesh includes m1 × m2 numerical nodes, hx and hy are steps by space in orthogonal directions.Time grid is defined by real numbers τ and T where τ is a time step interval and N = T/τ is a number of conducted time steps.Variable u n i, j is a value of unknown concentration CaSS(xi, yj, tn) in subspace.
Approximation of () leads to a following numerical scheme: where f n i, j is non-linear part which describes calcium release in model equation.
Dirichlet boundary conditions are described in a following form: Conjugated Gradients method is used to solve linear system based on these expressions.
The time course of the RyR-cluster opening is shown on the pic.2.In the beginning of the experiment a small number of RyRs transform to the opened state due to the tunneling process.Diffusion causes an increase of a local CaSS concentration of the nearest neighbors, switching on calcium dependent electronic stochastic transitions, thus "domino-like" CICR process occurs during the first period of Ca 2+ spark formation.
In our experiment Ca 2+ release at t=10 ms was sufficient to decrease CajSR below the threshold level KCa = 500 µM (pic.3).A reverse tunneling process from the opened to the closed state takes place, a probability of RyRs activation decreases at t>10 ms because of the lumen depletion and due to Ca 2+ diffusion over the borders of the dyadic space.
According to these observed effects we've made a conclusion, that in our model spark can terminate due to lumen depletion and a RyR's stochastic attrition also plays a major role.

Conclusions and Perspectives
Computer simulations within the framework of the united model of Ca 2+ dynamics in Ca 2+ RU is able to reproduce experimental data of calcium sparks observations [] in the cardiac cells.
Incorporation of the RyRs Electron-Conformational Model to the theory of Ca 2+ dynamics in RU allows us to describe the process of RyRs activation in details on the macromolecular level but in terms of only two parameters (KCa, CaSS_crit), responsible for RyRs activation by Ca 2+ ions from the dyadic and from the luminal site.As a summary, we've proved in our simulations that mechanisms of RyRs stochastic attrition and the process of the lumen depletion play the major role of sparks termination.
In the nearest future we plan to take into account conformational RyR-RyR coupling which can be described in terms of the Electron-Conformational Model and to explore the contribution of this phenomenon to the process of the sparks initiation-spread-termination.Previously ECM has shown a brilliant ability to describe thermal effects of the RyR conformational dynamics [8,9].Thermal fluctuation terms will be added to the system of equations of the United model to observe effects of the temperature impact on the Ca 2+ sparks dynamics.

Pic. 2 .Pic. 4 .
Simple Ca 2+ spark initiation-spread-termination process.Density plot illustrates a simulated time course of RyR's opening and spatial distribution of local Ca 2+ concentrations in the dyadic space.Blue circles correspond to open RyRs, white -closed.Pic. 3. Time course of calcium concentrations in the Ca 2+ -release unit compartments.Average Ca 2+ concentration in the dyadic space (CaSS) and lumen concentration (CajSR) during the spark initiation-spread-termination process are shown.Dashed line corresponds to critical value of CajSR .We observed a dependence of spark initiation and spread rate on the Ca 2+ diffusion coefficient in the dyadic space.As one can see on the pic. 4 the increase of diffusion coefficient accelerates RyRs opening in the cluster and the amplitude of Nopen(t).As a result the lumen depletion happens earlier, thus spark termination begins earlier.Time course of the number of opened RyR-channels in the cluster during the calcium spark for different values of diffusion rate d.