COSMAS: a lightweight toolbox for cardiac optical mapping analysis

Optical mapping is widely used in experimental cardiology, as it allows visualization of cardiac membrane potential and calcium transients. However, optical mapping measurements from a single heart or cell culture can produce several gigabytes of data, warranting automated computer analysis. Here we present COSMAS, a software toolkit for automated analysis of optical mapping recordings in cardiac preparations. COSMAS generates activation and conduction velocity maps, as well as visualizations of action potential and calcium transient duration, S1-S2 protocol analysis, and alternans mapping. The software is built around our recent ‘comb’ algorithm for segmentation of action potentials and calcium transients, offering excellent performance and high resistance to noise. A core feature of our software is that it is based on scripting as opposed to relying on a graphical user interface for user input. The central role of scripts in the analysis pipeline enables batch processing and promotes reproducibility and transparency in the interpretation of large cardiac data sets. Finally, the code is designed to be easily extended, allowing researchers to add functionality if needed. COSMAS is provided in two languages, Matlab and Python, and is distributed with a user guide and sample scripts, so that accessibility to researchers is maximized.


Algorithm description
The basic version of the presented algorithm is designed for detection of calcium transients in cardiac imaging data. It relies on finding the minima between adjacent calcium transients, reporting these as calcium transient boundaries. We first provide an informal graphical intuition of its principle when applied to calcium imaging. Second, the pseudocode is given. Third, we describe how the algorithm can be adapted for detection of action potentials in electrophysiological recordings or data from membrane potential imaging.
The code implementation of the algorithm (for Matlab and Octave) can be downloaded from the following GitHub address: https://github.com/jtmff/comb. This also includes a script which generates the figures used in this article.

Comb algorithm intuition for calcium imaging
The presented algorithm relies on prior knowledge of the basic cycle length (bcl) of the recording, i.e., how many ms apart are the starting times of calcium transients. The central visual intuition is that a 'comb' is positioned over the calcium fluorescence signal, with its teeth bcl ms apart 1 (Figure 1). The algorithm works in two stages: 1) comb positioning and 2) tooth refinement.
In comb positioning, the comb is slid over the signal, with the first tooth starting at 1 ms, 2 ms, ..., bcl ms 2 . In each position, the mean of signal under the comb teeth is computed, and the minimum over these means is taken as the optimal comb position. Example of three comb positions and corresponding means under the teeth are given in Figure 1. The strength of this approach is that it utilizes information from the whole recording at once and the comb position with minimum average value under its teeth naturally corresponds to locations very close to signal minima.
In the second step, tooth refinement, each near-minimum identified as the location of one of the comb's teeth is taken as a starting point of a local search. The search finds the minimum within ref inementW idth ms around the comb tooth, where ref inementW idth is a user-selected parameter. This step is designed for cases such as cardiac alternans or slightly irregular heartbeat, where the minima are not spaced fully equidistantly. The optimal value of ref inementW idth depends on the degree of irregularity of the recording 3 . 3 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted September 5, 2019. ; https://doi.org/10.1101/757294 doi: bioRxiv preprint

Comb algorithm pseudocode
A pseudocode is given in Algorithm 1, using Matlab-like notation. Therefore, the expression x : y is the vector [x, x + 1, ..., y], and x : y : z is [x, x + y, x + 2y, ..., z], where z is the largest number lower or equal to z. Furthermore, the expression a(b) where a, b are vectors corresponds to those elements in a that are at positions specified in b.
Input: trace (calcium or membrane potential), was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted September 5, 2019. ; https://doi.org/10.1101/757294 doi: bioRxiv preprint

