Computer algorithms for automated detection and analysis of local Ca2+ releases in spontaneously beating cardiac pacemaker cells

Local Ca2+ Releases (LCRs) are crucial events involved in cardiac pacemaker cell function. However, specific algorithms for automatic LCR detection and analysis have not been developed in live, spontaneously beating pacemaker cells. In the present study we measured LCRs using a high-speed 2D-camera in spontaneously contracting sinoatrial (SA) node cells isolated from rabbit and guinea pig and developed a new algorithm capable of detecting and analyzing the LCRs spatially in two-dimensions, and in time. Our algorithm tracks points along the midline of the contracting cell. It uses these points as a coordinate system for affine transform, producing a transformed image series where the cell does not contract. Action potential-induced Ca2+ transients and LCRs were thereafter isolated from recording noise by applying a series of spatial filters. The LCR birth and death events were detected by a differential (frame-to-frame) sensitivity algorithm applied to each pixel (cell location). An LCR was detected when its signal changes sufficiently quickly within a sufficiently large area. The LCR is considered to have died when its amplitude decays substantially, or when it merges into the rising whole cell Ca2+ transient. Ultimately, our algorithm provides major LCR parameters such as period, signal mass, duration, and propagation path area. As the LCRs propagate within live cells, the algorithm identifies splitting and merging behaviors, indicating the importance of locally propagating Ca2+-induced-Ca2+-release for the fate of LCRs and for generating a powerful ensemble Ca2+ signal. Thus, our new computer algorithms eliminate motion artifacts and detect 2D local spatiotemporal events from recording noise and global signals. While the algorithms were developed to detect LCRs in sinoatrial nodal cells, they have the potential to be used in other applications in biophysics and cell physiology, for example, to detect Ca2+ wavelets (abortive waves), sparks and embers in muscle cells and Ca2+ puffs and syntillas in neurons.


Introduction
Intracellular localized Ca 2+ releases from the sarcoplasmic reticulum (SR) via release channels (ryanodine receptors, RyRs) are crucially important events in cardiac cell physiology, as illustrated by the 'local control theory', describing excitation-contraction coupling in cardiac muscle [1]. Individual local Ca 2+ release events in ventricular myocytes, known as 'Ca 2+ sparks', have been extensively characterized [2]. They are activated by opening of L-type Ca 2+ channels via Ca 2+ -induced Ca 2+ release (CICR). Such sparks exhibit stereotyped spatiotemporal behavior, with a local rise followed by decay over a defined time course, having a duration of about 20 ms, and a size of about 2 µm. Traditionally, spark characteristics have been measured by confocal scan line techniques, and algorithms for their automatic analysis have been developed and characterized [3,4]. More recently, sparks have been studied in two-dimensions, necessitating the development of more advanced analysis programs, with the capability for analysis of spark behavior in both the x and y planes [5].
In cardiac pacemaker cells, local Ca 2+ releases (LCRs) differ from those in ventricular myocytes. They are spontaneous, rhythmical, arising stochastically during diastolic depolarization [6,7]. They can be also activated by T-type Ca 2+ channels [8] or L-type Ca 2+ channels [9], especially the low-voltage activated Cav1.3 isoform [10]. The LCRs interact with membrane bound Na + -Ca 2+ exchanger (NCX), contributing crucial inward current throughout diastolic depolarization (review [11]). These LCRs are highly heterogeneous, not stereotyped, appearing in the form of either short lived sparks or locally propagating waves, and ending by spontaneously fading and disappearing or by their fusion into the ensuing action potential (AP)-induced whole cell Ca 2+ transient. Despite their crucial physiological importance, their spatiotemporal structure has not been systematically studied, and no algorithm for the automatic and objective study of their complex characteristics has been devised.
There has been intensive study of LCRs over the last 15 years using confocal line scanning and, more recently, high speed camera recordings [12]. However, substantial progress has been hampered by the lack of availability of tools to rapidly and objectively analyze such complex recordings generated and thereafter yield full, detailed and reproducible characterization of the entire LCR ensemble. In addition to their spatiotemporal complexity, LCRs occur in morphologically highly heterogeneous, spontaneously contracting SA node cells (SANC). These characteristics have proved to be further major obstacles to automated analysis of LCRs. Programs previously developed for spark analysis in ventricular myocytes have proved to be inadequate in addressing the specific challenges of LCR analysis and characterization in SANC.
In this study, we have addressed these challenges by automating LCR analysis in two dimensional movies of Ca 2+ signals in SANC, through the development of two novel computerized algorithms. First, we developed a new ImageJ plugin that detects local curvature of the SANC, with subsequent affine transformation to render the cell visually stationary.
Thereafter, a second novel algorithm that was written in C++ objectively and automatically detects LCRs. Specifically, the second new algorithm enables complex spatiotemporal characterization of LCRs, including details of LCR birth, death, splitting, merger, and area path travelled within the cell, in addition to more established characteristics like size, intensity and duration. Such detailed characterization is crucial to give novel insights into the fine structure of the LCR ensemble and biophysical mechanisms of pacemaker cell function.

