Estimation of stochastic behaviour in cardiac myocytes: I. Ca 2+ movements inside the cytosol and sarcoplasmic reticulum on curvilinear domains

Abstract

Background

Since the discovery of Ca²⁺ sparks and their stochastic behaviour in cardiac myocytes, models have focused on the inclusion of stochasticity in their studies. While most models pay much attention to the stochastic modelling of cytosolic Ca²⁺ concentration the coupling of Ca²⁺ sparks and blinks in a stochastic model has not been explored fully. The cell morphology in in silico studies in the past is assumed to be Cartesian, spherical or cylindrical. The application on curvilinear grids can easily address certain restrictions posed by such grid set up and provide more realistic cell morphology. In this paper, we present a stochastic reaction–diffusion model that couples Ca²⁺ sparks and blinks in realistic shapes of cells in curvilinear domains.

Methodology: Transformation of the model was performed to the curvilinear coordinate system. The set of equations is used to produce Ca²⁺ waves initiated from sparks and blinks. A non-buffered and non-dyed version as well as a buffered and dyed version of these equations were studied in light of observing the dynamics on the two different systems. For comparison, results for both the Cartesian and curvilinear grids are provided.

Results and conclusions: A successful demonstration of the application of curvilinear grids serving as basis for future developments.

Introduction

Ca²⁺-induced Ca²⁺ release (CICR) is the repeated process in which influx of extracellular Ca²⁺ into the cytosol from L-type channels triggers a release of Ca²⁺ from the sarcoplasmic reticulum (SR) through ryanodine receptors (RyRs) on the z-lines. ¹ The CICR process is schematically illustrated in Figure 1. A more detailed review of the processes leading to contraction can be found in Arif et al. ² The Ca²⁺ release event appears as Ca²⁺ sparks^3–8 in the cytosol and Ca²⁺ blinks ⁹ in the SR. The successive activation of RyRs increases the local Ca²⁺ concentration in the cytosol. This process is called Ca²⁺ transient.^10–12 Ca²⁺ sparks observed in the heart muscle were found to behave stochastically ⁴ and that release occurs at discrete sites closely associated with transverse (t)-tubules at Z-lines ¹³ with regular spacing of 1.8–1.9  $\mu$ m longitudinally and more irregularly with a mean value of 0.79 $\mu$ m transversely. ⁸

Figure 1.

Schematic view of the influx and efflux of Ca²⁺ in a cardiac myocyte. Ca²⁺ entry through L-type channels on the sarcolemma triggers Ca²⁺ release from the SR. This is a repeated process of the CICR mechanism. The Ca²⁺ released from the SR increases the concentration in the cytosol. Ca²⁺ here travels to the sarcomeres and interacts with filaments and activates the contractile machinery. Upon relaxation, Ca²⁺ disassociates from the filaments and is pumped back into the SR by the SERCA pumps. Modified from Arif et al. ²

In cardiac ventricle cells, spontaneous propagating waves of Ca²⁺ may occur due to Ca²⁺ overload which is believed to occur under pathological conditions. This can affect the heart’s normal function and may lead to ventricular arrhythmias and heart failure.^14,15 Therefore, it has become increasingly important to understand the underlying mechanism of the CICR process and any experimental or mathematical studies would contribute towards this.

Several models have been used to simulate the Ca²⁺ sparks and waves. A number of them are partially based on the CICR process and the release from RyRs is described as either deterministic^16–26 or stochastic.^15,27–30 The models mainly consist of fractional differential equations (FDEs), partial differential equations (PDEs), and ordinary differential equations (ODEs). Although some of these models include buffers and dye and others do not, they all shape our understanding of the underlying mechanism of CICR in cardiac myocytes. Throughout the years models have been refined to account for more buffers and a more realistic representation of the calcium release units (CRUs).

Backx et al. described the CICR process with magnesium buffer reactions using a one-dimensional model that only takes into account longitudinal diffusion. ¹⁶ The Ca²⁺ release from the SR occurs at CRUs and is assumed to be a time dependent process with an exponential rise and fall of the Ca²⁺ flux from the SR. The occurrence was dependent on a threshold for Ca²⁺ concentration which leads to a deterministic description of this term in the model. A similar deterministic formulation was incorporated in the release flux by Soeller and Cannell which was then included in the total calcium current by volume integration. ¹⁷ The aim was to experimentally test the extraction of the release flux underlying calcium sparks including a more detailed dye and buffer binding structure.