Algorithm adaptation for membrane potential recordings
Whereas in calcium imaging data, the local minima generally correspond to boundaries between calcium transients, this may not hold in optical or electrophysiological recordings of membrane potential, where an action potential is followed by diastolic interval where the signal is relatively flat (Figure 2). A local minimum may lie relatively anywhere within the diastolic interval, which complicates the optimal comb placement. For such case, we suggest the following: 1. Flip the signal so that the peaks of action potentials point down. 2. Find minima using the baseline algorithm (the minima now correspond to peaks of action potentials). 3. Assume a sensible maximum duration of action potential upstroke (e.g., 20-30 ms) and place boundaries this long before the flipped signal minima. Therefore, one obtains signal boundaries that start right before the start of an action potential ( Figure 2).
One exception when this approach could fail is extremely rapid pacing, where the diastolic interval is virtually nonexistent. In such a case, putting an action potential boundary 30 ms before an action potential peak could, in theory, hit the repolarization phase of the previous action potential. However, when the diastolic interval is extremely short, the basic version of the algorithm may be used instead.
Another approach to detection of action potentials may be to use the comb on upside-down flipped vector of derivatives of the signal; in this case, the comb position corresponds to the locations of action potential upstrokes. Figure 2: Membrane potential recording illustration. Sample recording of rat action potentials with boundaries between these highlighted. In larger mammals, the flat diastolic intervals are typically longer than action potentials themselves.

5
. CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

Linear and nonlinear baseline drift, steps, and bleaching
Signal drift is a common finding in experimental recordings. It can result from gradual dye loading (leading to an increase in signal baseline over time), dye bleaching (reducing the signal baseline and amplitude over time), or more prosaic events such as the researcher bumping into the table with the imaging setup, or a droplet of solution passing over the imaged preparation, causing a transient change in signal intensity. The nonlinear drift is especially problematic for simpler algorithms for detection of calcium transients, given that it renders thresholding nearly useless unless the drift can be automatically learned and subtracted.
Both linear and nonlinear drift are processed correctly by the comb algorithm ( Figure  3A,B). In addition, we tested the performance on very abrupt step-like change in the signal baseline; this arises during optogenetic experiments, when the preparation is illuminated by stepwise pulses [2,3]. The algorithm functions well even when the step boundary falls within a calcium transient ( Figure 3C).
Ultimately, we simulated photobleaching by reducing the baseline and amplitude, which did not present an obstacle to the algorithm ( Figure 3D).
One type of drift-like behaviour which might prove problematic to the algorithm is periodic signal perturbation, e.g., arising from from contraction-induced motion artefacts in tissue mapping. In such a case, it is possible that there are signal minima that reflect tissue contraction rather than starts of cellular activation. Given that the frequency of contraction will be identical to the frequency of tissue activation, the artefact-induced minima may be detected as boundaries by the comb algorithm. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted September 5, 2019. ; https://doi.org/10.1101/757294 doi: bioRxiv preprint

Gaussian noise
Gaussian noise is frequently found in optical and electrophysiological recordings. It may confuse thresholding algorithms by producing falsely positive objects of short duration. Approaches based on upstroke velocity may be also confused if the noise amplitude is high enough to mimick a segment of upstroke. The comb algorithm presented here can deal with low to very high amount of noise ( Figure 4A-C), which follows from the global nature of the approach. If the recording is interpreted as the sum of true signal plus zero-mean noise, the noise component under the comb teeth is averaged over the teeth. Given that the noise is zero-mean, the longer the recording, the more teeth are present, making the average converge to zero, reducing the impact of noise on tooth positioning. However, this holds predominantly for initial positioning of the comb; the fine-tuning via local search around the initial tooth positions is vulnerable to noise to a greater degree. Therefore, not using excessively wide local search window is advisable for high-noise scenarios. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted September 5, 2019. ; https://doi.org/10.1101/757294 doi: bioRxiv preprint

