{'fraction_non_alphanumeric': 0.02620548765453841, 'fraction_numerical': 0.019442781163044626, 'mean_word_length': 4.59338061465721, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimation of stochastic behaviour in cardiac myocytes: I. Ca2+ movements inside the cytosol and sarcoplasmic reticulum on curvilinear domains

Background Since the discovery of Ca2+ 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 Ca2+ concentration the coupling of Ca2+ 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 Ca2+ 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 Ca2+ 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 2þ -induced Ca 2þ release (CICR) is the repeated process in which influx of extracellular Ca 2þ into the cytosol from L-type channels triggers a release of Ca 2þ from the sarcoplasmic reticulum (SR) through ryanodine receptors (RyRs) on the z-lines. 1 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. 2 The Ca 2þ release event appears as Ca 2þ sparks [3][4][5][6][7][8] in the cytosol and Ca 2þ blinks 9 in the SR. The successive activation of RyRs increases the local Ca 2þ concentration in the cytosol. This process is called Ca 2þ transient. [10][11][12] Ca 2þ sparks observed in the heart muscle were found to behave stochastically 4 and that release occurs at discrete sites closely associated with transverse (t)-tubules at Z-lines 13 with regular spacing of 1.8-1.9 lm longitudinally and more irregularly with a mean value of 0.79 lm transversely. 8 In cardiac ventricle cells, spontaneous propagating waves of Ca 2þ may occur due to Ca 2þ 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 2þ 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][17][18][19][20][21][22][23][24][25][26] or stochastic. 15,[27][28][29][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. 16 The Ca 2þ 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 2þ flux from the SR. The occurrence was dependent on a threshold for Ca 2þ 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. 17 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 2þ diffusion can be successfully reproduced. In these models, the diffusion of Ca 2þ sparks obeys Fick's Law. Later on, an anomalous subdiffusion model was proposed based on an FDE that represents the cytosolic Ca 2þ 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][19][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 2þ 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 2þ waves and the resulting phases of contraction using cylindrical coordinates. 25 More recently, stochasticity is introduced into models 27 in which the Ca 2þ 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 2þ wave in ventricular myocytes with heart failure. 15 These models were extended from Fick's Law-based Ca 2þ 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 2þ concentrations using only one equation to represent this coupled state. This equation is usually used to track the change of cytosolic Ca 2þ concentration in space and time. In order to measure or observe the change in sarcoplasmic Ca 2þ 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. 30 is based on two governing equations. One is based on anomalous subdiffusion that models the Ca 2þ sparks and waves in the cytosol and the other models the Ca 2þ 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 2þ 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 17 or cylindrical coordinates. 31 However, these can become impractical for more complicated shaped cells. Using a simplified domain to represent such complicated shapes can be useful, however, Pump Z-line Z-line Figure 1. Schematic view of the influx and efflux of Ca 2þ in a cardiac myocyte. Ca 2þ entry through L-type channels on the sarcolemma triggers Ca 2þ release from the SR. This is a repeated process of the CICR mechanism. The Ca 2þ released from the SR increases the concentration in the cytosol. Ca 2þ here travels to the sarcomeres and interacts with filaments and activates the contractile machinery. Upon relaxation, Ca 2þ disassociates from the filaments and is pumped back into the SR by the SERCA pumps. Modified from Arif et al. 2 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 2þ 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. 27 to include the Ca 2þ sparks coupled with Ca 2þ 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 2þ 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.

Mathematical model for Ca 2þ sparks and Ca 2þ blinks with buffer and dye
The set of differential equations (DEs) by Izu et al. 22,23,27 describing the cytosolic Ca 2þ 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 2þ diffusion with the sarcoplasmic Ca 2þ concentration. For simplicity, diffusion is restricted to follow Fick's Law and the Ca 2þ 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][19][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 ½ SR is the free Ca 2þ concentration in the SR; CaF ½ and CaB ½ are the concentrations of Ca 2þ bound dye and buffer; D Cx and D Cy are cytosolic Ca 2þ diffusion coefficients in the longitudinal and transverse directions; D Dx and D Dy are diffusion coefficients of the Ca 2þ bound dye in the longitudinal and transverse directions, respectively. Ca 2þ diffusion is reported to be anisotropic 34 with diffusion twice as fast along the sarcomere orientation than that in the transverse direction. 31 Standard values used satisfy experimental data that D Cy D Cx ¼ 0:5 and are given in Table 1 along with a short description. Fluxes due to Ca 2þ fluorescent indicator dye, J dye , and endogenous stationary buffers, J buffer , are k þ j and k À j are forward and reverse rate constants for dye (j ¼ F) and buffer (j ¼ B); F ½ T and B ½ T are total concentration of the indicator dye and buffers, values of which are given in Table 1. J pump , J Ca 2þ leak , J SR leak represent the pumping rate of SR Ca 2þ -ATPase, Ca 2þ leak from the cytosol and SR, respectively. The formulation of these terms are adopted from work by several others in the field [18][19][20]25,32,33 and shown to produce promising results. These terms are given as in which V pump , K pump and n pump represent the maximal rate, affinity constant and the Hill exponent for Ca 2þ uptake into SR, while k and k f are constants for Ca 2þ efflux from cytosol to extracellular medium and from SR into cytosol, respectively. The SR Ca 2þ release term, J CRU , is used by several others including Lu et al. 15 , Izu et al. 27 and Chen et al. 28,29 where the cytosolic Ca 2þ 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 Molar flux through the RyRs, r, is associated with the source strength, K S , and the three-dimensional source strength, r 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. 23 This approach is adopted here for a general two-dimensional case 35 (see Appendix 1 for the obtained r values for the two-dimensional case). The locations of the CRUs are identified through the Dirac The stochastic term, S, controls the firing of the CRU when a probability condition is satisfied where P C x; y; t ð Þ; K prob ; n prob ð Þ ¼ P max C n prob cyt =ðK n prob prob þ C n prob cyt Þ, P max is the maximum probability of Ca 2þ spark occurrence, K prob represents the Ca 2þ sensitivity factor and n prob is the Hill coefficient. At each time step, Dt, a uniformly distributed random number, u, is generated and compared to the probability function P calculated from the Ca 2þ concentration. A CRU fires if P=P max > u and thus S ¼ 1 otherwise (i.e. in the case where 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 2þ concentration (i.e. Ca 2þ ½ SR 1 ð Þ ), 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 2þ wave. These parameter values are given in Table 2.
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 36 is used as the numerical scheme. The ODE for the sarcoplasmic Ca 2þ 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.
In the curvilinear case, the diffusion equations (1) and (3) need to be transformed from the Cartesianbased system with the gridlines specified by x and y to the curvilinear coordinate system with gridlines specified by n and g prior to discretisation. Transforming equations (1) and (3) leads to @C cyt @t ¼ D Cx @n @x 2 þ D Cy @n @y 2 @ 2 C cyt @n 2 þ 2 D Cx @n @x @g @x þ D Cy @n @y @g @y  These equations are then solved by using the ADI method. 36 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 37 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 38 to evaluate the derivatives associated with the curvilinear coordinates on=ox, on=oy, og=ox and og=oy in equations (5) and (6).
The CRUs are distributed with l x ¼ 2 lm spacing in the longitudinal and l y ¼ 0:8 lm 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.
In the next section of this paper, temporal concentration profiles are provided for the cytosolic and sarcoplasmic Ca 2þ 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 2þ propagation is provided in each case. For the estimation of the velocity, the peak of the Ca 2þ wave is tracked from a point close to the centre to the end of either side of the domain. The distance travelled (in lm) 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 2þ waves propagate as sharp bands 39 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. 33 In the following results, Ca 2þ 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 2þ wave.
In all cases that are presented in this section, the initial cytosolic Ca 2þ concentration at rest is set to Ca 2þ ½ cyt 1 ð Þ ¼ 0:1 lM. The concentration in the SR compartment is 500-1000 lM 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 lM in the NBND system and 1000 lM 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 Dt ¼ 0:1. The spatial step sizes in the Cartesian case are set to Dx ¼ Dy ¼ 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).
Looking at the simulation time when the Ca 2þ 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 lm=s with amplitudes approximately ranging 17-20.04 lM in the Cartesian (Figure 7(a)) and 180.98 lm=s with amplitudes approximately ranging 19.23-23.79 lM in the curvilinear cases (Figure 7(c)).
There is a periodic exchange of Ca 2þ between the SR and the cytosol. These result in the mentioned Ca 2þ 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 lM while this drop in the curvilinear case is approximately between 575.9 lM and 585.9 lM (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 6 ¼ 0; J buffer 6 ¼ 0. Figures 8 and 9 show the Ca 2þ wave propagation in the presence of buffer and dye reactions. The mean wave velocity is observed to be approximately 167.15 lm/s in the Cartesian and 103.53 lm/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     The Ca 2þ 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 2þ . As Ca 2þ 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][4][5][6]8,9,17,[40][41][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 lM with a drop to 71.68-89.62 lM for SR Ca 2þ content in the Cartesian case ( Figure 11(a)) and that in the curvilinear case are 55.67-74.55 lM for spark amplitude and 26.98-80.11 lM for drop in SR Ca 2þ concentration (Figure 11(c)).

Discussion and conclusion
In this paper, the application of Ca 2þ waves on curvilinear grids is successfully demonstrated with comparisons against the Cartesian case. This enables the quantitative analysis on Ca 2þ 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 2þ 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 2þ 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 lm/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 lm/ 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 2þ 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 lm and transversely 0.4-0.8 lm 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. 43 A general buffer term was represented by a single ODE. However, it is known that Ca 2þ 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 2þ 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 2þ 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 2þ diffusion follows an anomalous subdiffusion. Fick's Law is known to produce Ca 2þ full-width at half maximum (FWHM) of 1 lm, half of that observed experimentally. Anomalous subdiffusion is known to address this issue and produce FWHM of closer to 2 lm. 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 2þ 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.
The numerical conversion for the source strength, K s , involves the Green's function obtained from the analytical solution of the diffusion problem for the two-dimensional case and convolving this with the Heaviside function. Some description on the derivation is given previously. 27 The deduced relations for each of the Cartesian and curvilinear cases are given below. For an anisotropic case where D Cx ¼ 0:3 lm 2 =ms and D Cy ¼ D Cz ¼ 0:15 lm 2 =ms, the following is obtained for the molar flux, r, that contains the source strength where E 1 x ð Þ ¼ Ð 1 z e Àt t dt is an exponential integral. Dx; Dy and Dz are distances between points in the x; y and z directions, respectively. In the case of the Cartesian coordinates, the molar flux is a scalar. However, in curvilinear coordinates, it is spatially dependent and becomes a vector.