Ethics Statement
The present study conformed to the Guide for the Care and Use of Laboratory Animals, The adequacy of anesthesia was monitored in rabbits until reflexes to ear and tale pinch and jaw tone were lost. The guinea pigs (Charles River Laboratories, USA) weighed 500 -650 grams and were acutely anesthetized with pentobarbital sodium using IP injection at approximately 150 mg/kg doses to an absence of toe pinch reflex and eye membrane retraction.

SANC Preparations
Spontaneously beating SANC were isolated from guinea pig and rabbits using essentially the same procedure as previously described [13]. The heart was removed quickly and placed in the solution containing (in mM): 130 NaCl, 24 NaHCO 3 , 1.

2D Ca 2+ dynamics measurements in single SANC
Ca 2+ dynamics within isolated single SANC were measured by 2D imaging of fluorescence emitted by the Ca 2+ indicator fluo-4 using a high-speed Hamamatsu C9100-12 CCD camera with an 8.192 mm square sensor of 512 x 512 pixels resolution or a PCO.edge 4.2 CMOS camera with an 13.2 mm square sensor of 2048x2048 pixels resolution.
To resolve LCR dynamics we acquired images at a rate of 100 frames/second that was possible only using a part of the sensor (e.g. 160 x 512 pixels for Hamamatsu camera or 1280x1280 for PCO camera). In some recording with the PCO.edge 4.2 CMOS camera we used 2x2 binning. Both cameras were mounted on Zeiss Axiovert 100 inverted microscopes (Carl Zeiss, Inc., Germany) with a x63 oil immersion lens and a fluorescence excitation light source CoolLED pE-300-W (CoolLED Ltd. Andover, UK).
Fluo-4 fluorescence excitation (470/40 nm) and emission light collection (525/50 nm) were performed using the Zeiss filter set 38 HE. Cells were loaded with 1.5 M fluo-4AM (Sigma-Aldrich) for 10 minutes at room temperature. Fluo-4AM was subsequently washed out of the chamber, and Ca 2+ signals were measured within the ensuing 60 minutes at 35 o C ± 0.1°C.
To avoid phototoxicity, fluo-4 was excited only for short periods of time (<5 s). Data acquisition was performed using SimplePCI (Hamamatsu Corporation, Japan) or PCO camware 64 (PCO AG, Germany). In some experiments (n=10 cells, S1 Movie) we also recorded simultaneous membrane potentials by using perforated patch clamp as previously described [12].

Bathing solution and temperature control
All Ca 2+ and AP measurements were performed at 35°C+/-0.1°C (500 μl chamber

Software development
The SANC Analyser tracking contractile motion of the SANC was written in Java as a plugin for the ImageJ image analysis platform (NIH, Bethesda, MD), using the Netbeans IDE. The The website provides the following materials: 1) Our program's C++ code with all libraries (needed for program compilation in Visual Studio).
2) msi file for program installation in computers with Windows operating systems.
3) A detailed User Guide that describes software installation and program use to detect LCRs.