Salt-and-pepper noise
Salt-and-pepper noise, also known as impulse noise, presents as sharp and sudden jump in a pixel's intensity, which is typically sparse. This can be handled implicitly by the comb algorithm in many cases, such as Figure 5A. A case of salt-and-pepper which confuses the comb algorithm can be constructed, when a single extremely low value is present ( Figure 5B). This acts as an attractor which forces the comb to put one of its teeth at this location. If such artefacts are present in data, the comb algorithm could be modified to minimize median of values under comb teeth, rather than the mean. Alternatively, median filtering may be applied to the signal, as it is usually an effective countermeasure against salt-and-pepper noise. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted September 5, 2019. ; https://doi.org/10.1101/757294 doi: bioRxiv preprint

Calcium alternans
One phenomenon which is not an artefact, but which may nevertheless prove troublesome, is calcium alternans, the periodic oscillation of large and small calcium transient at rapid pacing [4,1]. The comb algorithm can easily detect separate transients when alternans is present ( Figure 6A), even when it is severe ( Figure 6B). Particularly the latter case of very severe alternans highlights why such data are hard to process e.g., via thresholding approaches, given that the range of thresholds which would detect the smaller calcium transients is narrow (and it may vary between different cells).
Particularly the case of severe alternans ( Figure 6B) is one of motivations for the inclusion of the parameter allowing local search around initial comb tooth position. Given that the large calcium transients are longer than the short ones, a completely rigid comb with teeth 85 ms apart could not find transient boundaries accurately. was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted September 5, 2019. ; https://doi.org/10.1101/757294 doi: bioRxiv preprint

Discussion
Here we have presented a simple 'comb' algorithm for detection of calcium transients or action potentials in recordings of regularly activated cardiac preparations. The algorithm uses prior knowledge of the activation frequency to find all objects of interest at the same time, utilizing global information to overcome potential pitfalls in experimental data. It has the following advantages: • It is resistant to noise and/or drift of signal, unlike simple alternative approaches based on thresholding or detection of rapid upstroke velocity. • It does not need any reference data (like template matching), nor annotations (like machine learning). The independence on precise shape of the segmented object makes it readily applicable for other objects of interest (such as traces of ionic current measurements), as long as these have reasonably obvious minima or maxima. • It is generally very simple, easy to implement in any language, also offering fast computational performance. In addition, it has only a single free parameter which furthermore has a very clear graphical intuition and is thus simple to set.
The fact that the algorithm requires an approximately fixed heart rate is a limitation, but it still covers a range of widely used scenarios in cardiac research. Fixed-rate pacing is often employed in experiments, e.g., to control for unwanted drug-induced changes in heart rate, or when alternans is elicited via the 'S1' protocol. In addition, the 'S1-S2' protocol used to study action potential duration restitution [5] can also be processed automatically. Furthermore, it is often the case that even an unpaced cardiac preparation has a stable-enough heart rate, making it suitable for processing by the comb algorithm. The accurate detection of action potential or calcium transient boundaries then greatly facilitates reliable extraction of biomarkers such as action potential duration, calcium transient duration, calcium transient amplitude (which then may be used to estimate the degree of calcium alternans), or area under a curve (e.g., calcium transient or ionic current measurement). The algorithm may be further adapted to achieve even greater flexibility: • The basic algorithm uses equidistant comb teeth. However, the general approach can be used for any pattern of pacing as long as this is known in advance. E.g., a pacing protocol that employs pacing at gradually increasing rate corresponds to a comb with gradually reducing distance between teeth. Alternatively, alternans measurements typically employ pacing using blocks of fixed-rate stimuli, with the rate increasing between the blocks. This again may be fitted with a comb with corresponding blocks of teeth, where each block has equidistant teeth, but the distance is changed between blocks. • The generation of non-equidistant comb teeth as in the previous point could be in theory automated via approaches detecting signal frequency over time. Additionally, ECG mea-11 . CC-BY 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted September 5, 2019. ; https://doi.org/10.1101/757294 doi: bioRxiv preprint surements are often employed as an additional feature in whole-heart recordings, which can be used to extract heart rate over time, and this could be used to generate an appropriate comb for a parallel imaging recording.