Another deterministic construction of the SR calcium release was represented in two-dimensional models^21–23 so the spatial and anisotropic properties of the cytosolic Ca²⁺ diffusion can be successfully reproduced. In these models, the diffusion of Ca²⁺ sparks obeys Fick’s Law. Later on, an anomalous subdiffusion model was proposed based on an FDE that represents the cytosolic Ca²⁺ concentration^24,26 and it was shown to reproduce the spatial characteristic of sparks more realistically.

There have also been deterministic models with no buffer reactions.^18–20 These models focused on representing the CICR process in a compartmental model with no dye and buffer reactions. No virtual point sources of CRUs are described. Instead the SR Ca²⁺ release is associated with a nonlinear term based on the cytosolic and sarcoplasmic concentrations applied to all spatial points of the computational domain. In 2009, this two-pool CICR model was integrated into a force-calcium relation by Tracqui and Ohayon to successfully reproduce Ca²⁺ waves and the resulting phases of contraction using cylindrical coordinates. ²⁵

More recently, stochasticity is introduced into models ²⁷ in which the Ca²⁺ channel opening depends on a probability function. Lu et al. integrated rogue RyRs into this model to simulate the effect of these receptors on the Ca²⁺ wave in ventricular myocytes with heart failure. ¹⁵ These models were extended from Fick’s Law-based Ca²⁺ diffusion to anomalous subdiffusion.^28,29 This has greatly improved our understanding of the calcium dynamics under pathological conditions. The majority of stochastic models focus on the coupling of cytosolic and sarcoplasmic Ca²⁺ concentrations using only one equation to represent this coupled state. This equation is usually used to track the change of cytosolic Ca²⁺ concentration in space and time. In order to measure or observe the change in sarcoplasmic Ca²⁺ concentration a separate equation is required for this which is yet to be fully explored under stochastic conditions. One recent work by Li et al. ³⁰ is based on two governing equations. One is based on anomalous subdiffusion that models the Ca²⁺ sparks and waves in the cytosol and the other models the Ca²⁺ blinks in the SR. One limitation to this work is the use of a Cartesian domain in representing the cell structure which can cause difficulty in modelling Ca²⁺ waves in more complicated cell shapes.

The main motivation of this work concerns the cell domain considered in the computational studies. The majority of models to date, focus on representing the cell morphology using simplified domains applied to Cartesian, spherical ¹⁷ or cylindrical coordinates. ³¹ However, these can become impractical for more complicated shaped cells. Using a simplified domain to represent such complicated shapes can be useful, however, insufficient for the accurate estimation especially at the boundaries. An alternative method that can address this problem is the use of finite element method (FEM) which has been implemented before for the modelling of contraction in isolated cardiac myocytes.^25,32 In this work, an alternating direction implicit (ADI) method with curvilinear coordinate transformation technique is used. ADI is a finite difference-based method and is more computationally efficient than FEM and the method itself is also easier to implement. Finite difference-based methods are a good alternative to FEM, however, it has not been used in this area in combination with a curvilinear transformation technique. This is the first work to have ever used a finite difference-based method with a curvilinear coordinate transformation technique for estimating stochastic Ca²⁺ movement in a cardiac myocyte with and without buffers and dye.

The aims of the present paper are (1) extension of the stochastic formulation by Izu et al. ²⁷ to include the Ca²⁺ sparks coupled with Ca²⁺ blinks and its application on curvilinear domains (this allows one to work on real cell morphology); (2) representation of two different systems: (a) buffered and dyed system to account for Ca²⁺ movement observed in in vitro studies and (b) non-buffered and non-dyed system to observe the difference between the two systems; (3) comparison of the Cartesian and curvilinear domains. The model and its extension to curvilinear grids are described in the next section. The application to curvilinear grids and comparison results between models are presented in the results section with conclusions and potential future improvements discussed in the final section of this paper.

Methods

Mathematical model for Ca²⁺ sparks and Ca²⁺ blinks with buffer and dye

The set of differential equations (DEs) by Izu et al.^22,23,27 describing the cytosolic Ca²⁺ concentration with buffer and dye binding properties is modified to include an additional term for the cytosolic leak^19,20,25 and a slightly different representation of the SR leak term.^19,20,25 An extension to this model here includes a coupling of the cytosolic Ca²⁺ diffusion with the sarcoplasmic Ca²⁺ concentration. For simplicity, diffusion is restricted to follow Fick’s Law and the Ca²⁺ concentration in the SR is restricted to time diffusion and so is not diffusive in space giving rise to a simple ODE. This is adopted from several other authors^{18–20,25,32,33} and was shown to produce promising results. Therefore, for the current study it will suffice to use these assumptions for computational simplicity. Considering the above, the reaction–diffusion equations are expressed as

