1. Introduction {#sec0001}
===============

Communication between adjacent neurons can be investigated through calcium (Ca) imaging. When this technique is applied, the cells are treated with a Ca-sensitive fluorescent dye. Consequently, events related to the release of calcium such as inter-neuronal communication can be tracked by fluorescence microscopy. This technique has proven its worth because with its help researchers have managed to explain the communication mechanism between muscle cells, heart cells, and neurons [@bib0001].

Besides providing insights into various biological processes, this method can help in the understanding of diseases such as cancer, diabetes, or autoimmune diseases [@bib0002].

The behavior of neuronal circuits depends not only on the properties of its cells, but also on the network topology [@bib0003]. Hence, the knowledge of the connections within a neuronal network and its spatiotemporal properties can provide an insight into gene activation patterns, physiological responses, and pivotal cellular decisions [@bib0004].

While the technique can be applied on a wide range of cells for a variety of purposes, the required quantitative analysis is challenging. Namely, the cell bodies, so-called regions of interest (ROI), need to be determined, the calcium spikes detected, and a network of interacting cells constructed. While this can be done manually, it is a resource- and time-consuming task, especially given the large amount of data that is generally involved in computations. Hence, automation of the analysis is deemed helpful.

Although there are software solutions that allow automating the data analysis from an image stack, these solutions are often based on expensive software packages such as MATLAB [@bib0005], [@bib0006], [@bib0007], [@bib0008] or automate only the cell segmentation and/or the identification of calcium spikes, and they lack the functionality to detect any networks of interacting cells [@bib0009], [@bib0010], [@bib0011], although the identification of such networks is of paramount importance when working with neurons [@bib0012]. Some data types need to be pre-processed before they can be analyzed [@bib0007], [@bib0013], and some tools are handled through scripting [@bib0010], which can cause difficulties for a novice programmer. Furthermore, to finalize the analysis, it often requires a decent amount of time and a detailed knowledge of a system in question [@bib0005].

To aid in the analysis of the data required for calcium imaging, we developed a free and open-source standalone software tool: CALcium IMaging Analyzer (CALIMA, <https://aethelraed.nl/calima>), with the interface GUI shown in [Fig. 1](#fig0001){ref-type="fig"}. To evaluate the reliability, speed, and functionality of CALIMA, we used the experimental data obtained from primary rat cortical neurons [@bib0014] and in-house SH-SY5Y neuroblastoma cultures [@bib0015]. The ROIs and Ca-spike detection outcomes were validated against the results generated by a trained biologist, and several falsely detected or undetected cells and spikes were observed. Moreover, the output of the proposed data analysis method was compared to the output obtained by using alternative segmentation strategies available within various software solutions, viz. Active Contours (FluoroSNNAP) [@bib0013], Otsu thresholding and intensity-based declumping (CellProfiler) [@bib0011], Thresholding and Watershedding (ImageJ) [@bib0016], and spatiotemporal independent component analysis (SIMA) [@bib0010]. The Ca-spikes detection results were compared to the results of an algorithm, which implies comparing spike shapes with a database values to determine the ROI activity (FluoroSNNAP) [@bib0013].Fig. 1Screenshot of CALIMA GUI showing the activity and rate of change for selected ROIs and cross-correlation between selected ROIs (left), a bar graph of the cross-correlation between selected cells (middle), a spatiotemporal map of the neuronal spikes (top right, circle size is proportional to the amount of activity of a ROI, while color indicates the average). The reconstructed neuronal network is overlaid on the footage (bottom right). Red arrows represent detected neuronal connections and their direction of communication.Fig 1

2. Theory {#sec0002}
=========

Several distinct stages are defined to process an image set using CALIMA as shown in [Fig. 2](#fig0002){ref-type="fig"}. As a starting step, the ROIs are detected based on the area covered by cells. Next, the calcium-based activity spikes are identified per each ROI. Lastly, the generated data are analyzed to reconstruct physically possible inter-cell connections. The sections below describe these steps in more detail. The process starts by loading a file or files of calcium imaging experiments. Images can be loaded in such formats as .png, .jpg, .bmp, .gif, or .tiff (no tiff-stacks). Loading video files is enabled using FFmpeg [@bib0017], and the tool supports many formats including .avi, .wmv, and .mp4. The ROIs and analysis results can be exported as images or, in the case when the measured cell brightness and the detected communication network are also available, as .csv files.Fig. 2Flowchart of the CALIMA algorithm.Fig 2

2.1. Cell detection {#sec0003}
-------------------

First, the user selects a representative frame of a video or an image set, which is used as a reference to detect the cells. To make the brighter cells distinguished against the dark background, the contrast is stretched first according to [Eq. (1)](#eqn0001){ref-type="disp-formula"}:$$\mathbf{I}_{x,y} = \frac{\mathbf{O}_{x,y} - \text{min}\left( \mathbf{O} \right)}{\max\left( \mathbf{O} \right) - \text{min}\left( \mathbf{O} \right)}.$$

Here the brightness of a certain pixel of the resulting image **I** at location (*x,y*) is determined as a function of the brightness of the same pixel on the original image **O** and the maximum and minimum brightness values that occur in the image. To filter out the weak inter-neuronal connections and noise, a Gaussian blurring filter is convolved with the frame image as described in [Eq. (2)](#eqn0002){ref-type="disp-formula"} to obtain the blurred version **B** of the image as described in [Eq. (3)](#eqn0003){ref-type="disp-formula"}:$$\mathbf{G}_{x,y}\left( \sigma \right) = \frac{1}{2\pi\sigma^{2}}\text{exp}\left( {- \frac{x^{2} + y^{2}}{2\sigma^{2}}} \right)$$$$\mathbf{B} = \mathbf{I}*\mathbf{G}\left( \sigma_{n} \right)$$

The user-defined standard deviation σ~n~ is used as the 2D standard deviation of the blurring filter. Subsequently, a difference of Gaussians (DOG) filter, the subtraction of the same image modified by two different Gaussian filters, is applied to the resulting image **I** as described in [Eq. (4)](#eqn0004){ref-type="disp-formula"} to separate the fluorescing cells from the darker background [@bib0018]:$$\mathbf{D}\left( {\sigma_{a},\sigma_{b}} \right) = \mathbf{B}*\mathbf{G}\left( \sigma_{a} \right) - \mathbf{B}*\mathbf{G}\left( \sigma_{b} \right)$$where the resulting image **D** is defined as a function of the user-defined standard deviations σ~a~ and σ~b~. In the resulting image, ROIs tend to be of a lighter shade than the background. Therefore, the image is compared to a user-defined brightness threshold *th* to determine, which pixels belong to the ROI in question. On the resulting binary image, a recursive function is used to identify clusters of adjacent pixels belonging to ROIs. To account for possible variations in the input data, of which the most important are zoom level, resolution, brightness and focus, all four parameters of the above-mentioned functions can be adjusted by the user, namely, the standard deviations of the three Gaussian filters and the threshold. For cases where the automatic detection does not suffice, it is also possible to add or remove ROIs manually. To speed up the ROI detection process, the aforementioned filters are separated and a separable convolution algorithm is implemented [@bib0019]. The processing speed of the convolution is further improved by using a multithreaded implementation where computation is distributed across several processor cores.

2.2. Spike detection {#sec0004}
--------------------

The brightness value of a ROI is related to the cell activity, but the base brightness varies greatly per each cell; for example, if two cells are both inactive, one may appear much brighter than another. To account for this effect, the rate of change in brightness per ROI is used for calcium spike detection instead of the brightness itself, as this is not affected by different mean brightness values [@bib0009]. This indicator is then used to detect calcium spikes via the z-scoring algorithm [@bib0020].

For each ROI and for each frame *n*, the rate of change in brightness **r**\[n\] = **b**\[n\] - **b**\[n-1\] is calculated. The mean brightness *µ~l~* and standard deviation *σ~l~* are calculated within a temporal window of a user-defined length *l*, according to [Eqs. (5)](#eqn0005){ref-type="disp-formula"} and [(6)](#eqn0006){ref-type="disp-formula"}. Another user-defined parameter *i* limits the influence of spikes on the calculation of µ~l~ and *σ~l~* according to [Eqs. (7)](#eqn0007){ref-type="disp-formula"} and [(8)](#eqn0008){ref-type="disp-formula"} to increase the robustness of the algorithm [@bib0021].

The z-score is then defined according to [Eq. (9)](#eqn0009){ref-type="disp-formula"}. Such frames, for which ***z*** is greater than the user-defined threshold *th~z~*, and **r** is greater than the user-defined threshold *th~r~*, are detected as calcium spikes in the resulting vector ***s***, as described in [Eq. (10)](#eqn0010){ref-type="disp-formula"}. The user can compensate for spontaneous activity by setting the required minimum number of spikes detected for a ROI to be included in the analysis.$$\mu_{l}\left\lbrack n \right\rbrack = \frac{1}{l}\sum\limits_{k = n - l - 1}^{n - 1}\mathbf{r}\left\lbrack k \right\rbrack$$$$\sigma_{l}\left\lbrack n \right\rbrack = \sqrt{\frac{1}{l}\sum\limits_{k = n - l - 1}^{n - 1}\left( {\mathbf{r}\left\lbrack n \right\rbrack - \mu_{l}\left\lbrack n \right\rbrack} \right)^{2}}$$$$\mathbf{r}_{f}\left\lbrack n \right\rbrack = \left\{ \begin{matrix}
{\mathbf{r}_{i}\left\lbrack n \right\rbrack,\mspace{6mu} if\mspace{6mu} s\left\lbrack n \right\rbrack = 1\mspace{6mu} or\mspace{6mu} if\mspace{6mu}\frac{- \mathbf{r}\left\lbrack n \right\rbrack - 2\mu_{l}\left\lbrack {n - 1} \right\rbrack}{\sigma_{l}\left\lbrack {n - 1} \right\rbrack} > z_{th}} \\
{r\left\lbrack n \right\rbrack,\mspace{6mu} otherwise\mspace{6mu}} \\
\end{matrix} \right.$$$$\text{Where}\;\mathbf{r}_{i}\left\lbrack n \right\rbrack = i \cdot \mathbf{r}\left\lbrack n \right\rbrack + \left( {1 - i} \right) \cdot \mathbf{r}_{f}\left\lbrack {n - 1} \right\rbrack,$$$$\mathbf{z}\left\lbrack n \right\rbrack = \frac{\mathbf{r}\left\lbrack n \right\rbrack - \mu_{l}\left\lbrack {n - 1} \right\rbrack}{\sigma_{l}\left\lbrack {n - 1} \right\rbrack}$$$$\mathbf{s}\left\lbrack n \right\rbrack = \left\{ \begin{matrix}
{1,\mspace{6mu} if\mspace{6mu} r\left\lbrack n \right\rbrack > th_{r},z\left\lbrack n \right\rbrack > th_{z}} \\
{0,\mspace{6mu} otherwise\mspace{6mu}} \\
\end{matrix} \right.\mspace{6mu}$$

2.3. Correlation maps {#sec0005}
---------------------

After identification of the ROI activity and peaks, the Pearson cross-correlation x~τ~ for a lag τ between the ROIs is determined according to [Eq. (11)](#eqn0011){ref-type="disp-formula"}, where µ~f~ and µ~g~ denote the mean values of ***f*** and ***g*** according to [Eq. (12)](#eqn0012){ref-type="disp-formula"}. The variables **f** and **g** can be either the brightness **b**, the brightness rate of change **r**, or the spike vector **s** of ROIs as described in [Section 2.2](#sec0004){ref-type="sec"}.$$x_{\tau}\left\lbrack n \right\rbrack = \sum\frac{\mathbf{f}\left\lbrack {n - \tau} \right\rbrack - \mu_{f}}{\left| \left| {\mathbf{f} - \mu_{f}} \right| \right|} \cdot \frac{\mathbf{g}\left\lbrack n \right\rbrack - \mu_{g}}{\left| \left| {\mathbf{g} - \mu_{g}} \right| \right|};$$$$\mu_{y} = \frac{1}{\left| \mathbf{y} \right|}\mspace{6mu}\sum\mathbf{y}\left\lbrack n \right\rbrack,\mspace{6mu}\mathbf{y} \in \mathbf{f},\mathbf{g}$$

The detected peaks, the rate of change, and the activity are shown directly in a graph for separate ROIs ([Fig. 1](#fig0001){ref-type="fig"}, left), whereas *x~τ~* between all combinations of selected cells can be presented in a form of bar graph for −5≤ τ ≤ 5, a correlation matrix or a correlation map ([Fig. 1](#fig0001){ref-type="fig"}, center). The correlation matrix displays *x~τ~* so that the cross-correlation between the first and second selected cells is at the position (1, 2), between the first and third at (1, 3), etc., with the color of each pixel selected from a provided colormap based on the value of *x~τ~* (leftmost color for *x~τ~* = −1 up to the rightmost color for *x~τ~* = 1). This provides the user with information regarding the overall correlation between all selected ROIs.

The correlation map shows a spatial map of all detected ROIs colored depending on *x~τ~*. The single-cell correlation map displays through colors which areas show activity similar to a selected cell, providing a clear view of which ROIs have the calcium fluorescence behavior similar to the one of the selected cell. Furthermore, a heat map can be created to see which ROIs and/or areas of the cell culture are the most active. This can be seen in a spatial map of all detected ROIs colored depending on the value of **b, r** or **s** (the leftmost colormap value for 0, the rightmost color for the maximum value of **b, r** or **s**).

Another map is drawn to provide information regarding the spatiotemporal characteristics of the ROIs ([Fig. 1](#fig0001){ref-type="fig"}, top right). It displays a spatial map of all ROIs, for which calcium spikes were detected, drawn as circles. The diameter of these circles is linearly dependent on the number of detected spikes for that ROI, whereas the color depends on the average spike frame *n~avg~* = ∑*s*\[*n*\] · *n*, ranging from red for the first frame to blue for the last one.

2.4. Network analysis {#sec0006}
---------------------

As a final step, the inter-neuronal connections are reconstructed. The cross-correlation between ROIs is used as a measure of the likeliness of the connection presence. A minimum correlation factor, ranging from 0.0 to 1.0, is given by the user so that only those areas that produce similar spiking patterns are considered as being functionally connected. Furthermore, a maximum connection length of the neurites can be given as input for anatomical connections (along with the image width). In that case, only connections shorter than the upper limit are considered. This feature allows the user to either include or exclude anatomical (physical) connections of the network in the analysis and provides a means to separate between functional and anatomical connections. Functional network is obtained by setting connection length to max value so that the network is based using only the correlation of activity time series. A stricter anatomical network can also be obtained by restricting the connection length value between two ROI preventing the program from reconstructing networks with connections of such length that they are not physically possible. The same accounts for the maximum delay between an incoming spike and the cellular response, which can be tuned after setting the footage frames per second (FPS). The connections are presented on a map of the cell culture ([Fig. 1](#fig0001){ref-type="fig"}, bottom right), showing which ROIs are likely to be connected in physiological response. Lastly, the produced graphs and representations of the connections between the cells can be exported as a csv file for further external analysis.

3. Methods {#sec0007}
==========

The performance of CALIMA was measured on three different datasets, two of which were produced in-house, as described in [Section 3.1](#sec0008){ref-type="sec"}. Both cell detection ([Section 3.2](#sec0009){ref-type="sec"}) and combined spike detection were tested and benchmarked against several alternative cell segmentation and spike detection strategies as described in [Section 3.3](#sec0010){ref-type="sec"}.

3.1. Cell culture and PDMS preparation {#sec0008}
--------------------------------------

The neuroblastoma cell line SH-SY5Y [@bib0022] was used to characterize the performance of CALIMA software on polydimethylsiloxane (PDMS) substrates. PDMS and a cross linking agent were mixed in a cup in a ratio of 10:1 and degassed for approximately 20 min. After degassing, the liquid PDMS was poured onto a glass slide directly from the cup. Thereafter, the PDMS was cured on a hotplate at 95 °C for 10 min and coated with 50 µg/ml fibronectin solution in phosphate-buffered saline (PBS) for 45 min prior to cell seeding. SH-SY5Y neuroblastomas were cultured in DMEM/F-12 media (1:1) supplemented with 10% fetal bovine serum (FBS) and 1% pen/strep, and they were grown in an incubator at 37 °C, 5% CO~2~. When cell confluency was reached, trypsin (x1) was used to harvest the cells and thereafter, was centrifuged at 900 rpm for 5 min. A concentration of 200,000 cells/ml was used throughout the seeding experiments. Once the cells were seeded, they were left to adhere for 24 h following exposure to retinoic acid, thereafter, for 3 days at 10 µM in DMEM/F-12 media to differentiate the cells into neurons according to the manufacturer protocol [@bib0023]. Following 3 days of differentiation, a standard culture medium using DMEM/F-12 medium was added. For calcium staining experiments, the Fluo-4 Calcium Imaging kit (molecular probes, life technologies) was used according to the manufacturer\'s protocol and as described by [@bib0024]. Video footage of the samples were recorded using EVOS FL microscope (ThermoFisher, Waltam, MA, USA) at 1 frame per 10 s, and for 10 min with frames of 1280 by 960 pixels. When neurons maturate they become electrically active and start to blink in synchrony (coherent neuronal activity patterns) [@bib0025], a phenomenon that can be recorded by calcium imaging experiments [@bib0026]. One dataset (Dataset 1) was taken from an early differentiated neuron culture (10 days in culture), and it did not exhibit this phenomenon. Another dataset (Dataset 2) was taken from a late differentiated neuron culture (28 days in culture); for this one the effect of synchronous neuron firing was observed.

In addition to the data obtained as described, we analyzed a validated example, Dataset 3, which was provided by the FluoroSNNAP software (spontaneous activity of cultured primary rat cortical neurons transduced with AAV2-GCaMP6, 696 by 520 pixels, 10 Hz measurements for 2 min, downloaded from [@bib0014]).

As a result of the experimental procedures, it was observed that Dataset 1 consisted of zoomed-out densely packed barely active cells of various types and some cellular debris and Dataset 2 of zoomed-in densely packed active cells, likely connected in a network. Dataset 3 showed zoomed-in cells communicating in a network [@bib0013] located at a distance from each other. Examples of single frames from each dataset are shown in [Fig. A.1](#fig0005){ref-type="fig"}.

3.2. Cell detection {#sec0009}
-------------------

The aforementioned neuronal cultures were with CALIMA for estimation purposes. The parameters used for cell detection varied between the datasets, the detailed description is provided in [Table B.1](#tbl0002){ref-type="table"} for reference.

The cells as detected by the softwares were validated against a manual segmentation done by a trained biologist. It was checked how many of the total number of areas that were found manually (P) were detected successfully. An overlap between a detected ROIs and a manually detected cell were designated True Positives (TP). ROIs that overlapped with multiple cells were not included in the TPs and were counted separately. Likewise, if ROIs did not match a manually detected cell or if multiple ROIs were assigned to one cell, the redundant ROIs were designated False Positives (FP). Cells not detected manually were included in the False Negatives (FN).

Using this classification, the correctness of the detection algorithm was assessed using three different measures: the true positive rate (TPR) or sensitivity by [Eq. (13)](#eqn0013){ref-type="disp-formula"}, the positive predictive value (PPV) using [Eq. (14)](#eqn0014){ref-type="disp-formula"}, and the recall defined as per [Eq. (15)](#eqn0015){ref-type="disp-formula"}.$$\text{TPR} = \text{Sensitivity} = \frac{\text{TP}}{P}$$$$\text{PPV} = \frac{\text{TP}}{\text{TP} + \text{FP}}$$$$\text{recall} = \frac{\text{TP}}{\text{TP} + \text{FN}}$$

3.3. Benchmark ROI detection {#sec0010}
----------------------------

The performance of the cell detection on the same datasets was benchmarked against the alternative segmentation strategies. The active contours algorithm (FluoroSNNAP [@bib0013]) was studied for its ability to capture the crude shapes. The initial estimate made with the use of this algorithm did not require any user input.

Otsu thresholding (CellPRofiler [@bib0011]) was validated because of its ability to separate foreground and background pixels. A global threshold was chosen accepting ROIs diameter from 3 to 150 pixel units (pu) (Dataset 1), 25--250 pu (Dataset 2), and 10--45 pu (Dataset 3). A three-classes threshold was chosen, middle intensity pixels were assigned to the background and the threshold correction values were set to 1.3, 0.6, and 1.0 respectively for three datasets in question. In Dataset 1, a 50-pixel window adaptive threshold was chosen, the other datasets were processed with a global threshold. Merged cells were detected and separated based on pixel intensity.

We tested a combination of a thresholding and water shedding (ImageJ version 1.51k [@bib0016]) because of its speed. The image was converted to 8 bit to enable water shedding later on, the threshold was set to 15, 30, and 25 respectively (Datasets 1, 2, and 3). After water shedding was applied to separate the ROIs, only ROIs in the range of 50--500 pixels (Datasets 1), more than 80 pixels [(2)](#eqn0002){ref-type="disp-formula"}, and more than 10 pixels [(3)](#eqn0003){ref-type="disp-formula"} were considered by using the Analyze Particles option.

Lastly, we performed the spatiotemporal independent component analysis (SIMA [@bib0010], [@bib0027] because this algorithm combines both spatial and temporal information to determine the cell\'s locations. The number of principal components used was chosen to be 500 for Dataset 1, 160 for Dataset 2, and 100 for Dataset 3, which is slightly higher than the estimated number of active cells [@bib0006]. The relative weight between the spatial and temporal independent component analysis was set to 0.85 (dominant spatial information) regarding the not-so-active Dataset 1 and the sparsely populated Dataset 3. In Dataset 2, this weight was set to 0.15. Cells overlapping more than 50% were merged and the ROIs were smoothed and filtered using the standard implementation of a method SparseROIsFromMasks proposed in SIMA [@bib0010].

All data were processed using an HP EliteBook 850 G2 with 16 GB physical memory running a 64-bit Windows 10 operating system to compare the results obtained using each method. The time required for estimating the segmentation was recorded for each segmentation strategy. Setting the parameters required about the same amount of time for each method; thus, it was ignored.

3.4. Benchmark spike detection and network analysis {#sec0011}
---------------------------------------------------

The spikes detected by CALIMA were validated against the ones obtained through manual detection and compared to a waveform-based detection algorithm (FluoroSNNAP [@bib0013]). Dataset 3 was used to validate both the spike detection and the network reconstruction, as the optimal settings for spike detection were known and it had already been validated to contain neurons communicating in a network [@bib0013]. To avoid considering changes in the segmentations, the segmentation mask defined manually was used in both software tools before validating the spikes and network. The first 15 ROIs (indexes assigned by FluoroSNNAP) were used to check the spike detection.

The spikes were manually checked regarding height and shape, and sets of the cells matching the criteria were stored. In CALIMA, the spike detection was based on the mean values per ROI and the parameters were set to 10, 5, 0.50, and 0.100 for *l, z~th~, m*, and roc~min~ respectively. The relatively high value of *z~th~* is a direct result of the low amount of noise in Dataset 3. In FluoroSNNAP, the settings recommended by the authors were used: choosing a 10 s window and using the lower 50% of samples to determine the ∆ F/ F~0~, using their standard library to match and detect the Calcium events.

To reconstruct a network, CALIMA used the Pearson coefficients in combination with the inter-cellular distance. ROIs with a cross-correlation coefficient above 0.25, a distance lower than 208 µm (a full image width estimated at 1040 µm), and a response time of 500 ms or less to each other were considered in this case. FluoroSNNAP used a partial-correlation method with an alpha level of 0.001, as this method represents the inter-cell connections, which are the closest ones among all methods offered by the software [@bib0013].

4. Results {#sec0012}
==========

Using the methods as described in the previous section, the datasets described in [Section 3.1](#sec0008){ref-type="sec"} were processed using CALIMA and other alternative software tools. The results are shown in [Table 1](#tbl0001){ref-type="table"} and described in the following sections.Table 1Comparison of the assets and performance of different software tools on the given datasets.Table 1Tool nameCALIMAFluoroSNNAPCellProfilerImageJSIMASegmentation strategyDifference of GaussiansActive contoursOtsu-thresholdingThresholding & watersheddingSpatiotemporal ICASpike detectionYesYesYesNoNoNetwork analysisYesYesNoNoNoSensitivity range \[%\]80--8621--6551--8539--775--28Positive prediction value range \[%\]60--8015--9327--844--873.8--49Recall \[%\]88--9839--6979--9941--9441--100Merged cells (Dataset 1,2,3)45,14,1413,71,40,72,66,87,4159,143,10Mean raw calculation time for segmentation \[s\]5115401139

4.1. Cell detection {#sec0013}
-------------------

The results of the cell detection accuracy across the three datasets are given in [Table 1](#tbl0001){ref-type="table"} (a more detailed analysis per dataset can be found in [Table C.1](#tbl0003){ref-type="table"}.), and the ROIs found in dataset 3 are represented in [Fig. 3](#fig0003){ref-type="fig"} (ROIs found for datasets 1 and 2 can be found in [Figs. A.2](#fig0006){ref-type="fig"} and [A.3](#fig0007){ref-type="fig"}). CALIMA demonstrates high scores for all three datasets with a sensitivity of 80--86% and a PPV of 60--80%. Though CALIMA performs well on the densely packed Dataset 2, the detection method comes at a cost. CALIMA detects only a part of the ROIs rather than the entire shape. Nevertheless, the raw calculation time is small, which helps a user to test different settings and find the optimum quickly.Fig. 3Detected ROIs for Dataset 3 using various segmentation strategies. ROIs, manual (a), ROIs, CALIMA (b), ROIs, Active Contours (FSNNAP) (c), ROIs, Otsu Thresholding + brightness-based declumping (CellProfiler) (d), ROIs, Thresholding + Watershedding (ImageJ) (e), ROIs, spatiotemporal Independent component analysis (SIMA) (f). The ROIs found manually are depicted in blue, the true positives in green, and the false positives in orange. The ROIs containing multiple cells are shown in purple and the false negatives in red. Note that an exception is made for the stICA method, as it allows regions to overlap. To avoid confusion, the segmentation outlines found by this method are shown as well.Fig 3

The active contours algorithm appears to perform better in detecting the full extent of the ROIs but can detect only a smaller number of them. The Otsu-thresholding performs well in terms of sensitivity, speed, and determining the ROIs\' actual shape to the cost of obtaining some redundant ROIs at the borders. Furthermore, when the cells are packed too densely (Dataset 2), the automated declumping has problems. The Thresholding and Watershedding strategy has difficulties with densely packed cells as well; however, it is quite sensitive in other cases. The separation of the cells appears to work only for small clusters of cells, as indicated by a large number of merged cells in Dataset 2. The stICA method appears to scatter its ROIs, resulting in obtaining many isolated small and overlapping regions. However, stICA shows the best performance among all considered algorithms in terms of finding the ROI outlines, and it is able to detect shapes as small as neurites. An option to automatically detect cells of different shapes will be considered as well for CALIMA in the future [@bib0028].

4.2. Spike detection and network analysis {#sec0014}
-----------------------------------------

Both the z-score algorithm and the database comparison demonstrated appropriate scores in terms of sensitivity; FluoroSNNAP was able to identify 96.5% of the peaks on average (ranging between 90 -- 100% per cell) and CALIMA identified 98.5% (84--100%). The database comparison for FluoroSNNAP had 32 false positives spread over 9 cells, while CALIMA discovered 14 in 6 cells. Most false positives of the database comparison (16) were found in a single cell indicating the inhibition of a brightness drift and a large number of non-AP spikes. The Z-score\'s false positives were spread more evenly across the cells. The number of false negatives of FluoroSNNAP was 18 across 9 cells and for CALIMA 19 across 9 cells. The Z-score algorithm missed most of peaks (6) in a single cell; this observation indicates fast spiking activity.

As can be seen from [Fig. 2](#fig0002){ref-type="fig"}a and b, the networks reconstructed using the 2 different methods differ significantly from each other ([Fig. 4](#fig0004){ref-type="fig"}c and d). In addition, the binary adjacency matrices also differ, with much less ROIs being linked with one another. The most typical difference is that the connections in CALIMA are shorter, while FluoroSNNAP allows observing connections over larger distances. Since a typical neurite length is in the order of 100′s µm [@bib0029], the user-defined maximum length of neuronal connections allows CALIMA to filter out these longer connections, resulting in a physiologically more sensible communication network including both functional and anatomical parameters.Fig. 4The reconstructed networks and binary adjacency matrices for dataset 3. Reconstructed network, ROI Pearson cross-correlation + distance filter (CALIMA) (a). Reconstructed network, ROI partial correlation (FSNNAP) (b).The ROIs are shown as dots and the connections as lines. Binary adjacency matrix of CALIMA (c) and FSNNAP (d). The matrices represent the ROIs, those that are connected (in white) and those that are not.Fig 4

5. Discussion {#sec0015}
=============

CALIMA was tested on three real-life datasets; consequently, it demonstrated a high sensitivity for the ROIs and allowed performing quick and highly confident cell detection on a large variety of data. When compared with other segmentation strategies, CALIMA performed better in most cases (see [Tables 1](#tbl0001){ref-type="table"} and [C1](#tbl0003){ref-type="table"}).

The automatic mode of the active contours algorithm (FluoroSNNAP) allows detecting only a small number of ROIs. This algorithm performs better if the user provides growth regions. Otsu\'s method (CellProfiler) enables multiple-class thresholds that can be combined with an adaptive window. Though the algorithm detects most cells correctly as well as identifies the ROI shape in detail, it does not perform well when cells are densely packed. As a result, the algorithm might merge a number of cells. The 'declumping' algorithm would have to be fine tuned by the user to correct this. The thresholding strategy (ImageJ) works well to detect large clusters of cells, but it is inaccurate when it makes segmentation of these clusters into separate cells as a next step. The cell irregular shapes are the likely cause of this. The algorithm has a problem with detecting cells near the cellular debris, which is because these cells are merged in a layer of moderate brightness, which cannot easily be handled by using a single threshold. Lastly, the stICA method (used in SIMA) does not show good results in terms of the sensitivity of its performance, or the redundant areas it discovers. Dataset 1 does not show sufficient activity for the algorithm to find the separate cells, and in Dataset 2 the cells are packed too densely. Therefore, it was concluded that stICA, while being able to detect the cellular shape up to the detail level of the neurites, is not suitable to analyze such cultures as ones presented in Datasets 1 and 2 and requires a low culture cell density to successfully analyze the footage, which is not the case for CALIMA.

Using the correct parameters, CALIMA\'s Z-score algorithm is able to detect a large number of calcium-spikes correctly. A sample set of 15 cells demonstrated a sensitivity range of 86--100%. These values are comparable to the results observed for the active contours algorithm (FluoroSNNAP). However, it was concluded that the FluoroSNNAP algorithm exhibits difficulties when separating between action potentials and spike-mimicking noise. Meanwhile, the z-score algorithm of CALIMA has problems with identifying all spikes in cells showing near-constant activity. Although these are rare events, the excess of spikes disrupts the estimation of the standard deviation that is required by CALIMA to detect these spikes correctly.

The networks reconstructed by both CALIMA and FluoroSNNAP differ at various points; the lengths of the connections as well as the actual connections differ. This is probably a result of the different reconstruction strategies implemented. CALIMA focuses on the physical limitations of the cells when identifying the network structure and tries to establish the actual connections. FluoroSNNAP rather focuses on the transfer of information within the network and links all blinking cells together. As a result, CALIMA reconstructed network is a physiologically more sensible communication network including both functional (similar synchronous blinking pattern) and anatomical parameters (physical distance between neuron).

6. Conclusion {#sec0016}
=============

Although calcium imaging is widely used in biological research, there are relatively few free and automated tools available for performing such analysis. These tools often enable only part of the analysis. Here we proposed the standalone free and open-source software tool CALIMA, and we demonstrated that it could serve as a valuable addition to the field of calcium imaging. Starting from crude calcium imaging data, CALIMA can detect cells, determine their activity, and reconstruct a communication network within minutes.

As it was proved that CALIMA is able to detect ROIs, find Ca-peaks, and establish a communication network, we consider CALIMA to be a valuable addition to the existing set of Ca-analysis tools. Because the DOG algorithm detects only a small part of most ROIs, future versions of CALIMA will contain a region-growing step to complete identification process for those areas for which only a part has been found. An option to automatically detect cells of different shapes will be considered as well. Furthermore, the algorithm detecting the physical network will be combined with information transfer theories to study the communication between cells in detail.

Conflict of interest {#sec0017}
====================

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. This statement has also been added at the end of the manuscript.

Appendix A. The datasets {#sec0020}
========================

[Figs. A.1](#fig0005){ref-type="fig"}--[A.3](#fig0007){ref-type="fig"}Fig. A.1The averaged frame of each dataset. Dataset 1: A slightly active culture containing a mix of various cell types, different cell sizes and aggregates. The bar represents 400 µm (a). Dataset 2: An active culture consisting of densely-packed cells. The bar represents 200 µm (b). Dataset 3: An active culture consisting of separated cells. The bar represents 250 µm(estimated) (c).Fig A1Fig. A.2Detected ROIs for Dataset 1 using various segmentation strategies. The bar represents a width of 400 µm. ROIs, manual (a), ROIs, CALIMA (b), ROIs, Active Contours (FSNNAP) (c), ROIs, Otsu Thresholding + brightness-based declumping (CellProfiler) (d), ROIs, Thresholding + Watershedding (ImageJ) (e), ROIs, spatiotemporal Independent component analysis (SIMA) (f). The ROIs found manually are depicted in blue, the true positives in green, and the false positives in orange. The ROIs containing multiple cells are shown in purple and the false negatives in red. Note that an exception is made for the stICA method, as it allows regions to overlap. To avoid confusion, the segmentation outlines found by this method are shown as well.Fig A2Fig. A.3Detected ROIs for Dataset 2 using various segmentation strategies. ROIs, manual (a), ROIs, CALIMA (b), ROIs, Active Contours (FSNNAP) (c), ROIs, Otsu Thresholding + brightness-based declumping (CellProfiler) (d), ROIs, Thresholding + Watershedding (ImageJ) (e), ROIs, spatiotemporal Independent component analysis (SIMA) (f). The ROIs found manually are depicted in blue, the true positives in green, and the false positives in orange. The ROIs containing multiple cells are shown in purple and the false negatives in red. Note that an exception is made for the stICA method, as it allows regions to overlap. To avoid confusion, the segmentation outlines found by this method are shown as well. The reconstructed networks and binary adjacency matrices for dataset 2 are also shown. Reconstructed network, ROI Pearson cross-correlation + distance filter (CALIMA) (g). Reconstructed network, ROI partial correlation (FSNNAP) (h).Fig A3

Appendix B. Parameters used for the cell detection algorithm {#sec0021}
============================================================

[Table B.1](#tbl0002){ref-type="table"}.Table B.1The variables used for cell detection for the various datasets.Table B1σnσaσbthDataset 14140.001Dataset 26.2130.0001Dataset 36.512.70.0002

Appendix C. Results of the cell detection accuracy for different software tools {#sec0022}
===============================================================================

[Table C.1](#tbl0003){ref-type="table"}.Table C.1The total number of cells found in the dataset, the number of true positives, false negatives, merged ROIs, and the number false positive ROIs per segmentation strategy. The detection sensitivity, positive predictive value and recall are shown in the last three columns. The best score per dataset is marked bold.Table C1Dataset1TotalTPFNMergedFPDetec RatePPVRecallDifference of Gaussians (CALIMA)48841726452770.850.60.94Active contours (Fluorosnnap)48818429113460.380.80.39Otsu-thresholding + declumping (CellProfiler)488387101010530.790.270.79Thresholding + watershedding (ImageJ)48828519761870.60.390.59SpatioTemporal ICA (SIMA)4881341951592900.270.320.41Dataset2TotalTPFNMergedFPDetec RatePPVRecallDifference of Gaussians (CALIMA)1501201614160.80.880.88Active contours (Fluorosnnap)1503148711750.210.150.39Otsu-thresholding + declumping (CellProfiler)15072177150.50.830.99Thresholding + watershedding (ImageJ)15059487160.390.790.94SpatioTemporal ICA (SIMA)150701431750.050.041Dataset3TotalTPFNMergedFPDetec RatePPVRecallDifference of Gaussians (CALIMA)7560114120.80.830.98Active contours (Fluorosnnap)754922440.650.920.69Otsu-thresholding + declumping (CellProfiler)756456430.850.60.93Thresholding + watershedding (ImageJ)755813490.770.870.82SpatioTemporal ICA (SIMA)7549125510.650.490.98

Appendix D. Supplementary materials {#sec0019}
===================================

Image, application 1

This work was financially supported by the European Research Council (ERC), MESOTAS grant no. 280281, MESOTAS-SIEVE, grant no. 713732 and ZONMW Off road grant no. 40-08125-98-160. The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material associated with this article can be found, in the online version, at doi:[10.1016/j.cmpb.2019.104991](https://doi.org/10.1016/j.cmpb.2019.104991){#interref0003}.

[^1]: Both authors contributed equally to this work.
