Spectral Entropy Based Neuronal Network Synchronization Analysis Based on Microelectrode Array Measurements

Synchrony and asynchrony are essential aspects of the functioning of interconnected neuronal cells and networks. New information on neuronal synchronization can be expected to aid in understanding these systems. Synchronization provides insight in the functional connectivity and the spatial distribution of the information processing in the networks. Synchronization is generally studied with time domain analysis of neuronal events, or using direct frequency spectrum analysis, e.g., in specific frequency bands. However, these methods have their pitfalls. Thus, we have previously proposed a method to analyze temporal changes in the complexity of the frequency of signals originating from different network regions. The method is based on the correlation of time varying spectral entropies (SEs). SE assesses the regularity, or complexity, of a time series by quantifying the uniformity of the frequency spectrum distribution. It has been previously employed, e.g., in electroencephalogram analysis. Here, we revisit our correlated spectral entropy method (CorSE), providing evidence of its justification, usability, and benefits. Here, CorSE is assessed with simulations and in vitro microelectrode array (MEA) data. CorSE is first demonstrated with a specifically tailored toy simulation to illustrate how it can identify synchronized populations. To provide a form of validation, the method was tested with simulated data from integrate-and-fire model based computational neuronal networks. To demonstrate the analysis of real data, CorSE was applied on in vitro MEA data measured from rat cortical cell cultures, and the results were compared with three known event based synchronization measures. Finally, we show the usability by tracking the development of networks in dissociated mouse cortical cell cultures. The results show that temporal correlations in frequency spectrum distributions reflect the network relations of neuronal populations. In the simulated data, CorSE unraveled the synchronizations. With the real in vitro MEA data, CorSE produced biologically plausible results. Since CorSE analyses continuous data, it is not affected by possibly poor spike or other event detection quality. We conclude that CorSE can reveal neuronal network synchronization based on in vitro MEA field potential measurements. CorSE is expected to be equally applicable also in the analysis of corresponding in vivo and ex vivo data analysis.

Synchrony and asynchrony are essential aspects of the functioning of interconnected neuronal cells and networks. New information on neuronal synchronization can be expected to aid in understanding these systems. Synchronization provides insight in the functional connectivity and the spatial distribution of the information processing in the networks. Synchronization is generally studied with time domain analysis of neuronal events, or using direct frequency spectrum analysis, e.g., in specific frequency bands. However, these methods have their pitfalls. Thus, we have previously proposed a method to analyze temporal changes in the complexity of the frequency of signals originating from different network regions. The method is based on the correlation of time varying spectral entropies (SEs). SE assesses the regularity, or complexity, of a time series by quantifying the uniformity of the frequency spectrum distribution. It has been previously employed, e.g., in electroencephalogram analysis. Here, we revisit our correlated spectral entropy method (CorSE), providing evidence of its justification, usability, and benefits. Here, CorSE is assessed with simulations and in vitro microelectrode array (MEA) data. CorSE is first demonstrated with a specifically tailored toy simulation to illustrate how it can identify synchronized populations. To provide a form of validation, the method was tested with simulated data from integrate-and-fire model based computational neuronal networks. To demonstrate the analysis of real data, CorSE was applied on in vitro MEA data measured from rat cortical cell cultures, and the results were compared with three known event based synchronization measures. Finally, we show the usability by tracking the development of networks in dissociated mouse cortical cell cultures. The results show that temporal correlations in frequency spectrum distributions reflect the network relations of neuronal populations. In the simulated data, CorSE unraveled the synchronizations. With the real in vitro MEA data, CorSE produced biologically plausible results. Since CorSE analyses continuous data, it is not affected by possibly poor spike or other event detection quality. We conclude that CorSE can reveal neuronal network synchronization based on in vitro MEA field potential measurements. CorSE is expected to be equally applicable also in the analysis of corresponding in vivo and ex vivo data analysis.