$$\frac{\partial C_{cyt}}{\partial t}=D_{Cx}\frac{{\partial }^2{C}_{cyt}}{\partial x^2}+D_{Cy}\frac{{\partial }^2{C}_{cyt}}{\partial y^2}+J_{dye}+J_{buffer}-J_{pump}-J_{C{a}^{2+}leak}+J_{CRU}-J_{SRleak}$$

(1)

$$\frac{\partial C_{SR}}{\partial t}=J_{SRleak}+J_{pump}-J_{CRU}+J_{dye}+J_{buffer}$$

(2)

$$\frac{\partial \left[CaF\right]}{\partial t}=D_{Dx}\frac{{\partial }^2\left[CaF\right]}{\partial x^2}+D_{Dy}\frac{{\partial }^2\left[CaF\right]}{\partial y^2}-J_{dye}$$

(3)

$$\frac{\partial \left[CaB\right]}{\partial t}=-J_{buffer}$$

(4)

where $C_{cyt}={\left[C{a}^{2+}\right]}_{cyt}$ is the free Ca²⁺ concentration in the cytosol; $C_{SR}={\left[C{a}^{2+}\right]}_{SR}$ is the free Ca²⁺ concentration in the SR; $\left[CaF\right]$ and $\left[CaB\right]$ are the concentrations of Ca²⁺ bound dye and buffer; $D_{Cx}$ and $D_{Cy}$ are cytosolic Ca²⁺ diffusion coefficients in the longitudinal and transverse directions; $D_{Dx}$ and $D_{Dy}$ are diffusion coefficients of the Ca²⁺ bound dye in the longitudinal and transverse directions, respectively. Ca²⁺ diffusion is reported to be anisotropic ³⁴ with diffusion twice as fast along the sarcomere orientation than that in the transverse direction. ³¹ Standard values used satisfy experimental data that $\frac{D_{Cy}}{D_{Cx}}=0.5$ and are given in Table 1 along with a short description. Fluxes due to Ca²⁺ fluorescent indicator dye, $J_{dye}$ , and endogenous stationary buffers, $J_{buffer}$ , are

$$J_{dye}=-k_F^{+}{C}_{cyt}\left({\left[F\right]}_T-\left[CaF\right]\right)+k_F^{-}\left[CaF\right]$$

$$J_{buffer}=-k_B^{+}{C}_{cyt}\left({\left[B\right]}_T-\left[CaB\right]\right)+k_B^{-}\left[CaB\right]$$

Table 1.

Standard parameter values.^{22,23,27–29}

Parameter	Value	Units	Notes
$D_{Cx},D_{Cy}$	0.3, 0.15	$\mu {\text{m}}^2/\text{ms}$	Ca2+ diffusion coefficients
$D_{Dx},D_{Dy}$	0.02, 0.01	$\mu {\text{m}}^2/\text{ms}$	Diffusion coefficients of dye
$l_x,l_y$	2, 0.8	$\mu \text{m}$	Spacing of CRUs in x and y
${\left[C{a}^{2+}\right]}_{cyt\left(\infty \right)}$	0.1	$\mu \text{M}$	Initial Ca2+ concentration
$F$	96500	${\text{C mol}}^{-1}$	Faraday’s constant
$V_{pump}$	0.208	$\mu \text{M}/\text{ms}$	Maximum SR pump rate
$K_{pump}$	0.184	$\mu \text{M}$	SR pump Michaelis constant
$K_{prob}$	15	$\mu \text{M}$	Ca2+ sensitivity parameter
$n_{pump},n_{prob}$	3.9, 1.6	–	Hill coefficients
$T_{open}$	10	ms	Duration of current flow through CRU
$P_{max}$	0.3	${\text{per CRU ms}}^{-1}$	Maximum probability of Ca2+ spark occurrence
${\sigma }_3$	$I_{\text{CRU}}/2F$	$\text{mole}/\text{ms}$	Three-dimensional Ca2+ source strength
$k_f^{+},k_B^{+}$	80, 100	${\left(\mu \text{Ms}\right)}^{-1}$	Forward rate constants for dye and buffer reactions
$k_f^{-},k_B^{-}$	90, 100	${\text{s}}^{-1}$	Reverse rate constants for dye and buffer reactions
${\left[F\right]}_T,{\left[B\right]}_T$	50, 123	$\mu \text{M}$	Total concentration (bound + unbound) of fluorescent dye and buffer

CRUs: calcium release units; SR: sarcoplasmic reticulum.