LCRs in 2D movies
Rhythmical stochastic LCR signals in spontaneously beating SANC emerge between APinduced whole cell Ca 2+ transients (see example in S1 Movie). We show the simultaneous patch clamp/2D-Ca 2+ signal data in order to emphasize the importance of having a fast, reproducible and accurate program for the analysis of LCRs in 2D so that their biophysical behavior can be linked to the behavior of the membrane potential and membrane ion channels.

SANC Analyser
Tracking Ca 2+ signals in a particular cell location is problematic when the cell moves, e.g. due to its contraction. Such movement leads to a so-called 'motion artifact', because each fixed pixel in a stack of images does not exactly report the same location or cell neighborhood both within the same beat, and from beat-to-beat in recordings of several consecutive beats.
To solve this issue, we developed a plug in for ImageJ, 'SANC Analyser', which artificially rendered the complex contracting cell stationary. This allowed identical geographical locations to be followed accurately on a within beat and between beats basis. An example of such analysis and image transformations is illustrated in S2 Movie (a contracting cell before analysis) and S3 Movie (stationary cell after transformation).
Our SANC Analyser is a set of algorithms that tracks the contractile motion of the SANC, and using the tracked points, affine transforms the image series so that the cell appears motionless. Tracking is based on the hypothesis that the local curvature of the cell's long axis is conserved under deformation/contraction, i.e. a bend in the cell is a fixed point which can be tracked during contraction.
A "spine" along the cell's long axis is determined by the centre of mass. With the cell approximately parallel to the x axis of the image series i(x,y,t), the centre of mass () along the y axis is: (1) where b and  are the mean and standard deviation of the background fluorescence. (x,t) are the y coordinates of the spine (Fig 1A), for which the local curvature (k) is, The partial derivatives of (x,t) are calculated by Savitsky-Golay filters [14].
As the spine is the centre of mass, its local curvature is a function of: 1) the overall curvature of the cell; 2) bumps in the cell's silhouette which are asymmetric with respect to the long axis and similar asymmetries in fluorescence intensity. The first two (which are really the same thing at different scales) should be conserved during contraction and are therefore a basis for tracking. Asymmetries in fluorescence intensity from LCRs will interfere with this if they are large enough. This can be seen in k(t,x), the local curvature "map" (Fig 1B). Horizontal bands of constant curvature are seen that bend when the cell contracts. Just prior to each contraction these bands become a little noisy due to large LCRs which precede and merge into the AP. This Each pair of points is the base of two isosceles triangles, one below and one above the spine, of user defined height. Each pair of neighboring triangles define an intervening triangle of opposite orientation. The grid is used to affine transform the series. Where the three points of a triangle form the matrix: the transform from {x, y} in a reference image (to which all images are to be transformed) to {x', y'} in an image to be transformed is: where M r is the reference triangle and M d the transforming triangle.

Notes for practical use of SANC Analyser
The automatic identification of contractions does require some degree of care and experimentation in terms of selecting the correct kind of contraction to make the original archetypal contraction when you hit the 'set contraction' button. For example, if too bright a contraction is selected with 'set contraction', then 'match contraction' will not find subsequent This phase of the analysis is very important, and strenuous efforts should be made to identify as many contractions as possible. If a substantial number of contractions cannot be automatically identified as detailed above, then the results from that particular cell should be discarded. For further detail please see our user program guide that is provided online together with the program code and plugin (see Methods).

