NeuroCa: integrated framework for systematic analysis of spatiotemporal neuronal activity patterns from large-scale optical recording data

Abstract. Optical recording facilitates monitoring the activity of a large neural network at the cellular scale, but the analysis and interpretation of the collected data remain challenging. Here, we present a MATLAB-based toolbox, named NeuroCa, for the automated processing and quantitative analysis of large-scale calcium imaging data. Our tool includes several computational algorithms to extract the calcium spike trains of individual neurons from the calcium imaging data in an automatic fashion. Two algorithms were developed to decompose the imaging data into the activity of individual cells and subsequently detect calcium spikes from each neuronal signal. Applying our method to dense networks in dissociated cultures, we were able to obtain the calcium spike trains of ∼1000 neurons in a few minutes. Further analyses using these data permitted the quantification of neuronal responses to chemical stimuli as well as functional mapping of spatiotemporal patterns in neuronal firing within the spontaneous, synchronous activity of a large network. These results demonstrate that our method not only automates time-consuming, labor-intensive tasks in the analysis of neural data obtained using optical recording techniques but also provides a systematic way to visualize and quantify the collective dynamics of a network in terms of its cellular elements.


Introduction
2][3] Therefore, deciphering the spatial and temporal patterns of neuronal activity in a large population is essential to understand the operating principles of neural circuits. 4Optical recording using activitydependent fluorescence sensors, particularly calcium indicators, is a powerful method due to its superior resolution in space by comparison with electrophysiological approaches. 5Nowadays, it becomes feasible to simultaneously capture the activity of hundreds and thousands of individual cells from dissociated cultures, 6,7 tissue slices, [8][9][10] or the brain of living animals, [11][12][13] providing the glimpses of collective neural dynamics.Although the optical acquisition of neural activity data has been remarkably advanced owing to the efforts in developing genetic sensors [14][15][16] and imaging techniques, 17,18 it remains challenging to analyze the collected image datasets due to their massive size and complexity.
The analysis of calcium imaging data requires spike sorting to isolate the signals of individual cells and quantitative representation of their activity patterns.However, computational tools are often unavailable that supports automated processing and quantitative analyses of the imaging data.As a result, typical approaches have relied on the manual annotation of cells [19][20][21] and qualitative comparison of their fluorescent signals 13,22,23 despite the consumption of considerable time and human labor.Recently, several computational methods were suggested to automate cell identification [24][25][26][27] and spike inference [28][29][30][31][32] from the imaging data.However, additional programming still requires for integrating and optimizing these algorithms for practical use in data analysis.Additionally, even though the activities of individual neurons are identified, visualizing and deriving meaningful features of their ensemble activity are still challenging due to the deficiency of suitable and standardized methods. 33,34Thus, developing a comprehensive software package for automated processing and quantitative analyses of calcium imaging data would be beneficial for large-scale neurophysiological studies.
In this work, we have developed a new open-source, standalone toolbox called NeuroCa.This program includes our new algorithms for cell identification and spike detection from the calcium imaging data, and also allows the quantification and visualization of neuronal activity patterns in a large network (Fig. 1).To detect individual cells, we applied a morphological feature extraction method based on the circular Hough transform to the image set.Our method enabled isolating the regions of individual cell bodies with high efficiency, resulting in the fluorescence signals of each cell as "multi-channel neural data."To extract the calcium spikes of each neuron, we devised a procedure that corrected the background trend of cellular signals by using curve fitting and subsequently detected calcium spikes by using deconvolution with data-driven kernels.Using dissociated cultures of cortical neurons, we demonstrated that the reconstructed spike trains of individual neurons could be utilized not only to calculate the quantitative measures of cellular activity, but also to estimate the functional connectivity of the neural circuits.In particular, we attempted to exploit neuronal calcium spike trains to infer the firing patterns in the synchronous activity of a large network for demonstrating the utility of our approach to dissect the spatial and temporal organization of neuronal circuits at cellular scale.