INTRODUCTION
Temporally correlated activity between neurons or neuronal networks in vivo and in vitro has been vastly studied in terms of event based synchrony, or synchrony between oscillations or rhythmic activities in different frequency bands. Salinas and Sejnowski (2001) argued that the presence of correlations between the activities of pairs of neurons, or synchrony per se, is not important in general, since they may arise from common inputs or synaptic interactions, or from overlapping perceptive fields, respectively; however, changes in the correlation structure of a neuronal network reflect changes in its functional connectivity. Previous studies have shown that the pattern of synchronization determines the pattern of neuronal interactions, and that the efficiency of transferred information is also modulated by synchrony (Buehlmann and Deco, 2010;Battaglia et al., 2012). Thus, assessing the relations of synchrony is essential not only for fully developed neuronal networks, such as in the brain, but also for the assessment of development and plasticity of cultured neuronal networks.
In the past years, several studies concentrated on quantifying and analyzing the network relations of cultured neuronal cells (Garofalo et al., 2009;Mack et al., 2014). Most of the studies utilized binary analysis based on events, particularly the occurrences of spikes and bursts. For example, in several studies, transfer entropy (TE), joint entropy, mutual information (MI), coincidence index, and event synchrony (ES) were employed on detected spikes to evaluate network relations (Quiroga et al., 2002;Garofalo et al., 2009;Ito et al., 2011); however, as stated by Buzsáki et al. (2012), network relations affect local field potentials (LFPs) as well.
For brain studies, synchronization, causality, phase and frequency coupling, or tracking previously defined rhythms, can reveal network interactions (Ginter et al., 2005;Buehlmann and Deco, 2010). A review of a few connectivity measures to assess neuronal activity has been presented by Bastos and Schoffelen (2016). The amount of propagating activity observed in different frequency bands (rhythmic activities) or associated with well-defined electroencephalogram (EEG) rhythms (delta, theta, alpha, beta, and gamma) is important in interpreting the results of such a study (Ginter et al., 2005). On the other hand, generally in cultured neuronal networks, the different frequency bands are not as distinguishable as in the brain studies, or the absence of well-defined rhythms makes the analysis more challenging. Even though LFPs (or raw recordings that may include both spikes and LFPs) potentially carry information on network relations, they are not commonly used for the network analysis based on microelectrode array (MEA) measurement data from cultured cells. MEAs are usually used to measure extracellular field potentials from electrically active tissues and cell cultures at network and cell levels (Thomas et al., 1972;Gross et al., 1977;Pine, 1980;Egert et al., 1998). MEA electrodes record field potentials, e.g., from the neurons in their vicinity, which can carry contributions from both extracellular action potentials (EAPs) from individual neuronal cells and lower frequency contributions originating from neuronal population activity. Neurons may temporarily arrange themselves into synchronous functional ensembles to perform a given task. These ensembles may be volatile and only exist for short periods of time before new ensembles with partially different subsets of neurons are formed. Connected neuronal ensembles are thought to operate at certain frequencies (for general references, see Buzsáki and Chrobak, 1995;Penttonen and Buzsáki, 2003;Buzsáki and Draguhn, 2004). Consequently, frequency domain analysis has potential to obtain novel information also from cell cultures (Jarvis and Mitra, 2001;Brown et al., 2004). Frequency spectrum analysis may also be a good alternative in cases with unreliable spike detection either due to low amplitude spikes in noise or conflicting results from different spike detection algorithms.
Drawing from above, we have hypothesized that also temporal correlations of the frequency spectrum distributions could reflect the network relations of neuronal populations (Kapucu et al., 2016a). Intuitively, this is motivated by the possibility that measurements from functionally connected neuronal populations may be quite different if only time domain properties were considered. For analyzing the functional connectivity of a network, techniques quantifying the spectral properties of neuronal ensemble activity provide promising alternatives to the methods assessing the couplings or correlations between specific rhythms or frequencies. Spectral entropy (SE) quantifies the regularity, or complexity, of a signal based on its frequency dynamics. SE is a frequency based realization of Shannon's entropy algorithm (Shannon, 1948), which was previously used for analyzing certain neuronal events, such as burst suppression, and for the EEG based assessment of the depth of anesthesia (Viertiö-Oja et al., 2004). In our previous work (Kapucu et al., 2015), we utilized SE and sample entropy to quantify in vivo and in vitro neuronal bursts according to their complexities, and demonstrated similarities in the complexity values of bursts from neighboring channels. Also in Kapucu et al. (2016a), we tested the feasibility of SE for the assessment of synchronization. However, the method was never validated with simulations and its true applicability was not demonstrated with larger real measurement datasets. An earlier study of the relations between the synchronization and the activity level, as well as the relations between synchronization and connectivity levels (Chawla et al., 1999), was also a motivation for the evaluation of CorSE with different levels of connected networks.
In this paper, we expand the original idea proposed by Kapucu et al. (2016a) and investigate the benefits of SE time course correlation analysis as a tool for analyzing synchronization by analyzing a larger set of data. Here, we also name the proposed method CorSE. Firstly, we illustrate CorSE with a toy simulation of neuronal ensembles (Montgomery, 2014). Next, we validate our method with simulated MEA data produced with computational integrate-and-fire model based neuronal networks with known connectivity levels. The results of CorSE are compared to and assessed together with the results from three existing event based synchronization assessment algorithms, ES (Quiroga et al., 2002), MI (Gray, 1990), and TE (Schreiber, 2000). Finally, we demonstrate the applicability of CorSE with MEA data measured from cultured rat cortical neurons with different activity levels, and with MEA data measured from a developing network of mouse cortical neurons.