Our algorithm was implemented within an open-source computer program
XYTEventDetector (see Methods). The user Interface for this program is shown in Figure 2 and a detailed user guide is provided online together with the program code and program installation file. Below we provide principles of LCR detection together with a short description of specific program interface parameters customizing the execution of these principles. The sequence of LCR detection is illustrated in Figure 3. Our algorithm isolates the cell cytoplasm from the background, and subsequently normalizes the cytoplasmic signals, eliminating the global signal of the AP-induced whole cell Ca 2+ transient. It then decreases the remaining recording noise using spatial filters that the user defines for given experimental conditions and recording noise, which will vary from cell-to-cell. Thereafter, any local signal that remains is considered a possible, or 'candidate', LCR. These possible LCRs are then further refined with threshold algorithms for signal amplitude and spatial size. To identify the pixels representing the cell cytoplasm, the program uses a custom filter, which we termed the 'Max Filter' (Fig 2, Fig 3A, Fig 3B). Once we have identified the pixels of cell cytoplasm F(x, y, t), we go on to normalize their signal (Fig 3C), taking 255 as the number of levels of gray. It is executed by the "Normalization" box in the program interface. This is the first step in separating LCRs from noise: If there are AP-induced Ca transients in the data, it is important to separate the true local rise of the LCR signal from the global transient. This is accomplished by the "Decay Filter" in the program interface. We take the global transient as the average signal over the entire cell, F a (t) at each time frame. Then we subtract F a (t) from each pixel F norm (x,y,t) value inside the cell area over the entire time series, and any pixel value that becomes negative is clipped to be 0 ( Fig 3D): Then to identify cell locations in which the signal changes over time, we apply a custom differential filter (Fig 3E, "Differential Filter" box in user interface). This is a composite spatiotemporal filter that integrates the signal spatially in all frames and then calculates frameto-frame differences in each pixel. The spatial averaging width is defined by interface parameter "Search Distance", that is the distance around each pixel used for averaging. By default, Search Distance is set to 3 and defines a grid area of 7x7 pixels with the tested pixel in the middle of this grid. This default setting was used in majority of our LCR detection analysis with the Hamamatsu C9100-12 CCD camera, defining the averaging area of 7*7*0.2539 2 = 3.16 µm 2 . Thus, this complex filter excludes pixel-to-pixel noise and it is also very effective in excluding locations with constant fluorescence due to heterogeneous indicator loading within the cell cytoplasm and organelles such as nuclei. This identifies locations of potential LCRrelated Ca 2+ dynamics and it is important for subsequent detection of LCR birth/death events.
The "Differential Filter" performs further important discriminations based on signal strength, spatial distribution, and its timing with respect to the beginning of the AP-induced Ca transient. This function can be customized by tuning parameters described below for specific recordings. We first renormalize the pixels using equations 6 and 7. For a pixel to be considered as a part of a newborn or existing LCR, the signals must surpass a certain threshold given by the parameter "SD Detection" in the user interface, which is set by the user as a factor of the standard deviation of the signal in the cell. By default, this parameter is set at one standard deviation. We also assume that if a pixel is inside a possible LCR then it is also part of the LCR. determines whether the pixel intensity drops low enough to be considered no a part of an LCR.
The corresponding input parameter is "SD Termination".
We have also two additional (optional) filters, a median and a mean filter, to get rid of any remaining noise, depending on the presence of 'salt-and-pepper noise' (sparsely occurring white and black pixels, also known as 'data drop-out noise', 'intensity spikes' or 'speckle').
So far by using different types of filtering, we effectively reduced noise and identified locations potentially linked to LCR activity, i.e. birth and death. The next step is to identify the new LCRs and to define the dynamics of existing LCRs (Fig 3F, Fig 3G). This procedure starts with the "Find LCRs" button in the interface. In addition to LCR detection sensitivity parameters that exist within the "Differential Filter" described above, there are 2 additional sensitivity parameters to further fine tune LCR detection for optimal performance in rejecting noisy, smaller LCRs or to tailor the program to the specific aims of a particular study (for example, to detect selectively only large LCRs).
The program searches for clusters, i.e. sets of pixels having common borders, which would satisfy specific conditions as described below. A newborn LCR, a cluster of pixels of any morphology, in a certain frame must be composed of a minimum number of pixels, as defined by the parameter "Size Threshold". In practice, this determines the minimum starting size of an LCR in pixels. Increasing the value gets rid of noise but risks failing to detect smaller LCRs. A newborn LCR must be also sufficiently bright (on a scale of arbitrary units 0-255 in our program). It is important to note that the two sensitivity tests act together with the Boolean function AND, i.e. a newborn LCR must be large and bright enough to pass these tests otherwise it is considered as noise and therefore discarded.
The program subsequently looks further to see if any of these candidates are representative of the dynamics of an existing LCR. We consider that if a candidate (or a cluster of candidates) has a common border with any existing LCR, it becomes a part of that LCRs.
Finally, a cluster becomes a newborn LCR if its size and total signal surpass respective userdefined threshold levels. These two sensitivity parameters rejecting small signal fluctuations are given in the program interface as "Size Threshold" and "Intensity Threshold", representing respectively the total number of pixels and the sum of intensities in the smallest LCR.
After a pixel is considered a part of an LCR, it will continue to be part of that LCR until the LCR is presumed terminated, for example by a collision or by merging into the Ca 2+ transient. In the case of LCR decay, a pixel is considered part of a possible LCR until the average signal in its series becomes less than a factor of the standard deviation that is also by default set at one in the interface parameter "SD Termination" described above.
With dynamics in play (Fig 4, S4 Movie), we assume that when an LCR separates, the separated groups of pixels are still part of the same LCR. When LCRs collide, we consider the largest one as the survivor, which incorporates the other pixels. The other LCR is considered terminated by the collision.