Cell Culture
Cortical neurons were cultured to construct functionally active neural networks in vitro.The cortical tissues were dissected from E18 SD rats (Koatech, Republic of Korea) and immersed into HBSS (14175, Gibco, California).After dissociating the tissue into the single cells, we centrifuged the suspension at 1000 rpm for 2 min.Supernatant was then gently removed, and the plating medium [Neurobasal medium (21103, Gibco, California) supplemented with B27 (17504-044, Gibco, California), 2 mM GlutaMAX-1 (35050, Gibco, California, 12.5 μM L-glutamate (L-Glutamic acid, nonanimal source, G8415, Sigma, Missori, and 1% (v/v) penicillin-streptomycin (15140, Gibco, California)] was re-filled to suspend cells.Cells were cultured on the substrate with the density of ∼1000 cells∕mm 2 , and half of medium was changed with the maintenance medium (as same as the plating medium without L-glutamate) twice a week.All procedures of cultivation were performed according to the approved animal use protocols of the KAIST Institutional Animal Care and Use Committee.
To introduce the calcium indicator to the samples, the stock solution was added to bACSF and gently mixed by pipetting (the final concentration of OGB-1 was 2.5 mM).After aspirating the whole medium in the culture, we added the diluted indicator solution to the sample for 30 min of dye loading.Next, the solution was aspirated again and washed with fresh bACSF two or three times.The sample was soaked in fresh bACSF for 30 min again for stabilization and then ready to be measured.
The imaging setup was composed of an upright microscope (BX51, Olympus, Japan) with a light-emitting diode (LED) source (SOLA SM, Lumencor, Oregon), a camera (sCMOS Neo, Andor Technology, UK), and a heating plate with a temperature controller (TC01, Multichannel Systems, Germany) for maintaining 37°C during imaging.We acquired images with the frame rate of about 32 Hz and the field-of-view (FOV) of about 800 μm × 700 μm.
The computer for acquiring real-time image data was composed of Intel® Core i5-2400 processor (3.1 GHz; motherboard: Asus P8H67) supporting SATA 6 GB∕s and 16 GB of memory.For achieving fast and stable data acquisition, we connected solid-state drives (SSDs; 840 Pro, Samsung, Republic of Korea).In all calcium imaging experiments, image data were collected via custom-made software based on Andor SDK3 (Andor Technology, UK).

Chemical Stimulation
Several agonists of neurotransmitter receptors, such as N-methyl-D-aspartic acid (NMDA; M3262, Sigma, Missouri), α-amino-3hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA; Asc-130, Ascent Scientific), and γ-aminobutyric acid (GABA; A2129, Sigma, Missouri), and an antagonist of GABAergic receptor (bicuculline; Tocris), were used to change neural activity.NMDA, GABA, and bicuculline were dissolved into bACSF with 20 μM of concentration and AMPA was dissolved into bACSF with 40 μM of concentration.After baseline recording with the culture immersed in pure bACSF, half of the bath solution was replaced with the same amount of drug-containing solution; therefore, the final concentration of NMDA, GABA, and bicuculline was 10 μM, and that of AMPA was 20 μM.All the bath solution was aspired for washing out, and fresh bACSF heated up to 37°C in advance was filled in the culture.

Cell Body Detection
To correct the imbalance illumination and enhance the contrast level, top-hat filtering and contrast adjustment were applied to the image data.Next, circular Hough transform (CHT) was applied to the corrected image to detect all circular elements (here, somata) as independent regions of interest (MATLAB function: imfindcircles).This step gave rise to the center and radius of each region of interest (ROI) that is available to construct the binary mask image.We used radius 1 to fill the ROI area instead of the original value to avoid the overlapping between different ROIs.Furthermore, to prevent missing cells that are not activated at the very first snapshot, we repeated this procedure in the first 200 images of a sequence (if the total number of images were less than 200, all the images were used) and overlapped each mask into one by "OR" operation.Using this final mask, we could obtain the spatial information (center and radius) of each circular element that outlined the somata as ROIs one-by-one and traced the mean intensity of pixels in each ROI over time as a calcium signal.

Photobleaching Correction
To calculate ΔF∕F, we estimated the baseline, F, from each fluorescence calcium signal by double curve fitting.First, the signal was fitted with the exponential decay function, and we built the error histogram between the fitted curve and the raw signal (first fitting).By fitting the Gaussian distribution to this histogram, we defined the range of noise values as the full-width half-maximum (FWHM) of the histogram.Then, only considering these noise values, we fitted the signal once again with a new exponential decay curve to find the baseline signal.Using this signal, we finally calculated ΔF∕F to obtain the cellular signals with zero baseline.

Calcium Spike Detection
To detect calcium spikes from each signal, we used the deconvolution method suggested by Yaksi and Friedrich 35 with slightly modified steps.First, we used a low-pass Butterworth filter (second order) to attenuate noisy fluctuation in the calcium signal (cutoff frequency: 2 Hz).Subsequently, the large peaks in this filtered signal were simply detected by finding the local maxima above the threshold (3 to 5 times the standard deviation), and the signal was segmented with the 5-s window from the peaks.All the segments were aligned to zero and averaged to create a representative form of calcium spikes for this signal.Then, we used this averaged spike to estimate the decay time [τ in Fig. 4(c), third graph] of the kernel, y ¼ expð−t∕τÞ.The signal was deconvolved using inverse filtering of this kernel, and the peaks of this trace above the threshold were detected to extract the timestamps.

Synchronous Activity Detection
We simply detected the timing of network-wide synchrony from the calcium spike trains of all ROIs in one FOV.The first step of this procedure was to obtain the network burst profile by calculating the ratio of activated ROIs to the total number of ROIs at each frame.The peaks of this profile above the threshold (here, we used 0.2-0.3 in most cases) were detected as "network burst points."At each network burst point, we defined the duration of this burst as between two frames at which all cells were silent just before and after the burst point.

Statistical Analysis
Using the temporal ROI sequences of a network burst, we implemented the nonparametric statistical analysis based on the Kendall's rank correlation for evaluating the association and consistency of the temporal orders. 8The coefficient, τ b , indicates the association between two ranking orders as follows: where n c and n d are the number of concordant and disconcordant pairs, respectively; and n is the number of the rank.As the firing order of several neurons is usually tied due to the limited temporal resolution, the additional terms, t i and s i , were complemented with the number of concurrently fired ROIs in i'th order.In this measure, τ b can span from 0 to 1; if the order of two sequences is exactly the same, τ b is 1.The pairwise comparison of all pairs of ROIs could be then clustered to identify the patterns of inter-burst propagation in the similar way to the correlation analysis.To create the matrix of τ b , we defined the distance, dτ, between two temporal orders in i'th and j'th bursts, as follows: where τði; jÞ is τ b of i'th and j'th bursts.

Procedure for Data Processing and Analysis
Our approach was aimed at converting the two-dimensional images of optical neural data to the collection of neuronal calcium spike trains ("data processing") and using these spike data for the various analyses of neural dynamics ("data analysis").To achieve the first goal of data processing, we developed a step-by-step procedure ("Basic image processing" in Fig. 1); the first step was to identify the regions of individual cell bodies ("Cell body detection" in Fig. 1).Then, the mean fluorescence intensity of each ROI over time (frame) was traced to extract cellular signals ("Calcium signal tracing" in Fig. 1).These signals were subsequently converted into the relative change of intensity after automatic baseline correction ("Background correction" in Fig. 1), and the transients of calcium signals were detected as "calcium spikes" whose peaks imply the timing of cellular firing ("Calcium spike detection" in Fig. 1).
The neuronal spike trains derived from the imaging data were then used to analyze neural activities with a multitude of quantitative measures ("Post analysis" in Fig. 1).The temporal patterns of neural activity were quantified by using mean firing rate (MFR), inter-spike interval (ISI), or the amplitude of calcium spikes ("Basic spike train analysis" in Fig. 1).The functional connectivity between neurons was estimated by calculating pair-wise cross-correlation ("Cross-correlation" in Fig. 1).Finally, we utilized neuronal spike trains to reveal the spatiotemporal firing patterns in the synchronous activity of a network.("Network burst analysis" in Fig. 1).

Sorting of Cellular Signals
To automatically identify each cell body from an image where hundreds or thousands of cells were captured, we used the elliptical morphology of cell bodies.Following the preprocessing of images to enhance the contrast between cellular regions and background ["Adjusted contrast" in Fig. 2(a)], we applied CHT 36 to the calcium imaging data in order to isolate circular elements that corresponded to the cell bodies ["Detected somata" in Fig. 2(a); for further details, see Methods section].Despite the coexistence of thin linear structures (neurites) with similar brightness, CHT could detect circular regions with their location and size, thereby isolating each cell body as an ROI [ROI; outlined with red circles in Fig. 2(a)].Subsequently, the mean fluorescence intensity of each ROI was traced in all the frames of the image set as the signal of each cell [Fig.2(b)].
Using CHT to the cell body identification requires the predetermination of two parameters: sensitivity and the range of radius.The sensitivity of CHT determines how much rounded objects would be detected; the larger this value is, the more circles CHT can detect including distorted ones.The second parameter, the range of radius, constrains the size of detectable objects.To optimize these parameters, we compared the number of cells detected by CHT with manual annotation.We used 10 different images of dye-loaded neural networks in dissociated cultures and manually counted the total number of cells in each image (∼300-500 cells).The true positive rate (TPR), indicating the ratio of automatically detected cells to the manually identified ones, gradually escalated as the sensitivity increased.The maximum TPR was achieved when the sensitivity was in the range of 0.90-0.98[Fig.2(c)].On the other hand, the false positive rate, calculated by the ratio of false alarms at each sensitivity to the maximum likelihood (when the sensitivity was 1), did not change at the sensitivity of less than 0.9 but drastically increased above this value [Fig.2(c)].Thus, the optimal sensitivity was set to be 0.9 for the accurate detection of cell bodies [Fig.2(c), indicated by a black dashed line].
With the optimal sensitivity of 0.9, we calculated TPR as the measure of detection accuracy in varying ranges of radius.All paired combinations from 3 to 31 pixels, which corresponded to about to 2-20 μm in diameter, were examined.Our examination showed that the best result of our algorithm was achieved at the range of 5-6 pixels for the minimum radius and 15-31 pixels for the maximum radius with the detection accuracy of >90% [Fig.2(d)].For the analysis of real data, we chose 5-15 pixels as the optimal range of the radius.This pixel-range corresponds to 3-10 μm in diameter for our microscopic setup (objective: 20x, image resolution: 1280 × 1080) such that the optimal range for other imaging data could be converted, according to the spatial resolution of the microscopic setup.
Next, using the optimized parameters, we assessed how many cells our algorithm can detect and how well it can separate a cellular cluster into single cells.For the quantitative evaluation, we measured the detection accuracy of our method and compared it with that of a simple thresholding method.The simple thresholding method isolated the connected pixels brighter than their surroundings as one ROI.As a result, CHT showed superior performance of individual cell separation than the thresholding method [Fig.3(a)].Compared to the manual detection, CHT could identify more cells (90.6 AE 1.7%) than simple thresholding [41.1 AE 4.8%; Fig. 3(b)].In addition, CHT could select individual cells as different ROIs, but the simple thresholding method could not separate them from the bright region in which cells were too close to each other or formed clusters [Figs.3(a) and 3(c)].More than 95% of the ROIs contained only one cell by applying CHT [Fig.3(c), right], but only half of ROIs did in case of the simple thresholding [Fig.3(c), left].Furthermore, the distribution of ROI diameters detected by using CHT also showed that it was comparable to the size of cell bodies [about 10 μm in diameter; 37 Fig. 3(d)].In addition, the minimum distance between two ROIs was 3.2 μm, implying the possibility of our method to distinguish two overlapping cells as different ROIs.As a consequence, all the results indicated that our algorithm based on CHT had the outstanding ability to not only identify most of the cells but also isolate each of them in calcium imaging data.

Calcium Spike Detection from Fluorescence Signals
After extracting the mean fluorescence intensity of each ROI from the image sequence, we applied two operations to the fluorescence signal in order to reconstruct the spike train of each neuron [Fig.4(a)].First, the raw signal was converted into the relative value of the baseline (ΔF∕F, where F is the baseline level).Second, the rapid transients of calcium signals were detected to infer the timing of the action potential firing.
In the first step, precisely estimating the trend of the baseline is essential to correct each cellular signal, as the background fluorescence intensity was not consistent but often decreased over time mostly due to photobleaching.To correct this inconsistent decay, we developed a method based on curve fitting for estimating the baseline from each cellular signal [Fig.4(b)].First, we fit an exponential decay function to the original fluorescence signal [Fig.4(b), top graph].The histogram of the errors between the original and fitted values was then fit to Gaussian distribution to determine the noise level of the signal; we decided the noise range as the FWHM of this distribution.By only using the noise values for the second fitting of a new exponential decay function, we were able to estimate the background trend (F) for each ROI.The calcium signal of each cell (ΔF∕F) was finally calculated as follows: where  subsequent exponential decay, 38 we designed a new simple method based on deconvolution, 35  Our calculation based on curve fitting successfully corrected the decreasing trend of background without the distortion of fast calcium spikes [Fig.4(d)].Furthermore, we validated the ability of our algorithm to detect calcium spikes by measuring the detection accuracy with respect to the manual annotation according to the spike amplitude [Fig.4(e)].Our method was capable of detecting more than 90% of the spikes that have the amplitude of ≥5%.These results demonstrated that our method based on the combination of noise-estimated curve fitting and deconvolution with a data-driven kernel was efficient to detect calcium spikes, particularly the small ones that may represent a single action potential. 8,29,39

Reconstruction of Neuronal Spike Trains from a Large-Scale Network
To demonstrate the applicability of our tool, we used the calcium imaging data of spontaneous neural activity measured from cultured neural networks (Fig. 5).In our experimental setup, we captured a large area (∼800 × 700 μm 2 ) of the network in which 500-1000 cells were simultaneously recorded [Fig.5 Although the neurons were cultured in the serum-free condition, glial cells often emerged in the mature networks (>2 weeks in vitro).Some of the ROIs showed slow transients of astrocytes, and others contained the mixed features of neuronal spikes and glia transients in some cases that two different cells were overlapped.Among 5989 cells in nine networks, neuronal, glia, and mixed signals were 62.0%, 6.21%, and 20.7%, respectively.The rest of the ROIs were silent (11.1%).The results indicated that our algorithm could identify not only neuronal cell bodies, but also glia cells with the mean discrimination accuracy of 79.3%, which was comparable to the previous methods. 24,33n addition, our algorithm was also able to identify silent cells, which were vital to represent neural activity in behaviors 40 and cognitive functions. 41

Post Analysis 1: Quantification of Neuronal Activities
Using our method, we quantitatively analyzed the responses of cultured neurons to the controlled extracellular environment by pharmacological treatment (Fig. 6).We used four different chemicals to control neural activity; NMDA and AMPA are the agonists of excitatory glutamatergic receptors, and GABA and bicuculline are the agonist and antagonist of inhibitory GABAergic receptors, respectively.As we reconstructed spike trains from the image datasets, we could directly calculate the MFR, ISI, the mean amplitude of calcium spikes, and the ratio of active cells as the measure of cellular activity before and after treatment.When the neural network was stimulated by additional NMDA, MFR, and amplitude increased, but ISI did not significantly change [Fig.6(a)].The results implied that NMDA increased the activity of neural networks as we expected.On the other hand, AMPA showed the opposite results; the level of all measures significantly dropped [Fig.6(b)].We speculated that the opposite responses of neural activity to NMDA and AMPA resulted from the dose-dependency.According to the previous work, 42 the firing rate of neural networks significantly increased at the low concentration of AMPA less than 1 nM, whereas it decreased at the high concentration that ranged from 5 to 100 μM.NMDA treatment also showed dose-dependency, but the range was different; the increment of the firing rate appeared at the concentration of a few μM, whereas its decrement emerged when the concentration was higher than 100 μM.The concentration of AMPA and NMDA used in this work was 20 and 10 μM, which both corresponded to the range of decreasing activity and that of increasing activity, respectively.GABA stimulation also decreased all measures of neural activity [Fig.6  (c)].Interestingly, the mean ratio of active cells after GABA treatment was 12.7 AE 2.3%, which was the similar range of inhibitory cells in cultured neural networks. 43Bicuculline treatment caused the decrement of MFR and ISI but significantly increased the amplitude of spikes [Fig.6(d)].These results were also observed in other studies 8,44 that described this phenomenon as interictal discharges or bursts.
The spike train data of neurons sorted by our method were also useful to analyze the synchronous activity of neural networks.We calculated the pairwise cross-correlation between two cellular trains and applied the hierarchical clustering method to show the synchrony level with respect to the chemical perturbation.The presence of more reddish colors in the correlation matrix after the bicuculline treatment than that of before ("base") implied that the more cells simultaneously fired [Fig.7(a)].Quantitatively, the mean correlation coefficient after bicuculline treatment was significantly higher than that of the before ("base"), also supporting that bicuculline-induced network synchronization [Fig.7(b)].Altogether, our results demonstrated that our spike train-based approach to the calcium imaging data analysis allowed quantification of individual neural activities and network synchrony for stimulus-response experiments with five distinct measures.

Post Analysis 2: Spatiotemporal Mapping of Neuronal Firing Patterns in Synchronous Network Bursts
To demonstrate the applicability of our toolbox to synchronous events in neural networks, we attempted to analyze the firing patterns of neurons in bursting networks.To detect each network burst from the spike train data [Fig.8(a), (i)], we counted the number of firing cells at each frame of an image sequence [Fig.8(a), (ii)] and divided it into the total number of cells to construct the ratio profile of firing cells over time [Fig.8(a), (iii)].By thresholding this profile, we detected peaks above the threshold as the burst point.Furthermore, the start and end points of each burst were defined as the time point at which all neurons were silent before and after the burst point.Using the real data from cultured networks, we detected the timing of network bursts and extracted all the spikes involved with , respectively, and comparable to the results using the cultured networks at the similar ages with multichannel electrophysiological tools. 45,46he location of neurons and their spikes participating in one network burst were then utilized to map the propagation of synchronous activity.As exemplified in Fig. 9, we could categorize the spikes in a single burst according to their timestamp.Figure 9(a) shows the difference in the onset timing of seven groups, indicating the capability of our spike detection method to resolve neuronal firing within a single burst.It should be noted that our data were collected with a frame rate of 30 Hz, which is much slower than the previous work for similar analysis. 8This shows that our frame rate was sufficient to segregate neurons into several groups based on their firing, and signal propagation map of a single burst could be visualized ("Pseudocolored map" in Fig. 9).
By collecting all the firing patterns of each network burst, we quantitatively analyzed the similarity of the sequences.As the patterns of signal propagation were described as a sequence of neuronal groups in our analysis, we used Kendal's τ b to adjust the tied ranks of cells in the same group.The pairwise comparison using Kendal's statistics and hierarchical clustering of the correlation matrix revealed the similar sequences of neuronal firing emerged from different network bursts [Fig.10(a)].
Prior to the clustering of the matrix, there was no apparent region that shared similar activity patterns, indicating that several patterns appeared rather than locally repeated [Fig.10(b)].To quantitatively evaluate the nonrandomness of firing patterns, we separately shuffled the order of ROIs in each burst and compared the renewed correlation matrix to the original one.As a result, the mean coefficient (τ b ) of the real data was significantly higher than that of the shuffled one, which was close to 0. The results implied that despite the entire regeneration of cellular connections in dissociated cultures, several consistent patterns of neuronal firing alternatively repeated in spontaneous, synchronous events.

Discussion
Here, we introduced NeuroCa, a toolbox that includes a series of algorithms for the automated analysis of calcium imaging data from large-scale neural networks (Fig. 11 shows the graphical user-interfaces modularized for data processing and post analysis).Our cell identification algorithm was efficient to detect cell bodies regardless of cell types or activity levels and was notably useful in separating individual ones that are very close to each other.We also devised computational strategies for adjusting the baseline level of fluorescence signals and estimating the spike waveform from the real data,   which facilitated the successful detection of small transients.Using this method, we were able to construct individual spike trains of ∼1000 cells in a submillimeter-sized network.
Our method combined the advantage of optical recording and spike train analysis, allowing us to not only easily extract the quantitative measures of neural activity but also to perform the spatiotemporal mapping of cellular firing in a population activity, especially synchronous bursts.Moreover, identifying cell bodies using a morphological filter in NeuroCa allows not only detecting most of the cells in a FOV with high accuracy, but also distinguishing individual ones from clustered regions.Essentially, CHT determines a circular object based on the intensity of local pixels, therefore being less sensitive to the imbalance background illumination than other image segmentation algorithms. 22,25,27,34In addition, much higher spatial resolution was achieved in our method; for example, our CHT-based method showed less than 5% of the error rate (the percentage of ROIs that contained multiple cells), which was two times less than the previous algorithm based on watershed transform (11%). 42The minimum cell-to-cell distance that our algorithm could separate was about 3 μm (centerto-center distance).Considered the diameter of cells as 10 μm, this feature implies the possibility of our method for isolating cells even if they are partially overlapped each other.Consequently, our method will be also applicable to the analysis of larger and denser neuronal population including tissue slices or a three-dimensional brain, where cells were far closer than dissociated cultures.
Adapting the deconvolution method with an additional step for baseline adjustment and template estimation from the real data, we were able to detect calcium spikes as small as one reflecting a single action potential firing with high accuracy.As some calcium imaging and modeling studies suggested the exponential form of calcium spikes in response to the action potentials, 38,47,48 it was reasonable to infer the timing of action potential firing by deconvolving the original fluorescent signal with an exponential kernel. 17,24,30,35,49The performance of this method critically relied on how flat the baseline of the signal was and how accurate the estimated exponential model was for each neuron.However, most practical cases include the photobleaching decay of fluorescent indicators in calcium signals, and the waveform of calcium spikes varied in different cells.Here, we exploited curve fitting in order to estimate a background trend and spike waveform for each cellular signals, thereby compensating for cell variance.This capability of our method will facilitate data analysis obtained from various neural networks regardless of cellular types or their temporal firing patterns.
Our analysis of calcium imaging data recorded from a large cultured network demonstrates the utility of NeuroCa to scrutinize the spatial and temporal patterns of neuronal activity with various quantitative measures.In the chemical perturbation study, the alteration of neural activity and synchrony was faithfully represented by five measures including MFR, ISI, calcium spike amplitude, the percentage of active cells, and pairwise cross-correlation. 11,13,33Furthermore, with the advantage of optical recording that provides the spatial information of individual cells, it was also possible to map the rapid propagation of neural activity during the short duration of synchronous events.Several previous methods were suggested to identify the coherent activity from calcium signals and provide the spatial organization of cellular assemblies, 34,50 but it is still difficult to uncover the temporal sequences of cellular firing within such short period events.Our results using NeuroCa, on the other hand, showed the spatial location of cells as well as their temporal orders, allowing us to obtain insights that are pertinent to hidden signaling pathways mediating the synchronous activity.
Our studies demonstrate that cultured neural networks revealed several consistent patterns of neuronal activity in synchronous bursting events despite the overall regeneration of cellular connectivity.In slice cultures that are composed of intact neural circuits of neurons, the sequence of cellular activation was reliably repeated in synchronous activity, 8 and several types of patterns emerged from one network over time. 1 Such repetition and diversity of consistent firing patterns have implied their role in serving as the substrate for information processing and storage. 1 However, in dissociated cultures, although such a tendency of repeated patterns across the subregions of a large network has been reported in several studies using microelectrode arrays, [51][52][53] the contribution of each cell to the activity patterns of synchrony remains unknown.Our results support these previous results and also extend the existence of sequential activation at the cellular level.
NeuroCa is essentially applicable to other types of data using different indicators (genetic 5,14 or voltage-sensitive dyes 54 ) or imaging tools (two-photon 55 or light sheet microscopic imaging 11,13,17 ).In addition, data from tissue slices or living animals were also available to be analyzed.However, for further applications, several additional algorithms will be necessary.First, NeuroCa has a limitation in detecting other neuronal areas such as dendrites 24,33,56 or neurophils, 57 which have been also of interest for investigating cellular signaling or intercellular interaction in localized regions.We anticipate that it could be resolved by combining additional morphological filters.Second, a sophisticated algorithm for separating cell types, such as neurons and glial cells, will be necessary.Third, movement correction and 3D-optimized cell identification will be required, especially for the data from living animals. 24euroCa will provide the opportunities to come close to the answer to some interesting questions in neuroscience: how individual neurons functionally participate in the various modes of brain activity and what is the role of functional interactions between neurons and glial cells in circuit dynamics?One possibility for such studies is to apply our method to in vivo systems.Recently, optical recording technologies have been exploited to the various types of living animals, such as C. elegans, 13 zebrafish, 11,12 rodents, 16,19,48,58 or primates. 59These studies have reported the impressive results of neural activity at population level and also suggested the great promise of optical neural recording technologies in functional mapping of the brain.In this stream of neuroscience research, a convenient and user-friendly tool, such as our NeuroCa, is certainly necessary for the analysis of large-scale neural data.We believe that our tool will contribute to innovative studies such as the BRAIN initiative, 60 considerably reducing the enormous loads of postexperimentation steps.NeuroCa is available for free public download (see Ref. 61).

Conclusions
We present a new computational approach for the automated analysis of calcium imaging data from a large neural network.Our methods facilitated the identification of individual cells from an image sequence and the detection of spikes from neuronal calcium signals.The calcium spike trains were useful to quantitatively analyze the overall activity of neuronal population.Furthermore, our method enabled us to map the spatiotemporal patterns of cellular firing in a large network when they synchronously fired.We believe that our approach will be a powerful tool that can contribute to the functional mapping of the brain.

Fig. 1
Fig. 1 Analytical procedure of calcium imaging data.(a,b) The image sequence is decomposed into the regions of individual cell bodies by our cell body detection algorithm, and (c) the fluorescence signal of each region of interest (ROI) was calculated by averaging the mean intensity of each region at each frame.From the calcium signals, (d) the baseline fluctuation was compensated by our curve fitting method and (e) calcium spikes were detected.(f, g) The constructed spike train data could be further used for the analysis of neural activity (f) and functional connectivity (g).(h) We have focused on investigating the neuronal firing patterns in synchronous network bursts that are spontaneously emerged from the dissociated neural networks.
F i was the mean fluorescence intensity of each ROI at the i'th frame [Fig.4(b), bottom graph].The next step was to detect spikes from the corrected signals.Considering the fluorescence signal as the superposition of calcium spikes that follow the instantaneous increase and

Fig. 2
Fig. 2 Cell body detection and parameter optimization.(a) The procedure of the cell body detection from the image sequences.(a) The image ("Original image") was preprocessed by tophat filtering and contrast level adjustment ("Adjusted contrast") and used to detect cells by applying circular Hough transform (CHT) ("Detected somata").(b) Examples of the raw fluorescence calcium signal traced from each identified cell in (a).These signals were traced from seven neurons selectively marked with the same color in (a).(c) True positive rate [(TPR); orange] and false positive rate; blue versus sensitivity (the parameter for determining the threshold of circularity) was plotted (mean AE SEM; n ¼ 10).The dashed line indicates the sensitivity of 0.9.(d) The color-coded image of detection efficiency depending on the combination of minimum and maximum radii.The scale bar in "Original image" and "Detected somata" of (a) indicates 100 μm and 30 μm, respectively.
which used a kernel of calcium spikes derived from the real data [Fig.4(c)].The signal was smoothed using a low-pass filter to attenuate the fast, noisy component [second Butterworth; Fig. 4(c), top graph].Then the large peaks were detected from the signal with thresholding [Fig.4(c), the second graph].By averaging these large spikes, we obtained a template of a spike for each calcium signal and fitted the exponential decay curve to construct a kernel [Fig.4(c), the third graph].With this data-driven kernel, the signal was subsequently deconvolved [Fig.4(c), the fourth graph].Finally, we detected the peaks of the deconvolved signal above the threshold that corresponded to calcium spikes [Fig.4(c), bottom graph; for further details, see Methods section].
(a)].Using our method, we were able to extract the activity of each cell in the form of spike trains [exemplified in Figs.5(b) and 5(b')].

Fig. 3 Fig. 4
Fig. 3 CHT identifies most of the cell bodies as a single ROI.(a) The same image processed with simple thresholding and CHT algorithm.Each area outlined a red curve in simple threshold or circle in CHT indicates each ROI.(b) The ratio of identified ROI numbers to the number of cells that were manually detected (mean AE SEM; n ¼ 10; Ã Ã Ãp < 0.001).(c) The cumulative distribution of the number of cells detected as one ROI (left) with simple thresholding (gray) and CHT (magenta) and the normalized histogram of CHT (right; mean AE SEM; n ¼ 10).(d) The distribution of ROI diameters in 10 field-of-views (gray) and average (black).The scale bar in (a) indicates 10 μm.

Fig. 5
Fig. 5 Our automated procedure extracts cellular signals and spikes from large-scale calcium imaging data (Video 1, MOV, 4.6 MB) [URL: http://dx.doi.org/10.1117/1.NPh.2.3.035003.1].(a) The first image of an image sequence recorded by calcium imaging.Colored circles outline each cell body detected by our method.(b) Each signal traced from 100 cells randomly selected from (a). (b') The inferred spikes (marked as ticks under the each signal) of four selected cells.The scale bar in (a) indicates 100 μm.

Fig. 6
Fig. 6 Quantitative analyses of the effect of chemical stimulation on cultured neural networks.Four different chemicals: (a) NMDA, (b) AMPA, (c) GABA, and (d) bicuculline was examined.Quantitative features of neural activity, such as mean firing rate (MFR), inter-spike interval (ISI), spike amplitude (amplitude), and the ratio of active cells that showed calcium signals among all detected ROIs, were measured (mean AE SEM; n ¼ 5; Ã Ã Ãp < 0.001; ns: not significant).

Fig. 7
Fig.7The effect of the bicuculline treatment on network synchronization.(a) In the case of the bicuculline treatment, the cross-correlation analysis between before (base) and after treatment (+BIC) was performed to investigate the network synchronization effect.(b) Mean correlation coefficients (MCC) of each case were used as a quantitative measure (mean AE SEM; n ¼ 5 FOVs; Ã Ã Ãp < 0.001; ns: not significant).

Fig. 8 Fig. 9
Fig.8Synchronous network burst detection from calcium spike trains.(a) The procedure for network burst detection.From the raster plot of the entire ROIs (i), we calculated the number of activate ROIs at each frame (ii) and obtained a network burst profile by normalizing each value to the total number of ROI (iii).The peaks of the profiles above the threshold (blue dashed line) were detected as the point of network bursts (indicated by blue arrowheads).The duration of each network burst was defined as the period between the frames when all ROIs were silent (indicated by red bold lines).(b) The representative examples of the real raster plot (top) obtained from Fig.5and the network burst profile (bottom).The points of network bursts and their durations were marked by blue arrowheads and red bold lines, respectively.(b') The magnified version of the ratio graph in the bottom of (b), showing the start point at which all neurons were silent before the burst point.

Fig. 10
Fig. 10 Consistent sequence of neuronal firing repeated appeared in network bursts.(a) Representative correlation matrix based on Kendall's τ b in five different networks.All matrices were sorted according to a hierarchical clustering.(b) Pairwise calculation of Kendall's τ b between all firing orders of each burst (top, left) was clustered to find consistent sequences (top, right).To statistically confirm the nonrandomness of firing patterns, the firing order of each burst was shuffled (bottom, left) and the correlation matrix was clustered as same as the original one (bottom, right).(c) When the coefficients of all possible pairs were averaged, the real data revealed much higher values than shuffled [mean AE SEM; n ¼ 15 ("Real") and 75 ("Shuffled"); Ã Ã Ãp < 0.001 under paired t test].