2.1.

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\mu \mathrm{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 $\sim 1000\mathrm{cells}/{\mathrm{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.

2.2.

Calcium Imaging

We used a commercial calcium indicator, Oregon green BAPTA-1 (OGB-1; $50\mu \mathrm{g}$; O6807, Molecular Probes®, California), for measuring neural activities. The stock solution of the calcium indicator was made of OGB-1 (50 mg), $48\mu \mathrm{l}$ of dimethyl sulfoxide (DMSO; D4540, Sigma, Missouri), and $2\mu \mathrm{l}$ of F127 solution [Pluronic® F-127; 20% (w/v) in DMSO; P2443, Sigma, Missouri]. We also used buffered artificial cerebrospinal fluid (bACSF) composed of ${\mathrm{NaHCO}}_3$ (25 mM), D-glucose (25 mM), NaCl (125 mM), KCl (2.5 mM), ${\mathrm{NaH}}_2{\mathrm{PO}}_4$ (1.25 mM), ${\mathrm{MgCl}}_2·6{\mathrm{H}}_2\mathrm{O}$ (1 mM), and ${\mathrm{CaCl}}_2·2{\mathrm{H}}_2\mathrm{O}$ (2 mM) in third distilled water and supplemented with 10 mM of HEPES. Prior to the preparation of imaging, all materials were heated in an incubator (37°C, ${\mathrm{CO}}_2$ 5%).

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\mu \mathrm{m}\times 700\mu \mathrm{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\mathrm{GB}/\mathrm{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).

2.3.

Chemical Stimulation

Several agonists of neurotransmitter receptors, such as N-methyl-D-aspartic acid (NMDA; M3262, Sigma, Missouri), $\alpha$-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA; Asc-130, Ascent Scientific), and $\gamma$-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\mu \mathrm{M}$ of concentration and AMPA was dissolved into bACSF with $40\mu \mathrm{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\mu \mathrm{M}$, and that of AMPA was $20\mu \mathrm{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.

2.4.

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.

2.5.

Photobleaching Correction

To calculate $\mathrm{\Delta }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 $\mathrm{\Delta }F/F$ to obtain the cellular signals with zero baseline.

2.6.

Calcium Spike Detection

To detect calcium spikes from each signal, we used the deconvolution method suggested by Yaksi and Friedrich³⁵ 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 [$\tau$ in Fig. 4(c), third graph] of the kernel, $y=\exp (-t/\tau )$. 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.

2.7.

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.

2.8.

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.⁸ The coefficient, ${\tau }_b$, indicates the association between two ranking orders as follows:

3.1.

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).

3.2.

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³⁶ 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)].

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 ($\text{mean}\pm \mathrm{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\mu \mathrm{m}$ and $30\mu \mathrm{m}$, respectively.

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 ($\sim 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\mu \mathrm{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\mu \mathrm{m}$ in diameter for our microscopic setup (objective: 20x, image resolution: $1280\times 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\pm 1.7\%$) than simple thresholding [$41.1\pm 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\mu \mathrm{m}$ in diameter;³⁷ Fig. 3(d)]. In addition, the minimum distance between two ROIs was $3.2\mu \mathrm{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.

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 ($\text{mean}\pm \mathrm{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; $\text{mean}\pm \mathrm{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\mu \mathrm{m}$.

3.3.

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 ($\mathrm{\Delta }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.

Fig. 4

Spike detection algorithm including baseline estimation and deconvolution with data-driven kernels. (A) To detect the calcium spikes from the fluorescence signal of each cell, we used a two-step procedure: (b) the correction of baseline level and (c) deconvolution with data-driven kernels. (B) In the first step for removing the baseline trend, the raw signal was fitted to the exponential function (red dashed line) at first (top graph). Then, the noise level was estimated by means of Gaussian distribution fitting (black dashed line) to the error histogram (second graph). The values within the full-width half-maximum (FWHM) were considered as noise. Using only noise values (red dots), the second fitting was carried out to find more precise baseline (third graph). This second fitted signal was finally concerned as the baseline (blue dashed line indicated as “F”). The relative ratio of fluorescence intensity change to the baseline was calculated (bottom graph). (C) From this corrected signal, we could detect calcium spikes. To reduce the noise, we smoothed the signal using low-pass filter (top graph), and detect only large peaks (marked with red circles) by simply thresholding (red dashed line) the filtered signal (second graph). Using a 5 s window (blue dashed rectangles), we collected the large spikes and found the parameter, $\tau$, to determine the data-driven kernel (third graph). Then we deconvolved the filtered signal using this kernel and detected the peaks above the threshold (red dashed line; fourth graph; red arrowheads: hits, black arrowhead: miss). These peaks were then concerned as the timing of spikes (bottom graph; red ticks: detected spikes, gray tick: missed spike). (D) Validation of the ability of our method. Raw neuronal calcium signals (I) were corrected by our curve fitting method (ii). (E) The detection accuracy of our method from corrected signals with respect to the manual inspection was measured. Each point represents the accuracy of spike detection in the amplitude range of $i\%$ to less than $(i+1)\%$.

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 ($\mathrm{\Delta }F/F$) was finally calculated as follows:

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 subsequent exponential decay,³⁸ we designed a new simple method based on deconvolution,³⁵ 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].

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 $\ge 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

3.4.

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 ($\sim 800\times 700\mu {\mathrm{m}}^2$) of the network in which 500–1000 cells were simultaneously recorded [Fig. 5(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. 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\mu \mathrm{m}$.

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,33 In addition, our algorithm was also able to identify silent cells, which were vital to represent neural activity in behaviors⁴⁰ and cognitive functions.⁴¹

3.5.

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.

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 ($\text{mean}\pm \mathrm{SEM}$; $n=5$; $***p<0.001$; ns: not significant).

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,⁴² 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\mu \mathrm{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 $\mu \mathrm{M}$, whereas its decrement emerged when the concentration was higher than $100\mu \mathrm{M}$. The concentration of AMPA and NMDA used in this work was 20 and $10\mu \mathrm{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\pm 2.3\%$, which was the similar range of inhibitory cells in cultured neural networks.⁴³ Bicuculline 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.

Fig. 7

The 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 ($\text{mean}\pm \mathrm{SEM}$; $n=5$ FOVs; $***p<0.001$; ns: not significant).

3.6.

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 each burst [Figs. 8(b) and 8(b’)]. Under our definition, we observed that all neurons were silent $133\pm 7.21\mathrm{ms}$ before and $230\pm 13.3\mathrm{ms}$ after the burst point, respectively ($\text{mean}\pm \mathrm{SEM}$, $n=161$ bursts in 11 networks). The mean duration and inter-burst interval of network bursts was $364\pm 15.4\mathrm{ms}$ ($\text{mean}\pm \mathrm{SEM}$, $n=161$ bursts in 11 networks) and $24.5\pm 2.04\mathrm{s}$ ($\text{mean}\pm \mathrm{SEM}$, $n=150$ bursts in 11 networks), respectively, and comparable to the results using the cultured networks at the similar ages with multichannel electrophysiological tools.^45,46

Fig. 8

Synchronous 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. 5 and 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.

The 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.⁸ This 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).

Fig. 9

Spatiotemporal mapping of signal propagation within a network burst. (a) The kinetics of calcium spikes in each group (categorized according to the timing of detected spikes). The spikes of each cell were marked as gray curves, and the average spike in each group was marked as a black curve. Two red dashed lines indicate the onset timing of calcium spikes in the first and last group. (b) Black dots indicate the neurons fired at each time point and gray dots silent ones. These time-series images were then merged into one by color-coding (“Pseudocolored map”); in this image, early firing neurons were marked as red colors and late firing neurons as blue colors.

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 ${\tau }_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 alternatively 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 (${\tau }_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.

Fig. 10

Consistent sequence of neuronal firing repeated appeared in network bursts. (a) Representative correlation matrix based on Kendall’s ${\tau }_b$ in five different networks. All matrices were sorted according to a hierarchical clustering. (b) Pairwise calculation of Kendall’s ${\tau }_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 [$\text{mean}\pm \mathrm{SEM}$; $n=15$ (“Real”) and 75 (“Shuffled”); $***p<0.001$ under paired $t$ test].