Output LCR parameters
The parameters outputted by the XYTEventDetector program are as follows:

Sample analysis
To illustrate the functionality and capabilities of our new method, we have performed a sample analysis in 16 cells, comprising 8 rabbit cells and 8 guinea pig cells, measured by two different cameras. The results of our analysis of key LCR parameters are given in S1 Table and S2 Table. The LCRs are also illustrated in 2 movies with LCR borders marked by different colors, S4 Movie and S5 Movie, for rabbit and guinea pig cells, respectively. The S1 and S2 tables also provide all respective program parameters used to analyze LCRs in each cell, thus providing a clue to parameter selection for best performance.
The program output generates individual LCR traces. In Figure 5 shows an example of LCR traces plotted together with the entire LCR ensemble signal (i.e. a formal sum of all LCR traces) and with the whole-cell Ca 2+ signal measured in the entire cell area. Small LCRs appear at the beginning of the pacemaker cycle, then individual LCR grow in size, with the concurrently growing LCR ensemble signal (red thick line).
The LCR propagation is captured by the program as seen in the movies (S4 Movie and S5 Movie), and is evident in Figure 5 by the increase in path area towards the cycle end. These movies show that some of the LCRs behave like small waves that travel a variable distance within the cell before either dying out or merging into the oncoming whole cell Ca 2+ transient, other (typically smaller) LCRs travel very little or not at all. guinea pig (see S1 Table and S2 Table). For example, in rabbit cells these numbers were in the range of 44 to 119, i.e. substantially larger compared to an average of 14 LCRs per cycle previously detected manually in rabbit SANC [12]. This increase was clearly due to detection of the aforementioned 'smaller' LCRs (first bin in histograms in Fig 6).
Also, the LCR duration on average was much shorter, about 13 ms in rabbit, compared  Histograms show substantial contribution of small, short-lived LCRs in both rabbit and guinea pig SAN. Data used in histograms are those presented in S1 Table 1. The LCR size was evaluated by output parameter Max Area scaled to μm 2 .

