Burst analysis tool for developing neuronal networks exhibiting highly varying action potential dynamics

This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Burst analysis tool for developing neuronal networks exhibiting highly varying action potential dynamics Kapucu, Fikret E.; Tanskanen, Jarno M. A.; Mikkonen, Jarno; Ylä-Outinen, Laura; Narkilahti, Susanna; Hyttinen, Jari A. K. Kapucu, F., Tanskanen, J., Mikkonen, J., Ylä-Outinen, L., Narkilahti, S., & Hyttinen, J. (2012). Burst analysis tool for developing neuronal networks exhibiting highly varying action potential dynamics. Frontiers in Computational Neuroscience, 6 (38). doi:10.3389/fncom.2012.00038 2012 METHODS ARTICLE published: 19 June 2012 doi: 10.3389/fncom.2012.00038 Burst analysis tool for developing neuronal networks exhibiting highly varying action potential dynamics Fikret E. Kapucu1,2,Jarno M. A. Tanskanen1,2,Jarno E. Mikkonen3, Laura Ylä-Outinen2,4, Susanna Narkilahti2,4 and Jari A. K. Hyttinen1,2* 1Department of Biomedical Engineering, Tampere University of Technology, Tampere, Finland 2Institute of Biosciences and Medical Technology, Tampere, Finland 3Department of Psychology, University of Jyväskylä, Jyväskylä, Finland 4NeuroGroup, Institute of Biomedical Technology, University of Tampere and Tampere University Hospital, Tampere, Finland Edited by: Misha Tsodyks, Weizmann Institute of Science, Israel Reviewed by: Antonio Novellino, ETT S.R.L., Italy Jason Weick, University of Wisconsin-Madison, USA *Correspondence: Jari A. K. Hyttinen, Department of Biomedical Engineering, Tampere University of Technology, P.O. Box 692, FI-33101 Tampere, Finland. e-mail: jari.hytti[email protected] In this paper we propose a firing statistics based neuronal network burst detection algorithm for neuronal networks exhibiting highly variable action potential dynamics. Electrical activity of neuronal networks is generally analyzed by the occurrences of spikes and bursts both in time and space. Commonly accepted analysis tools employ burst detection algorithms based on predefined criteria. However, maturing neuronal networks, such as those originating from human embryonic stem cells (hESCs), exhibit highly variable network structure and time-varying dynamics. To explore the developing burst/spike activities of such networks, we propose a burst detection algorithm which utilizes the firing statistics based on interspike interval (ISI) histograms. Moreover, the algorithm calculates ISI thresholds for burst spikes as well as for pre-burst spikes and burst tails by evaluating the cumulative moving average (CMA) and skewness of the ISI histogram. Because of the adaptive nature of the proposed algorithm, its analysis power is not limited by the type of neuronal cell network at hand. We demonstrate the functionality of our algorithm with two different types of microelectrode array (MEA) data recorded from spontaneously active hESC-derived neuronal cell networks. The same data was also analyzed by two commonly employed burst detection algorithms and the differences in burst detection results are illustrated. The results demonstrate that our method is both adaptive to the firing statistics of the network and yields successful burst detection from the data. In conclusion, the proposed method is a potential tool for analyzing of hESC-derived neuronal cell networks and thus can be utilized in studies aiming to understand the development and functioning of human neuronal networks and as an analysis tool for in vitro drug screening and neurotoxicity assays. Keywords: spike trains, action potential bursts, burst analysis, hESCs, human embryonic stem cells, developing neuronal networks, MEA, microelectrode array INTRODUCTION In this paper we study methods to assess the bursting behavior of developing human neuronal networks. Previously, it has been shown that human embryonic stem cell (hESC)-derived neuronal cellsarefunctionalatsinglecelllevel(Carpenter et al., 2001; Erceg et al., 2008; Lai et al., 2008; Daadi et al., 2009; Bissonnette et al., 2011; Kim et al., 2011) and can form spontaneously functional neuronal networks (Heikkilä et al., 2009). Compared to the more widely studied in vitro neuronal networks, that is, rodent primary cultures, it seems that networking mechanisms and behavior of hESC-derived neurons are more variable in their statistics from individual spikes to bursts (Heikkilä et al., 2009). This calls for new methods for the assessment of network development and functioning, since traditional burst detection algorithms are not in general capable of capturing the bursts and related features of such networks. To assess the functioning of neuronal networks in their different developmental stages and to observe their responses to different drugs, toxins, and chemicals, substrate integrated microelectrode arrays (MEAs) provide an in vitro platform to monitor the firing patterns and the network activity (Gross et al., 1977; Pine, 1980; Wagenaar et al., 2006; Illes et al., 2007; Heikkilä et al., 2009; Ylä-Outinen et al., 2010). Neuronal activity is normally described either by single cell firing called spikes or actual network activity manifested by more or less regular occurring short episodes of intense firing called bursts (Kandel and Spencer, 1961; Connors et al., 1982; Gray and McCormick, 1996). In this type of network activity, neurons are interacting and firing in an orchestrated manner. It is suggested that bursts reflect and influence the plasticity mechanisms and could be used for assessment of network activity (Lisman, 1997). It has been shown, for example, that cultured networks of rat cortical neurons exhibit a significant increase in spontaneous bursting during the development of new synapses and networks (Ichikawa et al., 1993; Maeda et al., 1995; Kamioka et al., 1996). Thus, analysis of bursting behavior is a way to assess the developing neuronal network properties. Frontiers in Computational Neuroscience www.frontiersin.org June 2012 | Volume 6 | Article 38 |1 COMPUTATIONAL NEUROSCIENC E Kapucu et al. Neuronal burst analysis tool Even though bursting is a very fundamental property of the neuronal networks, the definitions of bursts and burst detection methods, however, differ between studies. Some define bursts according to interspike interval (ISI) thresholds and the numbers of spikes in bursts which are set by visual inspection, such as utilizing a fixed ISI of 100 ms and a minimum number of 10 spikes in bursts (Chiappalone et al., 2005), or again utilizing a fixed ISI and a minimum number of spikes in bursts which are chosen according to experimental conditions and differentiate bursts from other activity based on the slopes in time-spike number curves (Turnbull et al., 2005). Others utilize calculated average ISIs of the measurements (Mazzoni et al., 2007), average firing rates and alternatively a fixed ISI threshold of 100 ms (Wagenaar et al., 2006), or logarithmic histogram of ISIs to calculate an ISI threshold for detecting bursts (Selinger et al., 2007; Pasquale et al., 2010). These methods, except that by Turnbull et al. (2005), are focused on analyzing the activity of neurons extracted from rat central nervous system such as rat cortical neurons (Chiappalone et al., 2005; Wagenaar et al., 2006; Pasquale et al., 2010)andrat hippocampal neurons (Mazzoni et al., 2007). Thus, they may be tuned to the type of the analyzed network. In the developing networks, the spiking and bursting may behave differently as the network is formed by active neuronal movement, process formation, and synaptic modulation. In fact, beside the frequent occurrence of “primitive” bursts which are formed by a few spikes, we also observed bursts with tens of spikes and bursts lasting from milliseconds to seconds while studying maturing hESC- derived neuronal networks (Heikkilä et al., 2009). The earlier mentioned most widely applied burst detection and burst analysis methods, however, ignore the primitive and unstable spike train and burst activity of hESC-derived neuronal networks. The spike trains or bursts of such networks are in this paper defined to be “unstable” if their statistics such as the number of spikes forming bursts, ISIs inside and between bursts, burst durations, etc., highly vary. In addition to the need for an applicable burst detection method for developing neuronal networks, it is necessary to obtain characterization measures for the practical analysis of the responses of the networks to different treatments, drugs, toxins, or chemicals. Several parameters such as overall spiking activity, burst frequency, and duration can be used in activity characterization (Bal-Price et al., 2010; Äänismaa et al., 2011; Defranchi et al., 2011; Hogberg et al., 2011). These parameters can be obtained from large data pools by an analysis tool which has no bias for a certain type of analyzed networks, such as fixed burst parameters, e.g., ISI or the number of spikes in bursts. Thus there is a need for methods that provide these parameters also intrinsically from developing networks. Here, we propose a burst detection method without any a priori fixed burst criteria, and demonstrate its applicability with maturing hESC-derived neuronal networks. To demonstrate the need for such methods we illustrate the dynamic nature of hESC-derived neuronal networks during maturation by spike activity maps for different measurement days. Thereafter, we shortly review the existing burst detection methods and compare their performances in the analysis of hESC-derived neuronal network recordings. We compare the applicability of the methods, and finally discuss potential uses of the hereby proposed method in assessing the characteristics of various neuronal networks. MATERIALS AND METHODS CELL CULTURES hESCs [cell line Regea 08/023, passages 36 (used in dataset-I), 42 (used in datasets-II and -III), 44 (used in dataset-IV), and 60 (used in dataset-V)] were differentiated into neuronal cells using the previously published method (Sundberg et al., 2009; Lappalainen et al., 2010) and plated on MEAs as described in Heikkilä et al. (2009). Briefly, 10–15 small aggregates dissected from neurospheres (50,000–150,000 cells in total) were plated on MEA dishes coated with polyethyleneimine (0.05% solution, Sigma-Aldrich, St. Louis, MO, USA) and subsequently with human laminin (20 μg/ml, Sigma-Aldrich). Medium containing basic fibroblast growth factor (4 ng/ml, FGF, Sigma-Aldrich) and brain-derived growth factor (5 ng/ml, BNDF, Gibco Invitrogen, Carlsbad,CA,USAorPeprotech,RockyHill,NJ,USA)was replaced three times a week. The cell seeding area in the MEA was either the normal, that is, 20 mm in diameter or the area was restricted to Ø 4 mm to reduce the amount of cells needed and to guide the cells to grow on top of the electrode area. All the MEAs with cells were kept in an incubator (+37◦C, 5% CO2, 95% air) prior to and between recordings. All recordings were made using MEAs and equipment’s purchased from Multi Channel Systems MCS GmbH (MCS, Reutlingen, Germany). hESC experiments were performed in the Institute of Biomedical Technology (University of Tampere, Tampere, Finland) that has the approval from the Ethics Committee of the Pirkanmaa Hospital District to culture the hESC lines. ELECTROPHYSIOLOGICAL RECORDINGS Electrical activities were recorded using MEAs with square arrays of 59 substrate-embedded titanium nitride microelectrodes (30 μm in diameter, 200 μm inter-electrode distance, model: 200/30iR-Ti-gr, MCS) and internal embedded reference electrodes. Signals were sampled at 20 kHz, and stored to a standard PC using the MC_Rack software (MCS). The culture temperature was maintained at +37◦CusingaTC02tempera- ture controller (MCS) during the measurements. Recordings were visually inspected for artifacts and the measurements or channels likely to contain artifacts were excluded from the further analysis. In this paper, for demonstrating the proposed analysis methods, we opted for the analysis of artifact free data only. Spike detection was carried out online by setting an amplitude threshold at six times the standard deviation of the signal noise level and the spike time stamps and spike waveform cutouts were stored in the MC Rack software. For method validation, we utilized five different data sets (Ds-I, Ds-II, Ds-III, Ds-IV, Ds-V) that altogether contained measured data from 27 MEAs (each containing 59 electrodes referred to hereafter as channels or ch). Each MEA was measured altogether three times (mes1, mes2, and mes3) and each measurement lasted approximately 300 s. The first measurement day was chosen according to the criteria that at least 10% of the channels in a MEA were active and in active channels at least 100 spikes Frontiers in Computational Neuroscience www.frontiersin.org June 2012 | Volume 6 | Article 38 |2 Kapucu et al. Neuronal burst analysis tool were found during the recording period. The first measurement day (mes1) varied between 6 and 22 days of culturing in MEA. Thereafter, the developing networks were measured at 4–7 days after the first measurement (mes2), and the third measurement day (mes3) was 2–4 days after the second measurement day (Table 1). SPIKE ACTIVITY COLORMAPS To observe how dynamic the firing was during the maturation process of hESC-derived neurons, a spike activity map (colormap) of all 59 channels for different measurement days were formed. In our colormaps, MEA channels are represented as an 8×8 matrix and the layout of the MEA matches the colormap matrix elements spatially. Since a MEA has 60 channels (including the reference channel), the corner elements of the colormap matrix have no values. Thus, the color shown in the corners of the subfigures of Figure 1 is set based on values interpolated from the neighboring channels. At first, spikes were counted separately for the different measurement days, yielding the total spike counts of each measurement for every MEA channel. Secondly, the logarithms of these spikecountswerecalculatedtobeabletoshowawiderange of values in the same colormap, and the values were mapped to color, the colors were interpolated between the channels, and the colormaps were contoured. ISI AND ISI HISTOGRAM BASED BURST DEFINITIONS Time intervals between consecutive spikes, i.e., ISIs are very commonly used in the analysis of neuronal recordings. ISI has also been used as one of criteria to detect the burst activity in some reported algorithms (Chiappalone et al., 2005; Wagenaar et al., 2006; Mazzoni et al., 2007). In these algorithms, an ISI thresholdisselectedorcalculatedandafixednumberofconsecutive spikes (three, four, or ten, depending on the algorithm) with ISIs less than the selected threshold are considered burst spikes. On the other hand, an ISI histogram can be easily formed after spike detection by counting the spikes and time binning ISIs. It can be calculated for, e.g., each channel, recording, measurement day or complete dataset. Although ISI histograms are capable of representing some characteristics of the firing activities, it is usually not rewarding to analyze them alone. Gradual decay in the ISI histogram, forming a tail after the peak of the ISI histogram and fluctuations of local extrema (local minima and maxima) are common observations in ISI histograms. The analysis of such histogram behavior is a promising method for network analysis, although some aspects of network characteristics are not easily observable in raw ISI histograms. LOGARITHMIC INTERSPIKE INTERVAL HISTOGRAM (logISIH) ALGORITHM One of the alternative solutions for analyzing ISI histograms is plotting the logarithmic ISI instead of plain ISI, which was proposed previously by Selinger et al. (2007) and employed in ratcorticalcellmeasurementanalysisbyPasquale et al. (2010). The method is based on plotting ISI histograms using logarithmic instead of linear scale. In the algorithm, ISI threshold is selected at the point where intra-burst ISI is most clearly separated from inter-burst ISI. Clear separation is indicated by the distinct principal and secondary peaks formed in the logarithmic histogram representing intra and inter-burst ISIs, respectively (see Figure 7A). Parameter “void,” which is described in detail by Selinger et al. (2007), is calculated to assess this separation. In logISIH algorithm, if the ISI threshold is lower than 100 ms, bursts are detected only according to this threshold; otherwise a fixed ISI threshold of 100 ms is employed in finding burst cores and the calculated threshold is employed to detect burst boundary spikes. Complete detected bursts are formed by combining the detected burst cores with the boundary spikes adjacent to each burst core. The algorithm has strict rules such as the occurrence of the first peak (principal peak) in the ISI histogram, which represents intra-burst ISIs, within a defined time window (here 100ms). Recordings should have at least one peak at an ISI less than 100ms in their logarithmic ISI histograms. In the case where the logarithmic histograms of the recordings have no good separation between inter-burst and intra-burst ISIs or there is no principal peak before 100 ms, another strict burst definition is employed, such as requiring 10 spikes in a row which have ISIs less than 100 ms, as in the study by Chiappalone et al. (2005). As a result, this algorithm was not very suitable for the analysis of our data as bursts with 10 or more spikes are hardly available. To make the algorithm more comparable to our algorithm and applicable to our data, we modified the algorithm to consider three spikes in a row instead of 10. CUMULATIVE MOVING AVERAGE METHOD FOR DETECTING BURSTS For analyzing the recordings for which no clear separation can be observed between the inter- and intra-burst ISIs, like the majority of the recordings obtained from maturing hESCs (see Figure 8A), we propose an adaptive method based on the cumulative moving average (CMA) of the ISI histogram. As an alternative for analyzing the raw ISI histogram, CMA of the ISI histogram allows ustoobservethecumulativeaveragespikecountuptoaparticular ISI. Generally, CMA of a data series smoothens short term fluctuations and highlights long term trends. In our particular case, CMA of the ISI histogram provides us the general change Table 1 | The measurement days and the number of included MEAs for different data sets. Ds-I (7 MEAs) Ds-II (6 MEAs) Ds-III (7 MEAs) Ds-IV (6 MEAs) Ds-V (6 MEAs) mes1 (days after plating) 20 7 22 6 10 mes2 (days after plating) 25 11 27 12 17 mes3 (days after plating) 28* 15 29 15 20 ∗Analysis based on six MEAs instead of seven. Frontiers in Computational Neuroscience www.frontiersin.org June 2012 | Volume 6 | Article 38 |3 Kapucu et al. Neuronal burst analysis tool FIGURE 1 | Colormaps of the spike activities. Spikes were counted from all three measurements time points (mes1, mes2, mes3) from the individual MEAs from three different data sets (A) Ds-V, (B) Ds-IV, and (C) Ds-III. Logarithmic values of the counts are mapped to colors to better facilitate observations, dark blue representing silent channels (0 spikes) and dark red the most active ones (1000 spikes). in the trend of the ISI histogram, and allows us to define an ISIvaluewhichwecanuseasathresholdtodefinethebursts in a particular recording. Thus, without considering any local changes we can identify the ISIs at which critical changes occur, i.e., the ISI at which the average spike count starts decreasing. The proposed algorithm has no strictly fixed parameters that would render any particular types of neuronal network behaviors not analyzable. Thus, every recording is evaluated based on its inner dynamics and bursts are detected for further network analyzes and characterization. The proposed CMA algorithm was implemented in Matlab environment for post recording analysis and consists of the steps described in the following three subchapters. Calculating ISI threshold for burst like patterns Let yi,i=1,...,N,withNthe total number of ISI bins, be the spike count in the ith ISI bin. The value of the cumulative sum of the histogram CHIat the Ith, I≤N,ISIbinisdefinedas CHI= I  i=1 yi(1) The corresponding CMA is given by CMAI=1 I I  i=1 yi(2) whose maximum, CMAm, is reached at the mth ISI bin, and m=arg max k=1,...,N1 k k  i=1 yi(3) This point represents the maximum that the average spike count reaches. ISI threshold for defining a burst might, for example, be selected at this maximum point after which CMA begins to decrease. However, adding a tolerance to this maximum, i.e., selecting the actual ISI threshold as α·CMAm,where 0<α<1, strengthens the burst detection. Here, αis selected according to the ISI behavior of the ISI histogram. Generally, the ISI histogram of a burst containing recording and its CMA curve exhibit a peak at lower ISI values and a tail at higher ISI values. Intra-burst ISIs are expected to be in the neighborhood of the peak because of the fact that intra-burst ISIs are shorter than the ISIs of individual spikes which don’t belong to a burst. If they exist, individual spikes are located in the tail of the histogram, whereas burst tails or pre-burst spikes (burst related spikes) are located in the histogram between ISIs of intraburst spikes and ISIs of individual spikes. For a simple bursting model, in which ISI values do not widely vary, ISI values will form an almost symmetrical distribution with a short tail, i.e., the skewness of the distribution is approximately zero, and most of the burst spikes are expected to be located in the vicinity of CMAm(c.f., Figure 2A).Thus,wecandefinetwoparameters:α1 which can be set close to one, and for the burst related spikes α2<α1. However, for most of the burst models ISI histograms lean to the right with a tail, i.e., they are positively skewed. In these cases, ISI values have large variance, and also the tail of the ISI histogram contains intra-burst ISIs (c.f., Figures 2B,C). Longer the tail, more skewed the histogram, and the smaller α can be set not to miss the intra-burst ISIs in the tail. Denoting the ISI at CMAm(3) by xm, the ISI threshold xt,xt>xm,for burst detection is found at the mid time point of the ISI bin for which the value of the CMA curve is the closest to α·CMAm. Frontiers in Computational Neuroscience www.frontiersin.org June 2012 | Volume 6 | Article 38 |4 Kapucu et al. Neuronal burst analysis tool FIGURE 2 | Simulation of different bursting models and their ISI histograms with the corresponding calculated CMA curves and skewness values. (A) A simple model whose ISI values don’t widely vary. ISI distribution is almost symmetrical with a short tail and an approximate skewness value of 0.2. On the left hand side, the ISI histogram is shown with gray bars and the corresponding CMA curve is shown with black line. On the right hand side, the bursts detected using an αvalue which corresponds to the skewness value are labeled by black lines. Black circles denote the burst spikes whereas red crosses represent the burst related spikes [similarly for (B) and (C)]. (B) A bursting model with a wider tail in its ISI histogram with an approximate skewness value of 1.2. (C) A bursting model with a wider ISI distribution than in the previous models. Consequently, its histogram has a longer tail and a relatively higher skewness value of approximately 3. FIGURE 3 | Selecting thresholds for burst ISIs by using CMA curve. (A) The ISI histogram (gray bars), cumulative histogram (dash-dotted), and the corresponding CMA curve (solid). Vertical axis is logarithmic. (B) Maximum value of CMA curve, CMAm, is reached at the ISI xmand the ISI threshold xtfor bursting was found at the ISI corresponding to the CMA value closest to α·CMAm. The ISI threshold for burst detection is marked with red dash-dotted vertical line. Vertical axis is linear and thus the cumulative histogram is not shown. In Figure 3, CMA of the ISI histogram and the calculation of the threshold are illustrated. After calculating the threshold, burst detection is employed, here with the requirement of at least three spikes in a row (triplet) with ISIs below the calculated ISI threshold. We consider that extracellular burst means more event than one, as also noted by Lisman (1997)andIzhikevich et al. (2003), thus we concentrated on triplets to be sure for this study. However, also two spikes in a row (duplet) instead of a triplet can be accommodated in the algorithm’s if the user prefers this option. Frontiers in Computational Neuroscience www.frontiersin.org June 2012 | Volume 6 | Article 38 |5 Kapucu et al. Neuronal burst analysis tool Calculating ISI threshold for burst related spikes The term burst related spikes is used here for defining the burst tails or pre-burst spikes which are located in the neighborhood of bursts and are an essential part of the bursting behavior. In a previous study, changes in the firing activity before and after bursts were observed and analyzed (Wang and Hatton, 2005). Burst tails have also been studied and utilized in classifying bursting behaviors (Wagenaar et al., 2006). It is common to observe the burst related spikes in our recordings as well. To calculate ISI threshold for the burst related spikes with α2, the first step of the CMA algorithm is repeated. After detecting the potential burst related spikes, the ones which are not following or followed by a burst are omitted. To automatically select αvalues, we can form a scale of αvalueswiththecorrespondingskewnessvalues.Thisrelationshipis depicted in Figure 4 with α1values used for detecting burst spikes and α2values used for detecting burst related spikes. The scale was formed by experimenting with the relation of αand skewness for our recordings. Extending bursts with the burst related spikes and merging close bursts The burst related spikes are included in their neighboring bursts in this step. Also, the bursts which are closer to each other than the threshold calculated in the second step are merged together. This step also corrects erroneous burst spike detections and misses caused by ISI variance inside bursts and is especially advantageous for the analysis of maturing networks, which frequently have high ISI variance and consist of both bursts and individual spikes. Figure 5 demonstrates data with high intra-burst ISI variability. Black circles are detected burst spikes whereas red crosses are the detected burst related spikes after the first step of CMA algorithm, whereas black lines indicate bursts after the extending and merging process. In Figure 5,theburstaround the 50th second has an intra-burst ISI variability of approximately 600 ms and the burst related spikes are detected during the burst. Extending bursts to the burst related spikes and merging close bursts, we get the satisfactorily detected burst marked with black line instead of the erroneously detected two separate bursts. ANALYSIS OF THE DETECTED BURSTS After detecting the bursts by CMA algorithm, all the bursts from five data sets were pooled according to their measurement days. Additionally, we pooled the bursts according cell seeding area size on the MEAs. For analyzing the bursts, we calculated three parameters from the detected bursts for each channel: average burst duration (ABD), total number of bursts detected (TNB), and average number of spikes per burst (ASpB). The relations of these parameters were plotted for different measurement days to analyze the changes during the development. RESULTS A colormap which shows the spike activities is a useful tool for viewing the channel dynamics of hESC-derived neurons during the maturation process. The activity colormaps for the three measurement days (mes1, mes2, and mes3) of three MEAs (N732, N752, and N878) from three data sets (Ds-III, Ds-IV, and Ds-V) are presented in Figure 1.Thedynamicsofthenetworkscanbe observed from the fading and rising of the firing activities in time. Before proceeding to compare the burst detection results using the previously published methods (Chiappalone et al., 2005; Pasquale et al., 2010) and our method, we investigated if it was appropriate to set one fixed threshold to define intra-burst ISIs for the burst detection algorithms. Since the principal peak in a logarithmic histogram represents mostly the intra-burst ISI values as mentioned earlier, the principal peak ISI values of our recordings were calculated to observe the feasibility of using a fixed threshold as in previously published algorithms. In Figure 6 are shown the logarithmic ISI histograms of the selected active channels whose principal peaks of the logarithmic histograms include at least 30 spikes to make the three different common types of logarithmic ISI histogram trends found in hESC-derived neuronal recordings clearly observable. Histograms with two well separated peaks, with only one peak and no local extrema, and with local extrema which cannot be separated well enough are typical in our recordings. Figure 6A demonstrates that both the principal ISI values and the histogram shape may be greatly varying. To further show the variability of the ISI value of the peak, logarithmic ISI histograms for all of the channels from every recording were calculated and the locations FIGURE 4 | Scale of αvalues and the corresponding ISI distribution skewness values. The scale is formed by experimenting with the relation of αand skewness values for our recordings. α1values were set for burst spikes as 1, 0.7, 0.5, and 0.3 for the skewness values less than 1, 4, 9, and more than 9, respectively. α2values were set for burst related spikes as 0.5, 0.3, and 0.1 for the skewness values less than 4, 9, and more than 9, respectively. Frontiers in Computational Neuroscience www.frontiersin.org June 2012 | Volume 6 | Article 38 |6 Kapucu et al. Neuronal burst analysis tool FIGURE 5 | Bursting data with high intra-burst ISI variability. Black circles are the detected burst spikes whereas red crosses are the detected burst related spikes after the first step of CMA algorithm. Black lines indicate the bursts after the extending and merging process. FIGURE 6 | The logarithmic ISI histograms and their principal peak locations for the data selected to illustrate the different commonly encountered cases. (A) Histograms with two well separated peaks (blue dashed-dotted and green dashed lines), with only one peak and no local extrema (black dotted line), and with local extrema which cannot be separated well enough (purple solid line), are typical for ISIs of hESC-derived neuronal recordings. (B) Principal peak locations of selected channels vary from approximately 30–2000 ms. of the principal peaks consisting of the minimum of 30 spikes are shown in Figure 6B. The peak locations have an approximate ISI range from 30 to 2000 ms. Figure 7A demonstrates the logarithmic ISI histogram of the data from one channel of a MEA (MEA N728, Ds-II, mes2) which exhibits clearly separated bursts shown in Figure 7B.Inthecase shown in Figure 7, it can be seen that the principal peak (red circle) and the secondary peak (red cross) are well separated and the ISI threshold for the burst boundary spikes is calculated at 1259 ms (blue dashed line). On the other hand, the location of the principal peak is almost at the proposed threshold for the burst cores at 100 ms. Accordingly, we experimented with the ISI thresholds for the burst cores at 100, 200, and also at 1000 ms, which seems more rational when observing the Figure 7A. The burst detection results labeled with “1” in Figure 7B represent the results by the logISIH algorithm with the threshold ISI of 100 ms as also previously shown (Pasquale et al., 2010). The results obtained by changing the threshold to 200 ms for the principal peak locations Frontiers in Computational Neuroscience www.frontiersin.org June 2012 | Volume 6 | Article 38 |7 Kapucu et al. Neuronal burst analysis tool FIGURE 7 | The logarithmic ISI histogram of a channel with well separated inter- and intra-burst ISIs. (A) Logarithmic histogram (black solid line) has two different peaks where red circle shows the principal peak and the red cross shows the secondary peak. Threshold for the burst boundary spikes is calculated to occur at 1259 ms (blue dashed line). (B) Burst detection results of logISIH method by employing the burst core and principal peak thresholds at 100, 200, and 1000 ms are labeled as 1, 2, and 3 respectively. The result of the CMA algorithm is labeled as 4. FIGURE 8 | The logarithmic ISI histogram of a channel with poorly separated inter- and intra-burst ISIs. (A) Logarithmic histogram (black solid line) has no well-definable principal or secondary peaks. (B) Burst detection results obtained with pre-defined thresholds. The results obtained by employing a threshold for intra-burst spikes of 100 ms with at least 10 spikes in a burst, 100 ms with at least three spikes, and 200 ms with at least three spikes, are labeled as 1, 2, and 3, respectively. The result of the CMA algorithm is labeled as 4. and for the burst cores are shown labeled “2.” The results obtained by changing ISI threshold for the burst cores from 100 to 1000 ms and also changing the threshold for the principal peak locations to 1000 ms are shown labeled “3.” The results given by our CMA algorithm are shown labeled “4” in Figure 7B.SkewnessoftheISI distribution was in this case found to be 4.7, which corresponds to α1=0.5andα2=0.3. Figure 8A shows an example of logarithmic ISI histogram of a recording in which the bursts are not clearly separable. As can be seen, results cannot be obtained by logISIH algorithm since the required criterion for the separation cannot be satisfied. Instead, we used fixed ISI threshold values and fixed number of burst spikes as the criteria to compare the methods. In Figure 8B,50sof the recording and the results are shown to better observe the burst Frontiers in Computational Neuroscience www.frontiersin.org June 2012 | Volume 6 | Article 38 |8