$k_j^{+}$ and $k_j^{-}$ are forward and reverse rate constants for dye ( $j=F$ ) and buffer ( $j=B$ ); ${\left[F\right]}_T$ and ${\left[B\right]}_T$ are total concentration of the indicator dye and buffers, values of which are given in Table 1. $J_{pump}$ , $J_{C{a}^{2+}leak}$ , $J_{SRleak}$ represent the pumping rate of SR Ca²⁺-ATPase, Ca²⁺ leak from the cytosol and SR, respectively. The formulation of these terms are adopted from work by several others in the field^{18–20,25,32,33} and shown to produce promising results. These terms are given as

$$J_{pump}=V_{pump}\frac{C_{cyt}^{npump}}{K_{pump}^{n_{pump}}+C_{cyt}^{n_{pump}}}$$

$$J_{C{a}^{2+}leak}=k{C}_{cyt}$$

$$J_{SRleak}=k_f{C}_{SR}$$

in which $V_{pump}$ , $K_{pump}$ and $n_{pump}$ represent the maximal rate, affinity constant and the Hill exponent for Ca²⁺ uptake into SR, while $k$ and $k_f$ are constants for Ca²⁺ efflux from cytosol to extracellular medium and from SR into cytosol, respectively.

The SR Ca²⁺ release term, $J_{CRU}$ , is used by several others including Lu et al. ¹⁵ , Izu et al. ²⁷ and Chen et al.^28,29 where the cytosolic Ca²⁺ wave propagation initiated from either the boundary of center was successfully reproduced and the statistics and properties were studied and compared against experimental properties or snapshots. This term is expressed as

$$J_{CRU}=\sigma {\displaystyle \sum }_{i,j}\delta \left(x-x_i\right)\delta \left(y-y_j\right)S\left(x_i,y_j,t;T_{open}\right),\sigma =K_s{\sigma }_3$$

Molar flux through the RyRs, $\sigma$ , is associated with the source strength, $K_S$ , and the three-dimensional source strength, ${\sigma }_3=I_{CRU}/2F$ . $K_s$ needs to be derived for a two-dimensional model (see equation (7) in Appendix 1).

The CRUs are placed at regular intervals of $l_x$ along the x-axis and $l_y$ along the y-axis. It is assumed that once a CRU opens it will remain open for 10 ms (i.e. $T_{open}=10$ ) after which it will switch to a closed state and will not reopen (a long refractory period). An adjustment to the source strength was previously done to approximate the molar flux of a point source. ²³ This approach is adopted here for a general two-dimensional case ³⁵ (see Appendix 1 for the obtained $\sigma$ values for the two-dimensional case). The locations of the CRUs are identified through the Dirac delta function, $\delta$ where there is a channel (i.e. $x-x_i=0$ implies that $\delta =1$ which denotes there is a channel otherwise $x-x_i\ne 0$ indicates no channel $\delta =0$ ). The stochastic term, $S$ , controls the firing of the CRU when a probability condition is satisfied