Biological Data
In vitro MEA experiments were performed with dissociated rat and mouse cortical cell cultures. The data from dissociated rat cortical cells was originally collected for a previous study (see the details given by Weihberger et al., 2013); animal treatment was according to the Freiburg University (Freiburg, Germany) and German guidelines on the use of animals in research. Briefly, rat cortical cells were obtained from prefrontal cortical tissue of newborn Wistar rats and plated on MEA plates, which consist of 60 titanium nitride electrodes of 30 µm diameter and 200 µm interelectrode spacing on an 8 × 8 rectangular grid with corner electrodes missing (model: 60MEA200/30iR, Multi-Channel Systems MCS GmbH, Reutlingen, Germany). Seeded cell density was approximately 1250 cells/mm 2 . After 4 weeks of culturing, the cultures were considered mature (Wagenaar et al., 2006) and recordings were conducted inside a dry incubator.
To assess in vitro network development over time, commercially available mouse cortical cells (A15586, Gibco, Thermo Fisher) were plated on MEAs similar to those described above. Briefly, the MEAs were coated with poly-L-lysine and laminin, and the cells were sowed as droplets on to MEA plates for culturing (Wagenaar et al., 2006). Electrophysiological data were recorded three times a week starting from 4 days in vitro (DIV) until the 29th DIV. Every recording lasted for 5 min.
The recordings from rat cortical cell cultures were analyzed using the different synchronization assessment methods considered in this paper. All in vitro data was first filtered with a 50 Hz notch filter and 7 Hz high pass filter to alleviate the powerline noise and low frequency fluctuations, respectively.
For the analysis methods that are based on spike time stamps, spikes were detected using two thresholding methods to evaluate the effects of different spike detection methods on the results of the synchronization assessment algorithms: spike detection thresholds were set to five times the standard deviation of the signal, or at five times the estimated standard deviation of the background noise of the band pass filtered signal as proposed by Quiroga et al. (2004). Here, the different thresholding methods are denoted by STD and eSTD, respectively. CorSE was employed to demonstrate the tracking of network development in the mouse cortical cell cultures. Since CorSE operates on the measured signals themselves (either raw or filtered), spike detection was not performed.
Two cultures were selected for the analysis according to the number of active locations where spiking activity was observed according to a house made rule of 50 spikes per 300 s (Kapucu et al., 2012(Kapucu et al., , 2016b so that the chosen MEAs had different numbers of active locations. The analyzed cultures are denoted by MEA1 and MEA2. The MEA1 had 32 or 37 active sites as calculated with STD and eSTD, respectively, and MEA2 had 15 or 12 active sites as calculated with STD and eSTD, respectively.

Toy Simulations of Correlated Time-Variant SEs
A toy simulation of a MEA signal was generated to illustrate and explain the proposed SE based synchronization assessment method CorSE. Here, MEA measured neuronal unit activities, i.e., EAPs, and the average activity of neuronal ensembles, i.e., LFPs, were simulated as cardinal Sine, i.e., Sinc waves, and oscillation of Sines, respectively (see Montgomery, 2014 for the Sine wave representation of LFPs), to generate a simulated population signal as seen via an electrode. In such a model, neuronal synchronization between two populations was defined as simultaneous activity of neuronal ensembles in a simple way: when a number of neuronal ensembles (illustrated with green, red, and blue color filled ellipses in Figure 1A) activated in a population, the same number of neuronal ensembles was also activated in the other population, i.e., the model assumes that all the active neuronal ensembles are connected to other ensembles in the other population. Three populations were simulated with two of them (populations 1 and 2) fully synchronous with each other, and the third population (population 3) independent of the other populations ( Figure 1A).
Simulated signals were generated in 1 s sections which were simply concatenated to form a simulation of a 3 min measurement. For each 1 s section, a random number of constituent Sine signals, and a random number of Sinc signals, were generated and summed together. The sampling rate for the generated signals was 1 kHz. The number of constituent Sine signals was evenly distributed between 5 and 10, and the number of constituent Sinc signals evenly distributed between 0 and 10. The oscillation amplitudes, frequencies, and phases were randomly selected and evenly distributed. To generate the two connected populations, whenever a Sine or Sinc appeared in one of them, a Sine or Sinc, respectively, appeared also in the other population, both with independently random amplitudes, frequencies, and phases. The signals for the unconnected population were generated independently of the other two populations.
The simulations were implemented in Matlab (MathWorks, Inc., Natick, MA, USA). A large data set of 1000 triplets of populations 1, 2, and 3 was generated, each simulated signal corresponding to a 3 min recording. To investigate the functioning of the methods with either LFPs or EAPs or both, signals with five different EAP-LFPs power ratios P EAPs/LFPs = 0, 10, 20, 50, and 100% were generated (P EAPs and P LFPs denote the total estimated powers of the summed constituent Sinc and Sine signals, respectively).
Statistical validation was approached by calculating all pairwise synchronizations between all 1000 simulated recordings in all three populations, and by calculating statistics for the signals from the two connected populations to express detected synchronization vs. the synchronization between signals from the unconnected populations.

Simulations with Integrate-and-Fire Model Based Neuronal Networks
The computational neuronal network simulation was based on the model introduced by Tsodyks et al. (2000). The networks consisted of integrate-and-fire neurons with short-term plastic synapses, and exhibited clear population bursts. The parameters employed for neurons and synapses were the same as in Tsodyks et al. (2000). To introduce spontaneous activity in the network, the neurons were driven with white noise current. Three-minute MEA measurements were simulated.
The simulated networks consisted of two populations. Each population consisted of 50 neurons, of which 40 were excitatory and 10 inhibitory. Each population was internally fully connected, but without autapses. Simulations were conducted at five inter-population connectivity levels, 0, 10, 20, 50, and 100%, between the two populations. Hundred percent connected populations correspond to one population twice the size of the original populations. Here, the percentages give the probabilities for one neuron to be connected to another neuron in the other population. Hundred pairs of populations were simulated for each inter-population connectivity level. The weights of all connections were tuned so that the mean spike rate in a population would not vary highly between the simulations with different levels of connectivity, resulting in the mean spike rate of 39-45 spikes/second in a simulation.
The NEST simulator (Gewaltig and Diesmann, 2007) was used to provide spike time stamps of the EAPs of the individual neurons. From the time stamps, artificial MEA recordings were constructed to simulate raw MEA recordings: the time resolution of the simulation was 1 ms, and a single spike was recorded for a one-millisecond time bin if any number of individual spikes appeared during the bin (Figure 2A). Thereafter, a Sinc function with random parameters, similarly as in the toy simulation, was formed for each spike and located in time with the maximum at the spike time point. The generated Sinc signals, peaking in general at different points in time, were summed to generated an EAP signal. Basically, a Sinc kernel can be used to reconstruct a sampled signal (Blanche and Swindale, 2006) or a population activity (Nawrot et al., 1999), but here we intended not to reconstruct the original recording from its samples, but instead to obtain a continuous function based on the simulated spike time stamps. We call these simulated signals artificial raw recordings (c.f., Figure 2B).
Artificial LFPs were added to the EAP simulations as described in Section CorSE, SE Based Synchronization Analysis (c.f., Figure 2B). Simulations were conducted at P EAPs/LFPs ≈ 20% and P EAPs/LFPs ≈ 50%, to roughly correspond to two difference scenarios of background noise and activity levels vs. action potential amplitudes. Thereafter, to make the simulations more realistic, white Gaussian noise (WGN) was added to the generated artificial raw signals for all cases simulated. WGN was added to obtain signal to noise ratios (SNRs) of 50% and 20%, as calculated by SNR = 100 · P EAPs /P WGN , where P EAP is the estimated power of a generated artificial recording, and P WGN is the estimated power of the added WGN.

Shannon Entropy
Entropy in the context of information theory was introduced by Shannon as a measure of uncertainty (Shannon, 1948). Shannon entropy is defined as where p i is the probability that an amplitude value occurs in the ith amplitude bin, and is given by the probability density function of the time series.

SE
We calculated SE as Shannon's entropy on power spectrum as described by Viertiö-Oja et al. (2004): SE was calculated by first obtaining the frequency spectrum of the time series x(n), sampled at discrete time points n, by fast Fourier transform X f at frequency points f (2) with bold denoting a vector.
The power spectrum is given by where X * f is the complex conjugate of X f . Here, power spectrum was estimated by Welch periodogram with a Hann window of a preset length of 0.5 s and 50% window overlap to provide a smoother transition between windows and to increase temporal resolution. Power spectrum was normalized with a constant C for the Nyquist frequency range at K frequency points f 1 , . . . , f k , . . . , f K so that the sum of the normalized power spectrum P norm f k equaled unity (4).
SE S was calculated from the normalized power spectrum as and S was normalized to reside between 1 and 0 by In the sequel, SE is always considered to be the normalized SE S norm .

Correlation of SEs
Here, we quantify the synchronization of signals by the correlation of the temporal changes of their spectral contents. This is obtained by calculating SEs in time windows, and the degree of common temporal changes for different sites in the neuronal network is assessed by calculating cross correlations: the cross covariance C S x S y (t, t + τ ) of the SEs S x and S y of the signals x and y, respectively, describes how well the SE of the signal y at time t + τ is correlated with SE of the signal x at time t. Here, the sample cross correlations were estimated using the crosscorr function of Matlab (Chatfield, 2003). With O the number of time windows in which SEs were calculated, and i the index of a 0.5 s long time window, S x,i , S y,i ; i = 1, 2, . . . , O , cross covariance at the lag l = 0 is given by where S x and S y are the sample means of the corresponding SEs. The cross correlation r S x S y at lag l = 0 is then estimated as where σ S x and σ S y are the standard deviations of the corresponding SEs. SE cross correlation (8) values at lag zero of the were used to assess the level of synchronization between pairs of channels, i.e., CorSE between signals x and y is given by CorSE xy = r S x S y .

Known Event Based Methods Used for Comparison
For comparisons with CorSE, we implemented three commonly used event based synchronization assessment algorithms: ES (Quiroga et al., 2002), MI (Gray, 1990), and TE (Schreiber, 2000). Comparisons were made based on the results of all the algorithms considered for the computational neuronal network simulations and for the rat cortical cell recordings. EAPs detectable in these signals were taken as the events for the three event based algorithms considered. The spikes, i.e., EAPs, in both the artificial recordings and in real MEA recordings were detected with two different threshold based spike detection methods, STD and eSTD, as described in Section Biological Data, resulting in sets of binary strings as the inputs to algorithms. With the artificial recordings, the effects of added LFPs and WGN on the spike detection, and the consequences on the event based synchronization assessment algorithms, were also evaluated. ES, introduced by Quiroga et al. (2002), measures synchronization based on quasi-simultaneous appearance of events. As an initial step, the ES algorithm finds from the time series x and y the maximum time period τ i,j between two consecutive spikes so that the two spikes occurring at times t x i and t y j in signals x and y, respectively, where i and j are spike indexes, can be considered simultaneous, by calculating the local spike appearances: Then, cross covariance C x|y is defined as the number of times a spike appears in x within τ i,j after a spike has appeared in y, as: where M x and M y are the total numbers of events for x and y , respectively, and Finally, synchronization is given by The algorithm in its original form did not have a requirement for the minimum number of events to consider; thus, two time series with even one simultaneous detectable event could be considered as fully synchronized. To circumvent this for MEA recordings, we employed the in-house criterion of minimum 50 spikes in 300 s, as in Kapucu et al. (2012) and Kapucu et al. (2016b), for the data to be analyzed. This eliminates unjustified synchronizations, which could be caused by coincidentally appearing rare events. We named the modified algorithm the corrected ES (cES) method. The next method we used for comparisons was the well-known mutual information (MI) (Gray, 1990), which has been widely employed in many studies to quantify dependencies between time series or specifically for quantifying synchronization (Garofalo et al., 2009;Salazar et al., 2012). MI (13) where e x i and e y i are ith single events occurring in the signals x and y, respectively, p(e x i , e y i ) is the joint probability density function, and p(e x i ) and p(e y i ) are the single probabilities. Mutual information is a symmetric measure, i.e., MI xy = MI yx .
The last algorithm is transfer entropy (TE), which extends the concept of MI to conditional properties by considering the history of the influenced information (Schreiber, 2000). In other words, TE y→x measures the increase in predictability of knowing the future and the past of x, once y is known. TE can be calculated as where in addition to the parameters defined above, p e x i+1 |e x i denotes the conditional probability of observing the state e x i+1 after e x i . TE (15) was calculated in one-millisecond signal bins for delays up to three bins, and the maximum value of TE was considered as given in Ito et al. (2011).
where d is the time delay. Since TE is asymmetric and we did not consider the directionality, we considered the maximum TE value obtained between the channel pairs, i.e., TE = max(TE y→x , TE x→y ).

Assessment of the Toy Simulations
To present the results of the toy simulation, it is necessary to define the detection of correct synchronization, false positive synchronization, and false negative synchronization, and to calculate their rates. Here, synchronization is detected correctly if CorSE between populations 1 and 2 is larger than CorSE between populations 2 and 3, and larger than between populations 1 and 3, i.e., when CorSE XY > CorSE XZ and CorSE XY > CorSE YZ . In this case, populations 1 and 3, and populations 2 and 3 are not deemed synchronized. False negative synchronization is detected when the true synchronization is missed, i.e., when either CorSE between populations 1 and 3, or between populations 2 and 3, is greater than or equal to CorSE between populations: CorSE XY ≤ CorSE XZ and/or CorSE XY ≤ CorSE YZ . False positive synchronization is detected when either populations 1 and 3, or populations 2 and 3 are deemed synchronized. In all cases of false positive detection, the true synchronization is also missed. Vice versa, in all cases of false negative synchronization, a false positive synchronization is detected. Thus, due to the system of populations in Figure 1, and the synchronization criterion applied, the false positive and false negative rates are equal.

Comparative Assessment of the Analysis Methods for MEA Recording Analysis
For comparing the synchronization values obtained with the different methods considered, we first calculated the synchronization between each electrode pair, and thereafter created adjacency matrixes of the synchronization values. In Figures 5, 6, the matrixes are arranged for better visualization according to the ascending channel numbers on a MEA plate, e.g., channel 12 is located at (1,2) on an 8 × 8 MEA layout. In each matrix, the values were normalized respect to the maximum of the matrix to form a color scale. This was done to make the locations of the maximally synchronized electrode pairs more easily comparable between the different methods.
Since the signals were recorded from unguided neuronal cells that freely form networks on the MEA plates, there was no ground truth available on the actual connections or synchronization between the neuronal populations. However, we compared the results from the different algorithms based on a fixed number of most synchronized channel pairs according to each algorithm, i.e., the strongest synchronizations; to illustrate the differences in the results of difference algorithms, the 40 strongest synchronizations were found using each algorithm for MEA1, and the 20 strongest synchronizations for MEA2.
Also, we evaluated the similarity of the synchronization based functional connectivity results of the algorithms. Here, we considered the measurement channels as nodes and synchronized channel pairs as links. The results are presented as the numbers of links and nodes common between the results from the different methods, when considering only 10, 20, ..., or 50 strongest synchronizations found with each algorithm, i.e., synchronization maps were drawn with each method, and similar findings between pairs of methods were reported by plotting the numbers of the found common nodes and the numbers of common links as functions of the number of the strongest synchronizations considered.

Toy Simulation Analysis Results
The Results from the analysis of the toy simulations with 1000 triplets indicate that CorSE was able to clearly distinguish the two populations with time correlated frequency distributions from the population pairs that did not have the correlated frequency distributions by design. In Table 1, results are shown for the model with LFPs, with EAPs, and for a model with both LFPs and EAPs with different P EAPs/LFPs s (see Section SE). It is seen from Table 1 that over 96% of the cases were correctly identified as synchronized, and the results show consistent performance regardless of the EAP-LFP power ratio considered. In summary, the results in Table 1 demonstrate that the method can detect synchronization in the case simulated.

Integrate-and-Fire Model Simulation Analysis Results
We assessed the simulated integrate-and-fire model based computational neuronal networks with CorSE and compared the results with the three different event based synchronization assessment algorithms described in Section Known Event Based Methods Used for Comparison. Firstly, in Table 2, we present the spike detection results employing two different spike detection methods, i.e., STD and eSTD, and with respect to the different 99.5 0.5 X, Y, and Z, denote the recorded signals from the populations 1, 2, and 3 respectively, with only X and Y with mutually correlated frequency distributions. For this analysis system, false positive rate is equal to false negative rate. EAP-LFPs power ratios considered, as well as different with levels of added WGN. Table 2 shows the approximate mean values in percentages for false negative (FN) detections, i.e., the missed spikes for all 1000 artificial recordings created from 100 simulations for each five connectivity levels. False positive detection, i.e., detected spurious spikes, were observed very rarely and could be considered negligible in all cases. Results are given for the both spike detection methods for P EAPs/LFPs ≈ 20 and 50%, as well as for SNRs of 20 and 50%. The results in Table 2 indicate that both additive LFPs and WGN greatly affect spike detection, which is natural; the phenomenon occurs in all thresholding based spike detection systems, since with increasing LFP and/or WGN power, more and more spikes fall below overall noise level and cannot be detected by thresholding, and subsequent analysis is done based on only the detectable spikes. In practice in general, the FN detections correspond to the action potentials from neurons far away from the measurement electrode, so that the action potential amplitudes fall below the general noise level; such neurons are naturally more abundant than the neurons in the close vicinity of the electrode. Next, we evaluated the synchronization values obtained by the different algorithms for the different connectivity levels. Exemplary raster plots for individual neurons and their corresponding artificial population activity for different connectivity levels are presented on the left-hand and right-hand panels of Figure 3, respectively. By visual inspection, inter-and intra-population synchronization can be observed in the lefthand panels of Figure 3. In Figure 4, we present the calculated synchronicities from the artificial recordings in the left-hand panels, and the corresponding sample signals in the right-hand panels. Figure 4A presents the synchronicities calculated from the artificial recordings with only EAPs (see the right-hand panel for an exemplary 2 s signal segment). Synchronization values calculated by all the algorithms increase with the increasing connectivity, except for the connectivity levels 50 and 100% for all other methods than CorSE. CorSE was the only method showing monotonous increase in the synchronization when the connectivity increased for the considered connectivity values.
Figures 4B,C presents the synchronicities calculated from the artificial recordings with EAPs and added fully synchronized LFPs. Comparisons between the results with the EAP-LFP power rations P EAPs/LFPs ≈ 20% vs. 50%, and P EAPs/LFPs ≈ 20% vs. 20% can be seen in Figures 4B,C, respectively. With the added artificial LFPs (see Figures 4B,C right-hand panels for exemplary 2 s signal segments with P EAPs/LFPs ≈ 50 and 20%, respectively), the synchronization levels detected with the event based algorithms were noticeably smaller than those detected with CorSE as predictable from the FN spike detection rates shown in Table 1. With LFPs or WGN (Figures 4B-E), it is observed that the increased connectivity did not always lead to increased observed synchronization. Results indicate that the spike detection performances have a strong influence on the synchronization results of the event based algorithms, as expected. In contrary, although synchronization values of the artificial recordings measured with CorSE changed with the superimposed LFPs, the correlation between the increased synchronization and the increased level of connectivity is somewhat preserved and the behavior between simulations remained similar. Figures 4D,E shows the synchronicities calculated from the artificial recordings with EAPs and added WGN. The results from the different SNR values, 50 and 20%, are presented in Figures 4D,E, respectively, whereas right-hand panels again show corresponding exemplary 2 s signal segments. Adding WGN greatly decreased the synchronization detected by the event based algorithms (c.f., Figures 4D,E). On the other hand, to compare, cES presented better performance in distinguishing the relation between the increasing levels of connectivity and synchronization. Synchronization values calculated by CorSE changed much less due to WGN than those for the event based methods. Even though CorSE was still able to distinguish the unconnected (0% level connectivity) populations from the connected populations, it cannot distinguish the different levels of connectivity; especially with the lower SNR.
In conclusion, CorSE was applicable in determining the level of synchronization according to the level of population connectivity not only for the simulated recordings which solely consisted of EAPs, but also in presence of LFPs. CorSE was also usable for functional connectivity assessment to distinguish the connected and unconnected populations in low SNR conditions caused by high WGN power.

MEA Recordings
We present the results of different synchronization assessment algorithms applied on rat cortical network measurement data (MEA1 and MEA 2) in different forms: First, we demonstrate the obtained synchronization values for the different methods as adjacency matrices (Figures 5, 6). Then, the most synchronized channel pairs found with different methods are presented (Figures 7, 8) on the MEA layout. Finally, we compare the most synchronized channel pairs found by different algorithms (Figures 9, 10). To see the effects of the different spike detection methods, the results employing both spike detection methods (STD and eSTD) are presented for the event based algorithms (Figures 5-10).
Adjacency matrices present an overview of the found synchronicities. For MEA1, the highest synchronization was found between the channels 56 and 67 (represented in Figure 5 as (5,6) and (6,7), respectively, in the sequel similarly) by CorSE, and the TE (with eSTD) and MI (with eSTD) methods (Figures 5B-D). The rest of the methods resulted in different maximally synchronized culture locations. For MEA2, the highest synchronization was found between the channels 73 (7,3) and 83 (8,3) by MI (with STD and eSTD) and cES (with STD and eSTD) methods (Figures 6D,E).
For further comparisons between the algorithms and to demonstrate the most synchronized channel pairs (strongest synchronizations), the 40 strongest synchronizations for MEA1, and the 20 strongest synchronizations for MEA2, are shown in Figures 7, 8, respectively. The different algorithms found most synchronized channel pairs differently. Moreover, the different spike detection methods change the results of the event based synchronization assessment algorithms as can be seen also from Figures 5, 6. Among all algorithms considered, the results of the TE and MI methods are most similar where CorSE has more common outputs with both the TE and MI methods compared to the output of the cES method. For example, there are channels with more links than the others; in other words, network locations acting as hubs (see Figures 7, 8), such as channels 56 and 67 found by the TE and MI methods (with both spike detection methods) and CorSE as the top two channels with the highest numbers of links for MEA1. The cES method found channels 55 and 27 (with eSTD), and 83 and 51 (with STD) as the top two hub channels. For MEA2, channel 75 was a channel that could be considered a hub found by the TE and MI methods (with both spike detection methods) and CorSE, but again not by the cES method.
To assess the common aspects of the results produced with the difference algorithms in more detail, the number of common nodes (channels) and common links (pairwise synchronizations) are plotted according to the number of strongest synchronizations (Figures 9, 10). The results show that for MEA1 and using eSTD, all the algorithms found approximately a similar number of common nodes (Figure 9C), alike for MEA2 with either STD or eSTD (Figures 10C,D). However, the numbers of found common links between the common nodes show more dependence on the method used (Figures 10A,B). The results from the TE and MI methods had the most common links, and CorSE found more common links with these two algorithms than the cES algorithm did.
Finally, we present the results of CorSE in unraveling a developing mouse cortical neuron network. Figure 11 presents the development of a network between the 13th and 29th DIV. Channel pairs with synchronizations exceeding an arbitrary threshold (CorSE > 0.5, selected for illustrative purposes) are illustrated. The first links were seen on the 13th DIV, and thereafter the network gradually expanded while also synchronizations between all the channels were getting stronger, with the strongest network synchronizations found on the 22nd DIV (see Figure 12). The mean values of synchronizations between all the channel pairs (in total 1770 links) were calculated during the development of the network (see Figure 12). The results show that the mean values of all the synchronizations were also correlated with the number of channel pairs with synchronizations greater than CorSE > 0.5; in other words, the overall network synchronization strength followed the same trend with the observed channel pairs. Synchronization between some channels varied for different measurement days. A noteworthy example for this case is the smaller network that appeared on the 25th DIV, could not be observed on the 27th DIV, and reappeared on the 29th DIV with a stronger synchronization (Figure 11).

DISCUSSION
We have previously shown the feasibility of our correlated spectral entropy method CorSE for the assessment of synchronization (Kapucu et al., 2016a). In this paper, we have not only presented the usability of the method with larger real data, but we have made simulations justifying our concept and algorithm. For that, we realized both a toy simulation and simulations based on a widely used computation neuronal network model. The results of the simulations showed that it was possible to distinguish and assess the synchronization and connectivity strengths in neuronal populations by assessing the In practice, in addition to biological sources, many other sources contribute to the measured electrophysiological signals. Any sources of noise that are not biological, e.g., artifacts and electrical interference, might influence the results of the algorithms employed. To assess this, we also performed the simulations with additive white Gaussian noise. Increasing the level of noise naturally had a negative effect on the results, as expected. Clearly, with the increasing noise level, more spikes were missed during spike detection, and the employed event based synchronization measures were more likely to fail. Concerning CorSE, increasing the level of noise had an effect on the uniformity of the spectral distribution, and thus on the synchronization measure. However, simulations showed that even with the increasing noise, it was still possible to distinguish between the connected and unconnected populations, even though the level of connectivity strength was not distinguishable anymore. Here, it may be noted that in all the event based synchronization assessment algorithms in the literature, anything other than detectable spikes is generally considered as biological "noise" (Obien et al., 2015), and a reasonable amount of information from LFPs is omitted. CorSE takes also LFPs into account in synchronization assessment. Thus, we observed the effects of different power ratios of EAPs and LFPs on the detected synchronization: the results showed that CorSE can assess the level of connectivity under the effects of synchronized LFPs. In contrary, since LFPs are considered as biological "noise" by the event based analysis, signals with different power ratios of EAPs and LFPs affect the spike detection accuracy, and thus the results of these algorithms.
For the rat cortical cells used in this study, the mutual information and transfer entropy algorithms showed the most corroborative results, as expected, since the two methods are based on the same theoretical grounds. CorSE exhibited common findings with the mutual information and transfer entropy algorithms, but the results from the corrected event synchronization method were generally different from the results by the rest of the algorithms. Since the corrected event synchronization method is spike rate adaptive by its nature, one possible explanation could be that its results were significantly affected by spike detection.
The similarities and differences of the considered algorithms can be summarized by joint evaluation of adjacency matrices (Figures 5, 6), by observing the differences in the most synchronized MEA channel pairs (Figures 7, 8), and in the common links and nodes (Figures 9, 10): CorSE, and the mutual information and transfer entropy algorithms all found the same link as the most synchronized link for MEA1, whereas for MEA2, the same most synchronized link was found by the mutual information and corrected event synchronization methods. Additionally, channels which had the highest number of links (hub-channels) were identified similarly by CorSE, the mutual information and transfer entropy methods. Consequently, it may well be that the most synchronized links could be found similarly by the different algorithms employed in this work, whereas the most differences might be found in the weaker synchronized links. In fact, the numbers of common links and nodes presented in Figures 9, 10 show similar results for the strongest 10 synchronizations, whereas the difference grows with the increasing number of strongest synchronizations. It is to be noted that since we assessed unguided neuronal cells which developed freely on the MEA plates, there is no ground truth about the network structures which the cell cultures formed. Thus, actual validation of the synchronies measured by different algorithms was impossible with the real data, but has been to an extent provided by the simulations. Still, the comparisons of the findings from the real data give an idea of the general usability of the algorithm, and of its biological plausibility.
Also the feasibility of CorSE to track neuronal development by means of synchronization is studied in this paper. We tracked developing network synchronization for several measurement days using an arbitrary synchronization detection threshold, here, CorSE > 0.5. This provided a clear view to the appearance of the functional network (Figure 11), although setting such a clear-cut threshold might have the effect of the temporary appearance and disappearance of some channels close to the detection threshold in the network development map (see Figure 11, 22-27 days in vitro). The results also correlated with the overall network synchronization behavior (Figure 12).
In conclusion, we have shown that CorSE is a promising tool to assess synchronization in neuronal networks. The method does not possess the shortcomings of event based methods resulting from possible poor event, e.g., spike, detection performance. Moreover, CorSE does not need specified effective frequency bands for the analysis. Our method would be useful especially for the acute analysis of possibly noisy recordings collected with MEAs for the experiments where fast processing is necessary: since it is based on the efficient fast Fourier transform and simple cross correlation, it can be used even for online processing (Semmlow and Griffel, 2014). In fact, for more than a decade, spectral entropy has been utilized in real-time electroencephalogram monitoring to quickly assess the of depth of anesthesia (Viertiö-Oja et al., 2004). We believe that methods FIGURE 11 | The development of a functional neuronal network in a mouse cortical cell culture unraveled by CorSE. The development of the network is shown between the 13th and 29th DIV showing the channels with pair-wise synchronizations CorSE xy > 0.5. not based on events, i.e., methods using all data recorded from the electrodes, such as CorSE, can help us in obtaining more information from the valuable neuronal network measurements, and provide robust synchrony measures, both in vitro and in vivo.
The Matlab code for CorSE has been developed to be applied straight forward on any time series data, and is publicly freely available in the Matlab Central File Exchange (https:// se.mathworks.com/matlabcentral/fileexchange/59626-spectralentropy-based-neuronal-network-synchronization-analysiscorse).

AUTHOR CONTRIBUTIONS
FK developed the method, implemented the program code, designed and made the simulations, contributed to the culturing FIGURE 12 | The number of strong synchronizations and the mean overall synchronization found by CorSE in a developing mouse cortical cell network. The number of synchronizations with CorSE xy > 0.5 (blue) (left vertical axis), and the mean synchronizations between all channels (1770 links in total) (red) found at the difference measurement days (right vertical axis).
of mouse cortical cells, analyzed the data, and wrote the manuscript. IV made the computational neuronal network simulations, contributed to culturing of mouse cortical cells, and to the writing of the manuscript. JM cultured rat cortical cells, collected the data, and contributed to writing of the manuscript. CL contributed to method development and implementation. KL contributed to culturing of mouse cortical cells and writing of manuscript. JT contributed to method development and to the writing of the manuscript. JH contributed to method development and to the writing of the manuscript. All the authors approved the final version to be published.