Noise considerations
Any program detecting events from noise balances between two extremes: 1) Detecting events for sure, i.e. with no false-positive events, but inevitably missing important events closer to the noise; 2) Detecting almost all events, including, small ones, but inevitably catching many false-positive events. Varying the parameters defining sensitivity of LCR detection yields a variety of outputs for LCRs, including the aforementioned extreme cases. Our program is a tool to detect LCRs, rather than a definitive detection algorithm covering all cells and all possible recording conditions. This level of flexibility enables researchers to tune the sensitivity of detection in their studies. Depending on the specific hypothesis to be tested, this detection sensitivity (and the rate of false positive events) can be shifted. For example, if a study is focused on small events, close to the noise, then a larger (but still small) fraction of false positive events is inevitable and must be evaluated and care taken by statistical methods to test the hypothesis. If a study focuses on larger events, then the false-free detection can be substantially improved by decreasing sensitivity level far from the recording noise.
To ensure that we find real LCRs rather than noise, we developed a test to evaluate the rate of false-positive events, which can be used to tune the program for the optimal performance (illustrated in Fig 7). We cut a "control" sub-stack of images within the original data but far from the cell area. This control stack has instrumental noise, but no LCRs. Then we run the same LCR detection algorithm in both the control stack and in the cell area. The program reports the cell area in pixels, making it easy to normalize the number of false positive events found in the control stack area to all events found in the cell area. Our representative measurements of LCR parameters in S1 Table are provided along with an indication of false positive rates in each measurement. With the specific detection parameters, we were able to keep the rate of false positive events low, 0.97% and 6.7% in two rabbit cells and 0% in 6 other cells tested. Thus, this test indicates that the program parameters can be tuned to detect with a high probability true LCR events, rather than noise. Examples are given for SAN from two species, rabbit on the left and guinea pig on the right. The LCRs were found first within the cell area and then using the same procedure in the recording area marked by the yellow square, far from the cell. The number of detected LCRs are shown together with total number of pixels of respective search areas.
It is important to note that our program represents an exploratory tool to detect Ca 2+ release events (or any signaling events) from noise in 2D recording. LCRs should be visible by eye in the movie after filtering and contrasting to ensure that the program will work. It is desirable to have a camera spatial resolution that is close to or surpass the optical limit for a light microscope which is around 250 nm.

Discussion
Here, we have developed a novel and unique set of computer algorithms that allow the automatic detection and analysis of LCRs in spontaneously contracting SA node cells. The capabilities of the programs can be potentially applied to any kind of a spontaneously contracting cell. To achieve this outcome, it has been necessary to overcome many technical challenges, including cell contraction motion artifacts and recording noise. We have also developed new terminology and algorithms for describing LCR path area, collisions, and separations, which represent important, previously undescribed characteristics of LCR behavior.
The entire set of algorithms was developed de novo, and has now been made available for further evaluation, development and customization to test particular hypotheses within the boundaries of specific recording conditions (resolution and signal to noise ratio).
The importance of these new tools for future research in local Ca 2+ events is that they allow rapid, objective, detailed, effective and unbiased detection of LCR events, thus providing new insights into their emergence and biophysical behavior. Although the focus of this paper has been on the development of these intricate methods, we have also clearly demonstrated several important fundamental characteristics of LCRs in SA node cells. Our results show that, in contrast to stereotyped Ca 2+ sparks in ventricular myocytes, the automatically detected LCRs are indeed extremely complex heterogeneous spatiotemporal events that propagate short distances, and often undergo mergers and separations (Fig 4).
Future studies using automatic LCR detection will identify the prevailing physiological conditions favoring LCRs being small, non-propagating spark-like events, versus large, complex widely-propagating events. Based on prior studies using confocal microscopy, the smaller events are expected to occur when cells decrease their pacemaker rate under cholinergic receptor (parasympathetic) stimulation [17], whereas larger and more synchronized events are expected when cells increase their rate in the presence of beta-adrenergic receptor (sympathetic) stimulation [18] [19]. Thus, our new algorithms will be important in uncovering the intricate and detailed mechanisms of autonomic nervous system control of automaticity via LCRs and phosphorylation [20]. Previous studies of LCRs have been incomplete and imprecise due to necessary manual detection and analysis of LCRs along confocal scan lines, or in 2D movies. Our quick test of the new programs detailed above has already revealed substantially more LCRs per cycle, including a subtype of small LCRs (Fig 6), which are hard to reliably resolve manually without clear objective criteria. This presence of small, short lived LCRs has been recently predicted by numerical modeling of Ca 2+ releases generated by RyR clusters of different sizes (including very small clusters) identified by immunofluorescence confocal studies in 3D [21].
The much larger numbers of automatically detected LCR events enabled by the new programs will facilitate highly rigorous statistical analysis of LCR characteristics, which is required to study their stochastic nature in fine detail, and also to link their complex spatiotemporal parameters to the hierarchical structure of release channel distribution found in pacemaker cells [21]. Another important future study will involve comparison of the experimentally detected LCRs with numerically simulated LCRs [21,22], utilizing the same automatic LCR detection programs in both cases. This will validate the recent numerical models of pacemaker cells featuring local Ca 2+ control [21].
The algorithms developed here are also expected to be helpful in studies of physiology and biophysics of Ca 2+ signaling in other cell types, including neurons [23] and muscle cells of non-cardiac origin, in which Ca 2+ signals are typically challenging to detect due to motion artifacts. Indeed, the diverse gamut of LCRs, Ca 2+ wavelets, and abortive Ca 2+ waves demonstrated here in SA node cells also emerge in ventricular myocytes [24] and Purkinje cells [25] under different experimental and physiological conditions, and would be also characterized by the new programs. Specifically, phosphorylation of sarcoplasmic reticulum Ca 2+ cycling proteins unleashes a high-power, rhythmical LCRs in ventricular myocytes linked to arrhythmias and relevant to biological pacemaker design [24], and would benefit from full characterization using our novel programs.