$$S\left(x_i,y_j,t;T_{open}\right)=\{\begin{array}{ll}1,\hfill & \frac{P}{P_{max}}>u\hfill \\ 0,\hfill & \frac{P}{P_{max}}\le u\hfill \end{array}$$

where $P\left(C\left(x,y,t\right),K_{prob},n_{prob}\right)=P_{max}{C}_{cyt}^{n_{prob}}/(K_{prob}^{n_{prob}}+C_{cyt}^{n_{prob}})$ , $P_{max}$ is the maximum probability of Ca²⁺ spark occurrence, $K_{prob}$ represents the Ca²⁺ sensitivity factor and $n_{prob}$ is the Hill coefficient. At each time step, $\text{Δ}t$ , a uniformly distributed random number, $u$ , is generated and compared to the probability function $P$ calculated from the Ca²⁺ concentration. A CRU fires if $P/P_{max}>u$ and thus $S=1$ otherwise (i.e. in the case where $\frac{P}{P_{max}}<u$ ) a CRU will remain closed, i.e. $S=0$ .

Cartesian and curvilinear domains

The computational domain for the cell is represented by two different coordinate systems, namely Cartesian and curvilinear. It is assumed that the SR spans the size of the respective computational domains used. Two different systems are solved on the Cartesian and curvilinear domains: (1) non-buffered and non-dyed (NBND) system (i.e. $J_{dye}=J_{buffer}=0$ ) and (2) buffered and dyed (BD) system. Certain parameter values such as the initial SR Ca²⁺ concentration (i.e. ${\left[C{a}^{2+}\right]}_{SR\left(\infty \right)}$ ), $k$ , $k_f$ , and $I_{CRU}$ required modification when moving from an NBND system to a BD system in order to achieve a strong propagation of Ca²⁺ wave. These parameter values are given in Table 2.

Table 2.

Parameter values for NBND and BD cases applied to the Cartesian and curvilinear grids.

	NBND		BD
Parameter	Cartesian/curvilinear	Source	Cartesian/curvilinear	Source	Units
${\left[C{a}^{2+}\right]}_{SR\left(\infty \right)}$	750		1000	30	$\mu \text{M}$
$k$	2	25	2.5	19,25	${\text{s}}^{-1}$
$k_f$	0.2	25	0.25	19,25	${\text{s}}^{-1}$
$I_{\text{CRU}}$	6		30		$\text{pA}$

Note: Some values used in the case of the BD system are within the range of the cited values.

NBND: non-buffered and non-dyed; BD: buffered and dyed.

A two-dimensional domain is considered with the Cartesian and curvilinear grids shown in Figure 2. The Cartesian problem (Figure 2, top) is solved with equations (1) to (4). Finite discretisation of the PDEs is used and an ADI method ³⁶ is used as the numerical scheme. The ODE for the sarcoplasmic Ca²⁺ concentration is solved by using a fourth order Adams-Bashforth (AB) method and the other ODEs are solved by using the traditional forward Euler method.

Figure 2.

Grids generated on the Cartesian (top) and curvilinear (bottom) coordinate system. In both cases, the system is discretised into 276 × 71 grid points. A cell size of 110 $\mu$ m length and 28 $\mu$ m width is considered. Data for the generation of curvilinear coordinates are extracted from successive snapshots of a real cell morphology provided by Tracqui and Ohayon at time $t_4$ . ²⁵ The CRDT algorithm used in the generation of the curvilinear coordinates ensures the smooth distribution of the grid points throughout the curvilinear domain.

In the curvilinear case, the diffusion equations (1) and (3) need to be transformed from the Cartesian-based system with the gridlines specified by x and y to the curvilinear coordinate system with gridlines specified by $\xi$ and η prior to discretisation. Transforming equations (1) and (3) leads to

$$\begin{array}{l}\frac{\partial C_{cyt}}{\partial t}=\left(D_{Cx}{\left(\frac{\partial \xi }{\partial x}\right)}^2+D_{Cy}{\left(\frac{\partial \xi }{\partial y}\right)}^2\right)\frac{{\partial }^2{C}_{cyt}}{\partial {\xi }^2}\\ +2\left(D_{Cx}\frac{\partial \xi }{\partial x}\frac{\partial \eta }{\partial x}+D_{Cy}\frac{\partial \xi }{\partial y}\frac{\partial \eta }{\partial y}\right)\frac{{\partial }^2{C}_{cyt}}{\partial \xi \partial \eta }\\ +\left(D_{Cx}{\left(\frac{\partial \eta }{\partial x}\right)}^2+D_{Cy}{\left(\frac{\partial \eta }{\partial y}\right)}^2\right)\frac{{\partial }^2{C}_{cyt}}{\partial {\eta }^2}+J_{dye}\\ +J_{buffer}-J_{pump}-J_{C{a}^{2+}leak}+J_{CRU}-J_{SRleak}\end{array}$$

(5)

$$\begin{array}{l}\frac{\partial \left[CaF\right]}{\partial t}=\left(D_{Dx}{\left(\frac{\partial \xi }{\partial x}\right)}^2+D_{Dy}{\left(\frac{\partial \xi }{\partial y}\right)}^2\right)\frac{{\partial }^2\left[CaF\right]}{\partial {\xi }^2}\\ +2\left(D_{Dx}\frac{\partial \xi }{\partial x}\frac{\partial \eta }{\partial x}+D_{Dy}\frac{\partial \xi }{\partial y}\frac{\partial \eta }{\partial y}\right)\frac{{\partial }^2\left[CaF\right]}{\partial \xi \partial \eta }\\ +\left(D_{Dx}{\left(\frac{\partial \eta }{\partial x}\right)}^2+D_{Dy}{\left(\frac{\partial \eta }{\partial y}\right)}^2\right)\frac{{\partial }^2\left[CaF\right]}{\partial {\eta }^2}-J_{dye}\end{array}$$

(6)

These equations are then solved by using the ADI method. ³⁶ The source strength becomes spatially dependent since the distance between grid points is not uniform. Equations (2) and (4) to (6) are then solved on the curvilinear grid (Figure 2, bottom).

Orthogonal curvilinear coordinates are generated using the software gridgen-c developed by Sakov and CSIRO Marine Research (2006–2017) (see acknowledgement). The code incorporates a CRDT algorithm by Driscoll and Vavasis for conformal mapping ³⁷ which smoothly distributes the grid points throughout the domain. This algorithm provides the x and y positions of the curvilinear grid which is then used along with a 10th order compact finite difference scheme ³⁸ to evaluate the derivatives associated with the curvilinear coordinates $\partial \xi /\partial x$ , $\partial \xi /\partial y$ , $\partial \eta /\partial x$ and $\partial \eta /\partial y$ in equations (5) and (6).

The CRUs are distributed with $l_x=2\mu \text{m}$ spacing in the longitudinal and $l_y=0.8\mu \text{m}$ spacing in the transverse directions. In the curvilinear case, the number of spacing between CRUs is determined by setting equal number of grid points between channels leading to a pointwise evenly spread channels in both directions. The distribution of these channels is given on both Cartesian and curvilinear grids in Figure 3.

Figure 3.

Distribution of the CRU channels on the Cartesian (left) and curvilinear (right) domains. The CRUs are uniformly distributed in the Cartesian case while in the curvilinear case, they are spaced out evenly with an equal number of grid points between them leading to the distribution to appear denser at the concave (top bending) and coarser at the convex (bottom bending) regions. Sensors are also included in the domains with their x and y coordinates in $\mu m$ and their sensor number below. The locations of the sensors are marked as black squares.

As boundary conditions zero-flux is assumed at the cell boundaries by imposing $\partial C_{cyt}/\partial x=\partial C_{cyt}/\partial y=0$ . Similar assumption is imposed on the equation for the Ca²⁺ bound dye, i.e. $\partial [CaF]/\partial x=\partial [CaF]/\partial y=0$ . These conditions are widely used in previous works^{16,21–23,27–29} and shown to provide good results.

In the next section of this paper, temporal concentration profiles are provided for the cytosolic and sarcoplasmic Ca²⁺ concentrations. Sensors are placed in both domains to record these concentration values. The exact location of these sensors is given in Figure 3. In addition to these results, the mean wave velocity of Ca²⁺ propagation is provided in each case. For the estimation of the velocity, the peak of the Ca²⁺ wave is tracked from a point close to the centre to the end of either side of the domain. The distance travelled (in $\mu \mathrm{m}$ ) is then divided by the time (in seconds) taken to reach the end. An average is taken from five simulations to provide a mean velocity value.

Results

In cardiac myocytes, Ca²⁺ waves propagate as sharp bands ³⁹ with velocities and frequencies that are high enough to allow for the simultaneous propagation of distinct fronts in a given individual cell. The latter type of wave has been called type 1 for convenience. ³³

In the following results, Ca²⁺ waves of type 1 is observed. Once sparks are initiated they develop in 10 ms and increase the concentration around the region. The probability of firing is proportionally dependent on the local concentration. Therefore, a rise in concentration increases the probability of SR release and eventually leading to the activation of neighbouring CRUs. If the concentration is sufficiently high more CRUs fire consecutively which creates a propagating Ca²⁺ wave.

In all cases that are presented in this section, the initial cytosolic Ca²⁺ concentration at rest is set to ${\left[C{a}^{2+}\right]}_{cyt\left(\infty \right)}=0.1\mu$ M. The concentration in the SR compartment is 500–1000 $\mu$ M and is 2–3 orders of magnitude greater than that in the cytosol. Therefore, in the simulations presented here the initial concentration in the SR is set to 750 $\mu$ M in the NBND system and 1000 $\mu$ M in the BD systems. Simulations begin with 5 × 4 CRUs force activated in the middle of the domain in order to set off initial sparks. The temporal grids are set up with time step size of $\text{Δ}t=0.1$ . The spatial step sizes in the Cartesian case are set to $\text{Δ}x=\text{Δ}y=0.4$ and that in the curvilinear case is spatially dependent so will vary throughout the domain.

NBND system on the Cartesian and curvilinear domains

In this section, results are presented for the NBND system. The set of equations (1) to (4) are solved on the Cartesian domain and equations (2) and (4) to (6) on the curvilinear domain. Buffer and dye reactions are neglected in the model, i.e. $J_{dye}=J_{buffer}=0$ , in order to distinguish the effects of such reactions on the wave propagation velocity and amplitude. Solving the NBND system on both the Cartesian and curvilinear domains with the CRU distribution given in Figure 3, we obtain the results in Figures 4 (Cartesian) and 5 (curvilinear). Images show the development of the initial sparks in 10 ms and its propagation towards and along the boundaries at the top and bottom until it reaches and propagates outwards from the boundaries at the ends (left and right).

Figure 4.

Ca²⁺ wave propagation on the Cartesian grid in the absence of buffer and dye reactions. Snapshots are recorded at times 10, 25, 50, 75, 125, 150, 175, 190 and 300 ms (left to right, top to bottom). The colour bar denotes the Ca²⁺ concentration in $\mu$ M.

Figure 5.

Ca²⁺ wave propagation on the curvilinear grid in the absence of buffer and dye reactions. Snapshots are recorded at times 10, 40, 70, 125, 175, 225, 275, 375 and 550 ms (left to right, top to bottom). The colour bar denotes the Ca²⁺ concentration in $\mu$ M.

Looking at the simulation time when the Ca²⁺ reaches the end boundaries it was observed that the wave propagates in less time compared to that in the curvilinear case. As can be seen in Figure 6, the peak of the wave exits the domain in 200 ms in the Cartesian (Figure 6(a)) and 400 ms in the curvilinear (Figure 6(b)) cases. The longitudinal mean velocity of the wave is calculated to be approximately 312.41 $\mu \text{m}/s$ with amplitudes approximately ranging 17–20.04 $\mu$ M in the Cartesian (Figure 7(a)) and 180.98 $\mu \text{m}/\text{s}$ with amplitudes approximately ranging 19.23–23.79 $\mu$ M in the curvilinear cases (Figure 7(c)).

Figure 6.

Cytosolic Ca²⁺ wave propagation without buffer and dye reactions. The spatial characteristics of the wave propagation are plotted on the (a) Cartesian domain at times 50, 100, 150 and 200 ms and (b) curvilinear domain at times 100, 200, 300 and 400 ms. Output is given from five separate simulations.

Figure 7.

Cytosolic (a, c) and sarcoplasmic (b, d) Ca²⁺ concentrations are plotted on the Cartesian (a, b) and curvilinear (c, d) domains. See Figure 3 in previous section for the location of the sensors. Output is given from five separate simulations.

There is a periodic exchange of Ca²⁺ between the SR and the cytosol. These result in the mentioned Ca²⁺ sparks and blinks, the simulation results of which are presented in Figure 7. There is a slight difference in the amplitude of the cytosolic concentration profiles when comparing the Cartesian and curvilinear problems. Similarly, there is not much variation in the sarcoplasmic concentration profiles when comparing the solution on the two different domains. In the Cartesian case, concentration in the SR drops to approximately 580.7–582.5 $\mu$ M while this drop in the curvilinear case is approximately between 575.9 $\mu$ M and 585.9 $\mu$ M (see Figure 7(b) and (d)).

BD system on the Cartesian and curvilinear domains

In this section, results are presented for the BD system. The same set of equations as in the previous section are solved but with reactions included, i.e. $J_{dye}\ne 0,J_{buffer}\ne 0$ . Figures 8 and 9 show the Ca²⁺ wave propagation in the presence of buffer and dye reactions. The mean wave velocity is observed to be approximately 167.15  $\mu$ m/s in the Cartesian and 103.53  $\mu$ m/s in the curvilinear case. Figure 10 shows the time courses of the wave at times $t=100,200,300$ and $400$ ms (a and c) and the normalised fluorescent indicator dye profile, $F/F_0$ (b and d). The normalised fluorescent indicator profile is recorded at the cell centre.

Figure 8.

Ca²⁺ wave propagation on the Cartesian grid in the presence of buffer and dye reactions. Snapshots are recorded at times 10, 40, 75, 110, 150, 200, 300, 350 and 700 ms (left to right, top to bottom). The colourbar denotes the Ca²⁺ concentration in $\mu$ M.

Figure 9.

Ca²⁺ wave propagation on the curvilinear grid in the presence of buffer and dye reactions. Snapshots are recorded at times 10, 40, 60, 125, 225, 400, 500, 660 and 900 ms (top to bottom left to right). The colourbar denotes the Ca²⁺ concentration in $\mu$ M.

Figure 10.

Ca²⁺ wave propagation in the presence of buffer and dye reactions. The spatial characteristics of the wave propagation are plotted at times 100, 200, 300 and 400 ms on the (a) Cartesian and (c) curvilinear domains. Output is given from five separate simulations.

The Ca²⁺ concentration increases rapidly initially thus a steep rise in the normalised fluorescent value. The concentration is high enough so that neighbouring sites will activate eventually leading to wave of Ca²⁺. As Ca²⁺ propagates from the region the concentration will decrease leading to a decrease in the normalised fluorescent values. There is a small delay from initially activated channels to the neighbouring channels becoming active which explains the next two phases in the profile of decreasing then increasing trend. $F/F_0$ profile peaks at 1.13 and 1.135 in the Cartesian and curvilinear cases, respectively. The shape observed around the peak and the pattern of rise and fall is also observed in experimental^{3–6,8,9,17,40–42} and other numerical studies.^17,21,24,26

A large variation is observed in the sparks and blinks (Figure 11). The amplitude of the sparks range 41.3–53.46 $\mu$ M with a drop to 71.68–89.62 $\mu$ M for SR Ca²⁺ content in the Cartesian case (Figure 11(a)) and that in the curvilinear case are 55.67–74.55 $\mu$ M for spark amplitude and 26.98–80.11 $\mu$ M for drop in SR Ca²⁺ concentration (Figure 11(c)).

Figure 11.

Time courses of the Ca²⁺ wave propagation located at sensors 1–5. Cytosolic Ca²⁺ concentrations are plotted on the (a) Cartesian and (c) curvilinear domains. Similarly, sarcoplasmic Ca²⁺ concentrations are plotted on the (b) Cartesian and (d) curvilinear domains. See Figure 3 in previous section for the location of the sensors. Output is given from five separate simulations.

Discussion and conclusion

In this paper, the application of Ca²⁺ waves on curvilinear grids is successfully demonstrated with comparisons against the Cartesian case. This enables the quantitative analysis on Ca²⁺ sparks/blinks and waves by taking into account more realistic cell morphology. The in silico study undertaken here included BD and NBND systems to observe the Ca²⁺ wave characteristics on the different systems. It has been shown in the case of the NBND system that although not much variation was observed in the amplitudes of the cytosolic and sarcoplasmic Ca²⁺ concentration profiles, the computed wave velocity in the Cartesian case was much higher (almost twice as much) than that in the curvilinear case (which is closer to the experimentally observed value of 100 $\mu$ m/s). After adding buffer and dye, the same parameters become invalid in this new system thus a search for new parameters which produces a meaningful solution is required. In both the BD and NBND cases, the wave velocity calculated on the curvilinear coordinate system provided a closer estimate to that observed experimentally (i.e. 100 $\mu$ m/s). In the BD cases, the fluorescent indicator dye profiles follow a pattern which is recorded experimentally. The variation in spark amplitudes as well as that in SR blinks can be attributed to the fact that the channel sites are dependent on the grid. The same study can be carried out with more information regarding their actual location which can then be applied to both the Cartesian and curvilinear domains.

The formation and development of Ca²⁺ sparks into waves appear differently in the Cartesian and curvilinear grids. The wave propagation region appears to be more elongated in the y-direction for the curvilinear cases. This may be due to either the distribution of CRUs being grid dependent or the partially contracted state of the cell morphology making the distribution of CRUs in the bottom half of the domain more dense and coarser in the top half compared to that in the central part of the cell. In either case, having a curvilinear representation of the problem is a major contribution towards a more realistic cell shape which can handle more sophisticated shapes than the one presented here. The underlying algorithms used in the generation of the curvilinear domain ensure the conformal mapping of grid points. Although some regions will be densely or coarsely populated with grid points the CRDT algorithm maintains a relatively smooth distribution. Therefore, the distribution of CRUs will be dense or coarse in the respective regions accordingly. Having a spatially dependent distribution rather than a grid dependent one will make the curvilinear model more realistic since it is known that CRUs are longitudinally 2 $\mu$ m and transversely 0.4–0.8 $\mu$ m apart. The spatial distribution of the channels can be obtained from images of the location of the individual CRUs from experiments. From this, one can directly use this information to set the exact locations on the grid or can derive a spatial distribution function of the channels. The curvilinear domain can easily be adapted to accommodate these points as well as include control functions to have more grid point populated areas possibly around the channels. ⁴³

A general buffer term was represented by a single ODE. However, it is known that Ca²⁺ can bind to buffers consisting of calmodulin, troponin C, binding sites on the SR and sarcolemmal (SL) membranes. Therefore, a more detailed model accounting for each of these buffers would improve this model. Additionally, Ca²⁺ bound buffers and dye in the SR have different reaction and diffusional properties than that in the cytosol. Therefore, an additional set of equations accounting for this would make this model more comparable to the real case scenario. For the sarcoplasmic Ca²⁺ concentration, an ODE has been used. However, diffusion is also present in this compartment so a PDE can be used to account for this.

Recently, it has been shown that cytosolic Ca²⁺ diffusion follows an anomalous subdiffusion. Fick’s Law is known to produce Ca²⁺ full-width at half maximum (FWHM) of 1 $\mu$ m, half of that observed experimentally. Anomalous subdiffusion is known to address this issue and produce FWHM of closer to 2 $\mu$ m. Here transformation to curvilinear coordinates is more straightforward with Fick’s Law-based diffusion. Bearing this in mind this simplification was maintained in order to demonstrate the application. The model presented in this paper can always be extended to an anomalous subdiffusion model. Although the current work presents a study of stochastic Ca²⁺ sparks and blinks on Cartesian and curvilinear grids it would be interesting to undertake a similar study under deterministic conditions where wave propagation is measured upon a trigger in subcellular release.