Study Limitations
SANC Analyser relies on the cell being orientated in the horizontal plane for its movement subtraction capabilities to be fully realized. It also relies on the cell being studied not undergoing any translational movement during the time of study, i.e. that cells being studied are well adhered to the dish in which they are being studied. Some SAN cells in particular are highly complex morphologically and have sea-weed like fronds that do not always adhere to the study dish, and as such LCR events in these complex appendages are not possible to quantify using the current program. SANC Analyser however does not rely on cells being long and thin like a typical SAN cell -more uniform shapes such as those typically seen with rounded stem cell-derived cardiomyocytes would also be appropriately handled by the program.
XYTEventDetector detected substantial number of small, short-lived LCRs (Fig 6). Due to the 10 ms sampling interval of our recording cameras, the duration of many small LCRs is 1 sample. While detecting these small LCRs is important in studies of biophysical mechanisms of SANC function, the exact shape and duration of these 1 sample events is imprecise and these data should be treated and interpreted with caution.
It is important to note that increasing just sampling rate of 2D imaging may not necessarily provide a better LCR detection and time profiling. In the present study, we also tested LCR recording with 5 times higher sampling rate (500 vs. 100 frames/s) using a modern PCO.edge 4.2 CMOS camera. However, many LCRs were not detected by our program at this extremely fast sampling rate because to catch LCRs with same success, the signal/noise should also be improved by a factor of 5. Thus, the most reliable detection requires that sampling interval should be close to the characteristic time for signal rise to peak, typically 5-10 ms for sparks and LCRs. However, we want to emphasize that since our algorithms reliably detect and accurately evaluate durations of longer events, this 10 ms limitation is due to temporary constrains of current technologies in recording stronger intracellular calcium signals at higher rates in 2 dimensions, rather than fundamental limits of our algorithms.
Data presented in S1 and S2 Tables were collected form a few cells in a few cycles to illustrate the function of our algorithms and typical parameters of our analysis. The values of specific LCR characteristics in the tables vary greatly and should be interpreted with cation.
Future extended studies are needed for further, more precise evaluation of LCRs in 2 dimensions.

Conclusions
Our study describes in detail novel, specific algorithms for automatic detection and analysis of LCRs occurring in 2D recordings in live, spontaneously beating pacemaker cells. One algorithm tracks points along the midline of the moving/contracting cell, using these points as a coordinate system for affine transform, producing a transformed image series where the cell does not move/contract. Thereafter, a second algorithm produces detailed descriptions of major LCR parameters such as period, signal mass, duration, and propagation path. As the LCRs propagate in live cells, the algorithm identifies complex and elegant splitting and merging behaviors, indicating the importance of locally propagating CICR for the fate of LCRs, and building up a powerful ensemble Ca 2+ signal. While the algorithms were developed to detect LCRs in sinoatrial nodal cells, they have the potential to be used in other applications in biophysics and cell physiology, for example, to detect Ca 2+ wavelets (abortive waves), sparks and embers in muscle cells [2] and Ca 2+ puffs and syntillas in neurons [16].

Funding Statement
The work was supported by the Intramural Research Program of the National Institute on Aging, National Institute of Health. SP was supported the Canadian Institutes of Health Research (MOP12874) and the Natural Sciences and Engineering Research Council (386877).
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.