# Predicting Research Trends in Artificial Intelligence with Gradient Boosting
Decision Trees and Time-aware Graph Neural Networks

Yichao Lu  Layer 6 AI
Toronto, Ontario, Canada
<EMAIL_ADDRESS>
###### Abstract

The Science4cast 2021 competition focuses on predicting future edges in an
evolving semantic network, where each vertex represents an artificial
intelligence concept, and an edge between a pair of vertices denotes that the
two concepts have been investigated together in a scientific paper. In this
paper, we describe our solution to this competition. We present two distinct
approaches: a tree-based gradient boosting approach and a deep learning
approach, and demonstrate that both approaches achieve competitive
performance. Our final solution, which is based on a blend of the two
approaches, achieved the $1$st place among all the participating teams. The
source code for this paper is available at
https://github.com/YichaoLu/Science4cast2021.

###### Index Terms:

semantic network, link prediction, network analysis, feature engineering,
graph neural network

##  I Introduction

A machine learning model that is able to predict the future research trends
can greatly facilitate advanced research in a variety of disciplines [1, 2].
To date existing literature on predicting future research trends with machine
learning remains quite limited, primarily due to the lack of a large-scale,
authoritative benchmark dataset.

The Science4cast 2021 competition 111https://www.iarai.ac.at/science4cast/
aims at addressing this challenge by providing a high-quality benchmark
dataset that is an order of magnitude larger than existing semantic network
datasets. The dataset used in the competition consists of a growing semantic
network with $8,766$ snapshots, each corresponding to a day between the
beginning of 1994 and the end of 2017. Each vertex in the semantic network
represents an artificial intelligence concept, and an edge is drawn between a
pair of vertices when the two concepts have been investigated together in a
scientific paper. Participants in this competition are asked to build a
machine learning model that can predict the future edges of this evolving
semantic network.

In this paper, we present our solution to this competition. We explored two
approaches: a tree-based gradient boosting approach and a deep learning
approach. We demonstrate that both approaches achieve highly competitive
performance. Our final solution, which is based on a blend of the two
approaches, further boosts the performance, and it achieved the $1$st place
among all the participating teams.

##  II Approach

In this section, we describe our solution to the Science4cast 2021
competition. We begin by introducing the task of the competition, followed by
a description of our dataset split strategy. Then we describe our feature
engineering process. After that, we present two distinct approaches: a tree-
based gradient boosting approach and a deep learning approach. Our final
solution is based on a blend of the two approaches.

###  II-A Task Formulation

Let $G_{t_{1}}=(V,E_{t_{1}})$ be an undirected graph, where $V$ is the set of
vertices and $E_{t_{1}}\subseteq V\times V$ is the set of observed links until
time $t_{1}$. Each vertex represents an artificial intelligence concept, and
edges between vertices are drawn when two concepts are investigated together
in a scientific paper. Among vertex pairs that are unconnected by time
$t_{1}$, we aim at predicting which edges between unconnected vertex pairs are
more likely to form by time $t_{2}$.

The following notation is used throughout the paper. Let $u,v\in V$, the
binary adjacency matrix of $G$ is denoted as $A$, where $A_{u,v}=1$ if
$(u,v)\in E$ and $A_{u,v}=0$ otherwise. For any vertex $u\in V$, let
$\Gamma(u)$ be the $1$-hop neighbors of $u$, and $|\Gamma(u)|$ be the number
of $1$-hop neighbors of $u$.

TABLE I: The list of features and their importance scores. Name | Description | Importance
---|---|---
rank_num_neighbors_diff_2_year_u | The rank of the number of new neighbors since $2$ years prior to $t_{1}$ for vertex $u$ | 217753490.51
rank_num_neighbors_diff_2_year_v | The rank of the number of new neighbors since $2$ years prior to $t_{1}$ for vertex $v$ | 208691110.58
rank_num_neighbors_diff_1_year_v | The rank of the number of new neighbors since $1$ year prior to $t_{1}$ for vertex $v$ | 124923394.68
rank_num_neighbors_diff_1_year_u | The rank of the number of new neighbors since $1$ year prior to $t_{1}$ for vertex $u$ | 108458621.54
rank_num_neighbors_diff_3_year_u | The rank of the number of new neighbors since $3$ years prior to $t_{1}$ for vertex $u$ | 73609968.95
jaccard_index | The Jaccard index between vertices $u$ and $v$ at time $t_{1}$ | 63507171.32
rank_num_neighbors_diff_3_year_v | The rank of the number of new neighbors since $3$ years prior to $t_{1}$ for vertex $v$ | 61741280.01
jaccard_index_2000 | The Jaccard index between vertices $u$ and $v$ since $2020$ | 46164742.85
pagerank_score_v | The PageRank score at time $t_{1}$ for vertex $v$ | 18960451.79
pangrank_score_u | The PageRank score at time $t_{1}$ for vertex $u$ | 14231787.78
v | The ordinal encoding for vertex $v$ | 12303301.32
u | The ordinal encoding for vertex $u$ | 12156829.26
jaccard_index_1_year | The Jaccard index between vertices $u$ and $v$ at at $1$ year prior to $t_{1}$ | 10217562.38
jaccard_index_diff_3_year | The difference between the Jaccard indices at time $t_{1}$ and at $3$ years prior to $t_{1}$ | 1164973.71
jaccard_index_diff_2_year | The difference between the Jaccard indices at time $t_{1}$ and at $2$ years prior to $t_{1}$ | 784825.26
pagerank_score_diff_2_year_v | The different between the PageRank scores at time $t_{1}$ and at $2$ years prior to $t_{1}$ for vertex $v$ | 347271.85
pagerank_score_diff_2_year_u | The different between the PageRank scores at time $t_{1}$ and at $2$ years prior to $t_{1}$ for vertex $u$ | 296226.81
rank_num_neighbors_v | The rank of the number of neighbor at time $t_{1}$ for vertex $v$ | 213892.33
rank_num_neighbors_u | The rank of the number of neighbor at time $t_{1}$ for vertex $u$ | 213258.87
rank_num_neighbors_2000_v | The rank of the number of new neighbors since $2000$ for vertex $v$ | 190724.69
rank_num_neighbors_1_year_v | The rank of the number of neighbor at $1$ year prior to $t_{1}$ for vertex $v$ | 184202.64
rank_num_neighbors_2_year_v | The rank of the number of neighbor at $2$ years prior to $t_{1}$ for vertex $v$ | 182078.28
rank_num_neighbors_1_year_u | The rank of the number of neighbor at $1$ year prior to $t_{1}$ for vertex $u$ | 171340.45
rank_num_neighbors_3_year_u | The rank of the number of neighbor at $3$ years prior to $t_{1}$ for vertex $u$ | 162496.88
rank_num_neighbors_2_year_u | The rank of the number of neighbor at $2$ years prior to $t_{1}$ for vertex $u$ | 155750.99
rank_num_neighbors_3_year_v | The rank of the number of neighbor at $3$ years prior to $t_{1}$ for vertex $v$ | 153530.45
rank_num_neighbors_2000_u | The rank of the number of new neighbors since $2000$ for vertex $u$ | 143596.41
jaccard_index_diff_1_year | The difference between the Jaccard indices at time $t_{1}$ and at $1$ year prior to $t_{1}$ | 120941.88
pagerank_score_diff_1_year_v | The different between the PageRank scores at time $t_{1}$ and at $1$ year prior to $t_{1}$ for vertex $v$ | 103631.39
pagerank_score_diff_1_year_u | The different between the PageRank scores at time $t_{1}$ and at $1$ year prior to $t_{1}$ for vertex $u$ | 100832.93
jaccard_index_3_year | Jaccard index between vertices $u$ and $v$ at $3$ years priors to $t_{1}$ | 98789.44
pagerank_score_diff_3_year_u | The different between the PageRank scores at time $t_{1}$ and at $3$ years prior to $t_{1}$ for vertex $u$ | 86917.48
pagerank_score_diff_3_year_v | The different between the PageRank scores at time $t_{1}$ and at $3$ years prior to $t_{1}$ for vertex $v$ | 83766.08
pangrank_score_2_year_u | The PageRank score for vertex $u$ at $2$ years prior to $t_{1}$ | 43588.77
pagerank_score_1_year_v | The PageRank score for vertex $v$ at $1$ year prior to $t_{1}$ | 42649.46
pangrank_score_1_year_u | The PageRank score for vertex $u$ at $1$ year prior to $t_{1}$ | 39638.87
jaccard_index_2_year | The Jaccard index between vertices $u$ and $v$ at $2$ years priors to $t_{1}$ | 26012.81
pagerank_score_2_year_v | The PageRank score for vertex $v$ at $2$ years prior to $t_{1}$ | 11689.22
pagerank_score_3_year_v | The PageRank score for vertex $v$ at $3$ years prior to $t_{1}$ | 9604.56
pangrank_scores_3_year_u | The PageRank score for vertex $u$ at $3$ years prior to $t_{1}$ | 7482.16
jaccard_index_diff_2000 | The difference between the Jaccard index at time $t_{1}$ and the Jaccard index since $2000$ | 2.85

###  II-B Dataset Split

In the competition test set, $t_{1}$ is set to 20171231 and $t_{2}$ is set to
20201231. Equivalently, the task asks for the prediction of which scientific
pairs of concepts will be investigated by scientists over $3$ years. Our
dataset split strategy is to create the training set and the validation set
that mimic the test set. This ensures a high correlation of performance of the
model on the validation set and the test set. Specifically, for the training
set we set $t_{1}$ to 20111231 and set $t_{2}$ to 20141231, while for the
validation set we set $t_{1}$ to 20141231 and set $t_{2}$ to 20171231.

For the training set and the validation set, we sample $N$ vertex pairs that
are unconnected by time $t_{1}$ and the labels are determined by whether edges
between unconnected pairs will form by time $t_{2}$. We set $N$ to be
$10,000,000$ both for the training set and the validation set. While a larger
training set typically leads to better performance, there exists a trade-off
between performance and efficiency since a larger training set also incurs
heavier computational cost. We find setting $N$ to $10,000,000$ results in a
good balance between performance and efficiency.

The competition dataset is highly imbalanced with only $0.178\%$ newly
connected edges ($3,724,986$ new edges out of $2,092,398,869$ unconnected
vertex pairs) formed between $2014$ and $2017$. To deal with the scarcity of
positive samples, we include all positive samples when constructing the
training set and the validation set. This results in $1,275,503$ positive
samples in the training set, and $3,724,986$ positive samples in the
validation set. The negative samples are then uniformly sampled from the
remaining unconnected vertex pairs. We experimented with alternative sampling
strategies for negative samples, e.g., sampling unconnected vertex pairs
according to vertex degrees to include more “hard negatives” in the dataset
[3], but did not see gains as expected. We also experimented with a data
augmentation approach, which exploits the order invariance of the task, i.e.,
the model should output the same link prediction score for the vertex pair
$(u,v)$ and the vertex pair $(v,u)$. Briefly speaking, we double the size of
the training set by swapping the order of the vertices in each pair. We found
that this data augmentation approach only brings about marginal gains in terms
of performance at the cost of much longer training time.

Figure 1: High-level Architecture of the Time-aware Graph Neural Network.

###  II-C Feature Engineering

Feature engineering plays an important role in the winning entries of various
machine learning competitions [4, 5, 6]. While deep neural networks have
become the de facto approach to representation learning in domains such as
computer vision [7, 8, 9] and natural language processing [10, 11], machine
learning models still rely heavily on feature engineering with domain
expertise when it comes to dealing with tabular data [12].

We conducted extensive feature engineering with the aim of capturing the
centralities of the vertices, the proximity between vertex pairs, and their
evolution over time. The final set of features that we used can be partitioned
into two groups: vertex features and pairwise features, which we describe in
what follows.

Vertex features are derived using the number of neighbors and the PageRank
score [13] of the vertex, which reflect the vertex centrality. In order to
model the evolution of the vertex centrality over time, we additionally
compute the centrality of the vertex (in terms of the number of neighbors and
the PageRank score) at $1$, $2$, and $3$ years prior to $t_{1}$. The change of
centrality of the vertex is captured by taking the difference between the
centrality measures at time $t_{1}$ and that at $1$, $2$, and $3$ years prior
to $t_{1}$. We also compute the centrality measures on a subset of the full
semantic network, where edges that were formed prior to $2000$ are discarded,
which encourages the centrality measures to be more focused on more recent
edges. For vertex features derived using the number of neighbors of the
vertex, we apply the rank transformation, which replaces the features with
their ranks, or average ranks in case of ties. This is because the number of
neighbors of the vertex is constantly increasing over time, and we need a
normalization procedure to adjust these features so that the features within
different dataset splits are comparable. Rank transformation works well
because it is a monotonic transformation that preserves the order among the
features and the distribution of the rank transformed features remains
unchanged over time. Vertex features derived from the PageRank score of the
vertex do not require the rank transformation since PageRank scores are
already normalized.

Pairwise features are derived using the Jaccard index [14], which is defined
as the number of common neighbors between $u$ and $v$ divided by the total
number of neighbors of $u$ and $v$:

$J(u,v)=\left\\{\begin{aligned}
&\frac{|\Gamma(u)\cap\Gamma(v)|}{|\Gamma(u)\cup\Gamma(v)|},&|\Gamma(u)\cup\Gamma(v)|>0\\\
&\text{NaN},&|\Gamma(u)\cup\Gamma(v)|=0\end{aligned}\right.,$ (1)

where we use NaN to indicate missing values. The number of common neighbors
between $u$ and $v$ can be obtained by computing the square of the adjacency
matrix:

$|\Gamma(u)\cap\Gamma(v)|=(A^{2})_{u,v}.$ (2)

A larger Jaccard index indicates a higher proximity between vertices, and thus
a higher chance of forming edges between them. Similar to vertex features, we
model the evolution of the proximity between vertex pairs by computing the
Jaccard indices at $1$, $2$, and $3$ years prior to $t_{1}$ and their
differences to the Jaccard index at time $t_{1}$.

The list of all features that we used in this competition is presented in
Table I. We additionally show the importance score for each feature and sort
the features accordingly. The importance score for each feature is obtained
from a trained Light Gradient Boosting Machine (LightGBM) [15] by computing
the total gains of splits that use the feature. We can see that the features
corresponding to the number of new neighbors of the vertices $1$, $2$, and $3$
years prior to $t_{1}$ are the most important ones. This is because these
features are very strong indicators of the current research trends. In
addition, the Jaccard index is also a very effective feature since it measures
the proximity between the vertex pairs. We can also see that the PageRank
score is more useful than the number of neighbors of the vertex as it is a
better metric of vertex centrality.

###  II-D A Gradient Boosting Approach

In our tree-based gradient boosting approach, we use the Light Gradient
Boosting Machine (LightGBM), which is a highly efficient implementation of the
gradient boosting decision tree (GBDT) algorithm. We experimented with other
GBDT implementations, such as XGBoost [16] and CatBoost [17], and found that
LightGBM achieves the best performance.

We train LightGBM to minimize the binary cross entropy loss between its
predictions and the ground truths. LightGBM has also implemented the
LambdaMART algorithm [18], which is a state-of-the-art Learning-to-Rank (LTR)
technique [19]. Training with LambdaMART results in better validation AUC than
training with the binary objective but does not work well in the blending
stage. This is because LTR models output uncalibrated scores. Therefore we did
not use it in our final solution. In order to combat overfitting due to the
scarcity of positive samples, we set num_leaves to 16 and set max_depth to 4
to limit the growth of the trees. We also set bagging_fraction to 0.8 and set
feature_fraction to 0.9, which means that we randomly select $90\%$ of the
data and $80\%$ of the features before training each tree. We train LightGBM
for $10,000$ iterations with a learning rate of $0.01$, and stop training when
the validation AUC does not improve in the last $100$ iterations.

Since decision trees are invariant to monotonic transformations of features
and can handle missing values without imputation, we directly feed raw
features into LightGBM. Note that we directly feed the ordinal encoding of the
vertices $u$ and $v$ into LightGBM, despite the fact that they are essentially
categorical features. A common practice to deal with categorical features in
GBDT algorithms is to one-hot encode the features. However, given the high
cardinality of the vertices, we found that ordinal encoding results in the
best performance.

###  II-E A Deep Learning Approach

In our deep learning approach, we employ a time-aware graph neural network to
learn vertex representations on dynamic semantic networks. The learned vertex
representations for the vertex pair are then concatenated with the engineered
features and fed into a multi-layer perceptron (MLP) which outputs the final
link prediction score; see Figure 1.

We propose to apply an exponential decay on the edge weights of the semantic
network so as to encourage the graph neural network to pay more attention to
more recently formed edges. This also improves the time-awareness of the graph
neural network since it can now recover the exact time that a particular edge
is formed from the edge weight. The adjacency matrix with exponentially
decayed weights $\hat{A}$ is defined as

$\hat{A}_{u,v}=\left\\{\begin{aligned} &e^{-0.0001(t_{1}-t_{u,v})},&(u,v)\in
E\\\ &0,&(u,v)\notin E\end{aligned}\right.,$ (3)

where $t_{u,v}$ denotes the most recent time prior to $t_{1}$ that an edge
between $u$ and $v$ was formed, and $t_{1}-t_{u,v}$ denotes the number of days
between $t_{u,v}$ and $t_{1}$.

Our graph neural network is composed of $3$ graph convolution blocks, where
each block takes the output of the previous block as its input except for the
first block. We concatenate the engineered vertex features with a trainable
vertex embedding matrix as the input to the first block of the graph neural
network. Each graph convolution block consists of a SAGEConv layer [20], a
Batch Normalization layer [21], a ReLU activation layer, and a skip connection
[22]:

$h^{l}=\text{ReLU}(\text{BatchNorm}(\text{SAGEConv}(h^{l-1},\hat{A})))+h^{l-1},$
(4)

TABLE II: Validation performance for LightGBM, the time-aware graph neural network, and the blend. Approach | AUC | Training time
---|---|---
LightGBM |  0.90293041 | 5 minutes
GNN | 0.90297814 | 97 minutes
Blend | 0.90314956 | 102 minutes

where $h^{l}$ denotes the vertex representations learned by the $l$-th graph
convolution block.

After obtaining the learned vertex representations from the output of the last
graph convolution block, we index the representations corresponding to
vertices $u$ and $v$ and concatenate them with the engineered pairwise
features. The concatenated features are then fed into a $5$-layer MLP with
densely connected layers [23]. We use ReLU activation units for intermediate
layers of the MLP, and use the sigmoid activation unit in the final layer,
which squashes the output value to the interval of $(0,1)$ to represent the
probability of existence of the new edge.

Deep neural networks are very sensitive to the scale of the inputs, and hence
the need for standardization. We preprocess the engineered features before
feeding them into the neural network by first filling all the missing values
with $0$ and then standardizing each feature through dividing the mean
subtracted values by the standard deviation. The training objective for the
neural network is the same as that for LightGBM, which is to minimize the
binary cross entropy loss. The model is optimized with the Adam optimizer [24]
with batch size $512$, learning rate of $1e-4$, and a weight decay of $1e-3$.
We apply a dropout of $0.2$ to prevent the model from overfitting and train
the model for $30$ epochs.

###  II-F Blending

Our final solution is based on a blend of the gradient boosting approach and
the deep learning approach. We use blending instead of other ensembling
techniques such as stacking since we observe that our models have very high
variances and relatively low biases. For blending, we begin by taking the
third power of the link prediction scores given by each model so as to amplify
the relative difference [25] and then compute the weighted average of the
transformed scores as the final link prediction score for each vertex pair.

##  III Experiments

###  III-A Dataset and Evaluation Metric

The dataset used in the competition is a dynamic semantic network
characterizing the content of the scientific literature in the field of
artificial intelligence between 1994 and 2017. The network consists of
$64,719$ vertices, each representing an artificial intelligence concept. Edges
between vertices are drawn when two concepts are investigated together in a
scientific paper, and the timestamp of the edge records the date of
publication of the scientific paper.

The evaluation metric used in the competition is AUC (area under the
curve–receiving operating characteristic), which is a widely adopted metric to
evaluate the performance of a binary classifier on different confidence
thresholds. ROC is a probability curve that visualizes change of the true
positive rate with respect to the false positive rate on different threshold
values. The AUC represents probability that a positive sample is ranked higher
than a negative sample [26, 27]. An AUC equal to 1 is ideal and represents
that the model can perfectly distinguish positive samples from negative
samples.

###  III-B Results and Analysis

TABLE III: Ablation results for the graph neural network. Approach | # SAGEConv | AUC
---|---|---
MLP | - | 0.89503802
MLP w/ vertex embedding | - | 0.89712115
GNN w/ binary adjacency matrix | 3 | 0.90203951
GNN w/ exponential decay | 1 | 0.90217735
GNN w/ exponential decay | 2 | 0.90254369
GNN w/ exponential decay | 3 | 0.90297814
GNN w/ exponential decay | 4 | 0.90278343
GNN w/ exponential decay | 5 | 0.90273941

The validation performance of LightGBM, the time-aware graph neural network,
and the blend is presented in Table II. We can see that the graph neural
network is able to achieve slightly better performance, while LightGBM is
significantly more efficient. The blend is able to further boost the
performance, albeit by a small margin.

We additionally performed an ablation study on the graph neural network; see
Figure III. We observe that directly applying the simple MLP on top of our
engineered features fails to achieve comparable performance with LightGBM and
the graph neural network. Equipping the MLP with a vertex embedding module
improves its performance since it can now learn vertex representations in a
data-driven manner. The vanilla graph neural network, which performs graph
convolution operations over the binary adjacency matrix, has a significant
performance advantage over the MLP since it is able to learn vertex
representations aggregated within a certain neighborhood region. Our
exponential decayed adjacency matrix further improves the performance by
enabling the model to learn time-aware vertex representations. We also studied
the impact of the number of graph convolution blocks on the performance of the
graph neural network. Empirically we found that the graph neural network
achieves the best performance when the number of graph convolution blocks is
set to $3$. This may indicate that the neighborhood information within 3-hop
is most useful in the semantic network of this competition.

The final leaderboard result
222https://www.iarai.ac.at/science4cast/competitions/2021-core/?leaderboard
for the top $5$ teams is presented in Table IV. Our solution achieved the
$1$st place among all the participating teams, with more than $1.3\%$ of
relative improvement over the $2$nd place team.

##  IV Discussion and Future Work

In this section we describe a few possible avenues of future improvements to
our solution.

When creating the training set, we fix $t_{1}$ to 20111231 and $t_{2}$ to
20141231. This limits the number of positive samples in the training set and
makes the model prone to suffer from overfitting. A possible way to obtain
more positive samples in the training set is to vary $t_{1}$ and $t_{2}$. We
can obtain a number of training sets with different $t_{1}$ and $t_{2}$, and
combine them into a larger training set with more positive samples. This can
potentially lead to further improvement to our solution since the performance
of machine learning models is largely dependent on the size of the training
dataset.

In addition, we demonstrate that using a weighted adjacency matrix can lead to
better performance of the graph neural network. In our solution, the weights
for the edges are defined using a handcrafted formula. An alternative approach
is to learn the edge weights with deep neural networks. For example, we can
employ the graph attention network (GAT) [28] to dynamically learn the edge
weights through an attention mechanism. Given sufficient amounts of training
data, we typically find handcrafted features to be outperformed by the
representations learned with deep neural networks.

Furthermore, in our graph neural network approach, the embedding vectors for
the vertices are fixed through time. This limits the capability of the vertex
embedding module to capture the evolution of the vertex over time.
Alternatively, we can try incorporating a time dependent embedding module into
our graph neural network [29].

TABLE IV: The final leaderboard result for the top $5$ teams. Rank | Team | AUC
---|---|---
1 | oahciy | 0.9283886
2 | Hash Brown | 0.9273865
3 | SanatisFinests2 | 0.9221236
4 | Bacalhink | 0.9185331
5 | nimasanjabi | 0.9184558

##  V Conclusion

We present our solution to the Science4cast 20201 competition, which is based
on a blend of a tree-based gradient boosting approach and a deep learning
approach. Our approach yields highly competitive results, achieving the $1$st
place out of all the participating teams.

## References

  * [1] Krenn, Mario, and Anton Zeilinger. ”Predicting research trends with semantic and neural networks with an application in quantum physics.” Proceedings of the National Academy of Sciences 117, no. 4 (2020): 1910-1916.
  * [2] Fortunato, Santo, Carl T. Bergstrom, Katy Börner, James A. Evans, Dirk Helbing, Staša Milojević, Alexander M. Petersen et al. ”Science of science.” Science 359, no. 6379 (2018).
  * [3] Lu, Yichao, Cheng Chang, Himanshu Rai, Guangwei Yu, and Maksims Volkovs. ”Learning Effective Visual Relationship Detector on 1 GPU.” arXiv preprint arXiv:1912.06185 (2019).
  * [4] Volkovs, Maksims, Himanshu Rai, Zhaoyue Cheng, Ga Wu, Yichao Lu, and Scott Sanner. ”Two-stage model for automatic playlist continuation at scale.” In Proceedings of the ACM Recommender Systems Challenge 2018, pp. 1-6. 2018.
  * [5] Volkovs, Maksims, Anson Wong, Zhaoyue Cheng, Felipe Pérez, Ilya Stanevich, and Yichao Lu. ”Robust contextual models for in-session personalization.” In Proceedings of the Workshop on ACM Recommender Systems Challenge, pp. 1-5. 2019.
  * [6] Lu, Yichao, Cheng Chang, Himanshu Rai, Guangwei Yu, and Maksims Volkovs. ”Multi-view scene graph generation in videos.” In International Challenge on Activity Recognition (ActivityNet) CVPR 2021 Workshop, vol. 3. 2021.
  * [7] Yu, Wei, Yichao Lu, Steve Easterbrook, and Sanja Fidler. ”Efficient and Information-Preserving Future Frame Prediction and Beyond.” In International Conference on Learning Representations. 2019.
  * [8] Yu, Wei, Yichao Lu, Steve Easterbrook, and Sanja Fidler. ”CrevNet: Conditionally Reversible Video Prediction.” arXiv preprint arXiv:1910.11577 (2019).
  * [9] Lu, Yichao, Himanshu Rai, Jason Chang, Boris Knyazev, Guangwei Yu, Shashank Shekhar, Graham W. Taylor, and Maksims Volkovs. ”Context-aware Scene Graph Generation with Seq2Seq Transformers.” In Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 15931-15941. 2021.
  * [10] Lu, Yichao, Ruihai Dong, and Barry Smyth. ”Coevolutionary recommendation model: Mutual learning between ratings and reviews.” In Proceedings of the 2018 World Wide Web Conference, pp. 773-782. 2018.
  * [11] Lu, Yichao, Ruihai Dong, and Barry Smyth. ”Why I like it: multi-task learning for recommendation and explanation.” In Proceedings of the 12th ACM Conference on Recommender Systems, pp. 4-12. 2018.
  * [12] Jannach, Dietmar, Gabriel de Souza P. Moreira, and Even Oldridge. ”Why are deep learning models not consistently winning recommender systems competitions yet? A position paper.” In Proceedings of the Recommender Systems Challenge 2020, pp. 44-49. 2020.
  * [13] Page, Lawrence, Sergey Brin, Rajeev Motwani, and Terry Winograd. The PageRank citation ranking: Bringing order to the web. Stanford InfoLab, 1999.
  * [14] Ahmad, Iftikhar, Muhammad Usman Akhtar, Salma Noor, and Ambreen Shahnaz. ”Missing link prediction using common neighbor and centrality based parameterized algorithm.” Scientific reports 10, no. 1 (2020): 1-9.
  * [15] Ke, Guolin, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. ”Lightgbm: A highly efficient gradient boosting decision tree.” Advances in neural information processing systems 30 (2017): 3146-3154.
  * [16] Chen, Tianqi, Tong He, Michael Benesty, Vadim Khotilovich, Yuan Tang, and Hyunsu Cho. ”Xgboost: extreme gradient boosting.” R package version 0.4-2 1, no. 4 (2015): 1-4.
  * [17] Dorogush, Anna Veronika, Vasily Ershov, and Andrey Gulin. ”CatBoost: gradient boosting with categorical features support.” arXiv preprint arXiv:1810.11363 (2018).
  * [18] Burges, Christopher JC. ”From ranknet to lambdarank to lambdamart: An overview.” Learning 11, no. 23-581 (2010): 81.
  * [19] Cao, Zhe, Tao Qin, Tie-Yan Liu, Ming-Feng Tsai, and Hang Li. ”Learning to rank: from pairwise approach to listwise approach.” In Proceedings of the 24th international conference on Machine learning, pp. 129-136. 2007.
  * [20] Hamilton, William L., Rex Ying, and Jure Leskovec. ”Inductive representation learning on large graphs.” In Proceedings of the 31st International Conference on Neural Information Processing Systems, pp. 1025-1035. 2017.
  * [21] Ioffe, Sergey, and Christian Szegedy. ”Batch normalization: Accelerating deep network training by reducing internal covariate shift.” In International conference on machine learning, pp. 448-456. PMLR, 2015.
  * [22] He, Kaiming, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. ”Deep residual learning for image recognition.” In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770-778. 2016.
  * [23] Huang, Gao, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q. Weinberger. ”Densely connected convolutional networks.” In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 4700-4708. 2017.
  * [24] Kingma, Diederik P., and Jimmy Ba. ”Adam: A method for stochastic optimization.” arXiv preprint arXiv:1412.6980 (2014).
  * [25] Lu, Yichao, Ruihai Dong, and Barry Smyth. ”Context-aware sentiment detection from ratings.” In International Conference on Innovative Techniques and Applications of Artificial Intelligence, pp. 87-101. Springer, Cham, 2016.
  * [26] Yuan, Zhuoning, Yan Yan, Milan Sonka, and Tianbao Yang. ”Large-scale Robust Deep AUC Maximization: A New Surrogate Loss and Empirical Studies on Medical Image Classification.” In Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 3040-3049. 2021.
  * [27] Ataman, Kaan, W. Nick Street, and Yi Zhang. ”Learning to rank by maximizing AUC with linear programming.” In The 2006 IEEE International Joint Conference on Neural Network Proceedings, pp. 123-129. IEEE, 2006.
  * [28] Veličković, Petar, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. ”Graph Attention Networks.” In International Conference on Learning Representations. 2018.
  * [29] Rossi, Emanuele, Ben Chamberlain, Fabrizio Frasca, Davide Eynard, Federico Monti, and Michael Bronstein. ”Temporal graph networks for deep learning on dynamic graphs.” arXiv preprint arXiv:2006.10637 (2020).


Burst detection methods
Ellese Cotterill, Stephen J. Eglen
Department of Applied Mathematics and Theoretical Physics
University of Cambridge
CB3 0WA UK

Abstract
‘Bursting’, defined as periods of high frequency firing of a neuron separated
by periods of quiescence, has been observed in various neuronal systems, both
in vitro and in vivo. It has been associated with a range of neuronal
processes, including efficient information transfer and the formation of
functional networks during development, and has been shown to be sensitive to
genetic and pharmacological manipulations. Accurate detection of periods of
bursting activity is thus an important aspect of characterising both
spontaneous and evoked neuronal network activity. A wide variety of
computational methods have been developed to detect periods of bursting in
spike trains recorded from neuronal networks. In this chapter, we review
several of the most popular and successful of these methods.

### Table of Abbreviations

CMA | Cumulative Moving Average
---|---
IQR | Inter-Quartile Range
IRT | ISI Rank Threshold
ISI | InterSpike Interval
LTD | Long Term Depression
LTP | Long Term Potentiation
MEA | MultiElectrode Array
MI | Max Interval
PS | Poisson Surprise
RS | Rank Surprise
RGS | Robust Gaussian Surprise

##  1 Introduction

Neuronal bursting, observed as intermittent periods of elevated spiking rate
of a neuron (see Figure 1), has been observed extensively in both in vitro and
in vivo neuronal networks across various network types and species [98, 19,
65]. ††margin:  These bursts can be isolated to a single neuron or, commonly,
occur simultaneously across many neurons, in the form of ‘network bursts’ [93,
95, 66, 4]

Figure 1: Example of bursting activity in a spike train recorded from mouse
retinal ganglion cells. Horizontal blue lines show the location of bursts.
Scale bar represents 1 s.

Bursting activity is believed to play a role in a range of physiological
processes, including synapse formation [52] and long-term potentiation [48].
Analysis of patterns of bursting activity can thus be used as a proxy for
studying the underlying physiological processes and structural features of
neuronal networks. A common method of studying bursting activity in vitro
involves the use of MEA recordings of spontaneous or evoked neuronal network
activity [50, 14, 69, 93]. This approach has been employed to study changes in
spontaneous network activity over development [95], and the effect of
pharmacological or genetic manipulations [23, 15].

Despite the importance of bursting and its prevalence as a feature used to
analyse neuronal network activity, there remains a lack of agreement in the
field about the definitive formal definition of a burst [21, 30]. There is
also no single technique that has been widely adopted for identifying the
location of bursts in spike trains. Instead, a large variety of burst
detection methods have been proposed, many of which have been developed and
assessed using specific data sets and single experimental conditions. As most
studies of bursting activity have been performed on experimental data from
recordings of rodent ††margin:  neuronal networks [14, 55], this type of data
has most often been used to assess the performance of burst detection
techniques [19, 55, 30].

Recently, it has been shown that networks of neurons derived from human stem
cells can be grown successfully on MEAs and exhibit spontaneous electrical
activity, including bursting [36, 33]. Human stem-cell derived neuronal
cultures have also been demonstrated to be a suitable alternative to rodent
neuronal networks in applications such as neurotoxicity testing [100]. This
has ††margin:  led to a demand for a robust method of analysing bursting in
these networks, which commonly exhibit more variable and complex patterns of
bursting activity than rodent neuronal networks [40] (see Figure 2). Recently,
some burst detection methods have been developed which specifically focus on
analysing bursting activity in these types of variable networks [40, 88].

Figure 2: Examples of spike trains from mouse and human neuronal networks.
Each row represents the spikes recorded from one electrode and the scale bar
represents 30 s. Recordings from human neuronal networks often exhibit more
variable and complex spontaneous activity patterns. ††margin:

##  2 Physiological significance of neuronal bursting

Neuronal bursting is a frequently observed phenomenon in MEA recordings of
cultures of dissociated neurons, as well as in numerous in vitro systems [95,
66, 98, 46]. In cultured rodent cortical networks, bursts, and in particular
synchronised ‘network bursts’ generally arise as a feature of the spontaneous
network activity after around 1 week in vitro [38]. Most studies observe that
these network bursts then increase in frequency and size before reaching a
peak around 3 weeks in vitro [91, 93, 18]. This peak in network bursting
activity generally corresponds to the period in which the synaptic density of
the network reaches its maximum [90, 38, 91]. This is followed by a period of
shortening of network bursts, which coincides with a stage of ’pruning’ or
reduction in dendritic spine synapses and maturation of excitatory connections
between neurons [18, 36, 35, 92]. As well as being correlated with neuronal
network development and maturation, bursting patterns of spontaneous activity
are also believed to play an important role in regulating cell survival. High
frequency bursting has been shown to increase neuronal survival in cortical
cultures, while suppression of spontaneous activity has been observed to
greatly increase rates of programmed cell death [29, 32].

Bursting has also been observed to be involved with a range of physiological
processes in mature neuronal networks. For example, bursting is believed to be
a more efficient method of information transfer between neurons than single
spikes. Central synapses in various brain regions have been shown to exhibit
low probabilities of neurotransmitter release in response to single
presynaptic spikes, making information transfer by single spikes unreliable
[8, 10, 2]. However, bursts of spikes can lead to ‘facilitation’, a process in
which a rapid succession of spikes leads to a build up of intracellular Ca2+
in the presynaptic terminal. This increases the probability of
neurotransmitter release and resultant production of EPSPs with subsequent
spikes [83, 43]. In addition to being involved in these mechanisms of short
term plasticity, bursting has also been implicated in long term potentiation
(LTP) and depression (LTD). For example, in the hippocampus, postsynaptic
bursting at temporally relevant intervals could produce long term synaptic
changes [68, 26, 82].

It has also been suggested that bursts of spikes transmit information with a
higher signal-to-noise ratio than single spikes [78]. Evidence of this has
been seen in a variety of brain regions, such as the hippocampus, where place-
fields have been shown to be more accurately defined by bursts than individual
spikes [64]. Bursting has also been shown to produce sharper sensory tuning
curves [13, 43] and more reliable feature extraction than single spikes [27,
78, 44].

The importance of neuronal bursting has also been demonstrated through its
association with a variety of behaviours in vivo, including visual processing,
reward and goal directed behaviour and sleep and resting conditions [13, 43,
84, 75, 74, 25, 6, 56, 98, 80]. Bursting of hippocampal place cells has also
been observed during exploration of new environments [63, 24]. The presence of
bursting in these, as well as other memory-related behaviours [11, 99],
suggests that bursting plays a specific role in memory and learning in the
adult brain [67].

Additionally, bursting activity has been seen to be altered in certain
pathological conditions [96, 37, 58, 79]. For example, increased bursting
activity has been observed in the basal ganglia of Parkinson’s patients, with
correlations between the level of bursting activity and the progression of the
disease [49, 60]. This suggests that the study of bursting activity could not
only reveal important features of normal brain function, but also how this is
altered in diseased states.

††margin:  ††margin:  *deleted this section*

##  3 Previous approaches to burst detection

Since the development of the first methods to identify bursting in neuronal
networks more than three decades ago, many techniques have been proposed.
These methods take a variety of approaches. ††margin:

###  3.1 Fixed threshold-based methods

The simplest approaches involve imposing thresholds on values such as the
minimum firing rate or maximum allowed interspike interval (ISI) within a
burst, and classifying any sequence of consecutive spikes satisfying these
thresholds as a burst. In well-ordered spike trains, these thresholds can be
set as fixed values by visual inspection [98, 19]. Other methods also
incorporate additional thresholds on relevant parameters such as the minimum
interval between two bursts and the minimum duration of a burst, to restrict
detected bursts to those with biologically realistic properties [59].
††margin:

###  3.2 Adaptive threshold-based methods

As opposed to having fixed threshold parameters that are chosen by the user,
other burst detection algorithms derive the values of their threshold
parameters adaptively from properties of the data, such as the mean ISI [16]
or total spiking rate [69]. Commonly, this involves the use of some form of
the distribution of ISIs on a spike train. For spike trains containing
bursting activity, the smoothed histogram of ISIs on the train should have a
peak in the region of short ISIs, which represents within-burst ISIs, and one
or more peaks at higher ISI values, representing intra-burst intervals. A
threshold for the maximum ISI allowed within a burst can be set at the ISI
value representing the turning point in the histogram [21].

Several other adaptive burst detection algorithms also use distributions
related to the ISI histogram to calculate thresholds for burst detection. [76]
and [65] argue that the histogram of log(ISI)s provides a better separation of
within and between-burst intervals, and use this histogram to set the
threshold for the maximum within-burst ISI at the minimum between the first
two well separated peaks. [39] use the histogram of discharge density rather
than ISIs for burst detection, ††margin:  while [40] derive the threshold
parameters for detecting bursts in their algorithm from the cumulative moving
average of the ISI histogram. ††margin:

###  3.3 Surprised-based methods

Another category of burst detection techniques are the surprise-based methods,
which use statistical techniques to distinguish periods of bursting from
baseline neuronal firing. The earliest of such methods was developed by [46],
and detects bursts as periods of deviation from an assumed underlying Poisson
process of neuronal firing. This method critically assumes Poisson-distributed
spike trains, which has been shown to be inappropriate for many common spike
trains, in particular because of the refractory period between spikes [12].
Despite this, the Poisson Surprise method has been one of the most widely used
burst detection methods since its development over thirty years ago ††margin:
(398 citations as of June 2018) and is still commonly used for analysing
bursting activity in experimental studies of numerous neuronal network types
[79, 70, 77]. More recently, other surprise-based burst detectors have been
developed that replace the assumption that baseline firing follows a Poisson
process with other assumptions about the underlying distribution of spikes
[42, 30]. ††margin:

###  3.4 Other methods

Other burst detectors take alternative approaches to separate bursting from
background spiking activity. [86] examine the slope of the plot of spike time
against spike number to detect bursts as periods of high instantaneous slope.
[54] require bursts to be separated by intervals at least two standard
deviations greater than their average within-burst ISIs, while [81] propose a
parameter-free burst detection method, in which sequences of spikes are
classified as bursts if the sum of their within-bursts ISIs is less than the
ISIs immediately before and after the burst.

Numerous studies have also used various forms of Hidden Markov Models to
analyse neuronal activity patterns [71, 17, 1]. These methods assume that a
neuron stochastically alternates between two or more states, characterised by
differences in their levels of activity. [85] apply this idea to burst
detection by modelling neuronal activity using hidden semi-Markov models.

###  3.5 Burst detection methods

In this section, we will outline a number of key existing burst detection
algorithms. Given the vast number of available burst detection techniques, the
following have been chosen for their relevance and popularity in the existing
literature ††margin:  , and represent examples of each of the approaches to
burst detection outlined above.

Abbrev. |  Method | Reference
---|---|---
Fixed threshold-based methods
MI |  MaxInterval | [59]
Adaptive threshold-based methods
logISI |  LogISI | [65]
CMA |  Cumulative Moving Average | [40]
IRT |  ISI Rank Threshold | [34]
Surprise-based methods
PS |  Poisson Surprise | [46]
RS |  Rank Surprise | [30]
RGS |  Robust Gaussian Surprise | [42]
Other methods
HSMM |  Hidden Semi-Markov Model | [85]
Table 1:  Burst detectors classified by their approach to burst detection.

MaxInterval method [59]
Bursts are defined using five fixed threshold parameters, shown in Figure 3.
The value of these parameters are chosen a priori and any series of spikes
that satisfy these thresholds is classified as a burst.

Figure 3: Illustration of the parameters used by the MaxInterval method.

LogISI method [65]
The histogram of log(ISI)s on a spike train is computed, using a bin size of
0.1 in $\log(\text{ISI})$ units. Let $C_{k}$ denote the ISI count in the $k$th
bin of this histogram, which corresponds to an ISI size of ISIk, and MCV
denote a pre-specified threshold value, known as the maximum cutoff value. The
location of the peaks of this histogram are found using a custom peak finding
algorithm described in [65]. The largest peak of the histogram corresponding
to an ISI less than or equal to MCV is set as the intraburst peak, $C_{IBP}$.
If no peak is found in the histogram with ISI${}_{k}\leq$ MCV, the spike train
is classified as containing no bursts.

In the case that an intraburst peak is present, the minimum value of the
histogram between the intraburst peak and each of the following peaks,
$C_{p_{i}}$ ($i=1,...,N$), is found. For each minimum, a void parameter is
calculated that represents how well the corresponding peak is separated from
the intraburst peak, as

$void(i)=1-\frac{C_{min_{i}}}{\sqrt{C_{IBP}\cdot C_{p_{i}}}}$

where $C_{min_{i}}$ is the minimum value of $C_{k}$ for $IBP<k<p_{i}$.

The smallest ISI${}_{min_{i}}$ for which $void(i)>0.7$ is set as the threshold
for the maximum ISI in a burst, $maxISI$ (see Figure 4). Any series of at
least three spikes separated by ISIs less than $maxISI$ are classified as
bursts. If no point with a void value above 0.7 is found, or if $maxISI>$ MCV,
bursts are detected using MCV as the threshold for the maximum ISI in a burst
and then extended to include spikes within $maxISI$ of the beginning or end of
each of these bursts.

Figure 4: Example of log-adjusted ISI histogram with the threshold for
intraburst and interburst intervals found using the logISI method.

Cumulative Moving Average (CMA) method [40]
This method also uses the histogram of ISIs on a spike train. The cumulative
moving average (CMA) at each ISI bin of the histogram is calculated. The CMA
of the $N$th ISI bin is defined as

$CMA_{N}=\frac{1}{N}\sum^{N}_{k=1}C_{k}\,,$

where $C_{k}$ is the ISI count in the $k$th bin. The skewness of the CMA
distribution is used to determine the values of two threshold parameters,
$\alpha_{1}$ and $\alpha_{2}$, based on the scale given in [40]. The maximum
of the CMA distribution, $CMA_{max}$, is found and the value of $maxISI$ is
set at the ISI bin at which the CMA is closest in value to $\alpha_{1}\cdot
CMA_{max}$ (see Figure 5). Burst cores are then found as any sequences of at
least three spikes separated by ISIs less than $maxISI$.

[40] suggest extending these burst cores to include burst-related spikes.
These are found using a second cutoff, set at the value of the ISI bin at
which the CMA is closest to $\alpha_{2}\cdot CMA_{max}$. Spikes within this
cutoff distance from the beginning or end of the existing burst cores are
classified as burst-related spikes. For this study, only the burst cores
detected by this method were examined, omitting any burst related spikes.

Figure 5: Example of ISI histogram with the threshold for intraburst and
interburst intervals found using the CMA method. Red line shows the cumulative
moving average of the ISI histogram.

ISI rank threshold method [34]
In the ISI rank threshold (IRT) method, the rank of each ISI on a spike train
relative to the largest ISI on the train is calculated, with $R(t)$ denoting
the rank of the ISI beginning at time $t$. The probability distribution,
$P(C)$, of spike counts in one-second time bins over the spike train is also
found. A rank threshold, $\theta_{R}$, is set to a fixed value, and a spike
count threshold, $\theta_{C}$, is calculated from $P(C)$. A burst is then
defined to begin at a spike at time $t$ if the rank of the proceeding ISI
satisfies $R(t)<\theta_{R}$ and the spike count in the following second,
$C(t,t+1)$, exceeds $\theta_{C}$. The burst continues until a spike is found
for which $C(t,t+1)<\frac{\theta_{C}}{2}$.

Poisson Surprise method [46]
The average firing rate, $\lambda$, on a spike train is calculated, and the
underlying activity on this spike train is assumed to follow a Poisson process
with rate $\lambda$. The Poisson Surprise (PS) statistic for any period of
length $T$ containing $N$ spikes is calculated as

$S=-\log P$

where

$P=\exp\left(-\lambda T\sum_{n=N}^{\infty}\frac{(\lambda T)^{n}}{n!}\right)$

is the probability that $N$ or more spikes occur randomly in a period of
length $T$.

A surprise maximization algorithm described in [46] is then used to find the
set of bursts that maximises the PS statistic across the entire spike train.
This involves initially identifying bursts as any sequence of three
consecutive spikes separated by ISIs which are less than half of the mean ISI
on the spike train. Spikes are then added to the end and removed from the
beginning of each of these initial bursts until the sequence of spikes with
the maximum PS statistic is found. Finally, any bursts which have a PS
statistic below a pre-defined threshold level are discarded.

Rank Surprise method [30]
The Rank Surprise (RS) burst detection algorithm is a non-parametric
adaptation of the Poisson Surprise approach. To implement this method, all
ISIs on a spike train are ranked by size, with the smallest ISI given a rank
of one. In the absence of any bursting activity, the ISI ranks should by
independently and uniformly distributed. For any period containing $N$ spikes
separated by $N-1$ ISIs with ranks $r_{n}$, …, $r_{n+N-1}$, the Rank Surprise
statistic is defined as

$RS=-\log(P(D_{N}\leq r_{n}+...+r_{n+N-1}))$

where $D_{N}$ is the discrete uniform sum distribution between 1 and $N$ and
$r_{n}$ is the rank of the $n$th ISI on the spike train.

Bursts are then chosen to maximise the RS statistic across the entire spike
train using an exhaustive surprise maximisation algorithm, outlined in [30]. A
fixed threshold for $maxISI$ is first calculated from the distribution of ISIs
on the spike train. The first sequence of at least three spikes with ISIs less
than $maxISI$ are found and an exhaustive search of all of the subsequences of
ISIs within this period is performed to find the subsequence with the highest
RS value. If this value is above a fixed minimum significance threshold,
chosen a priori, it is labelled as a burst. This process is repeated on the
remaining ISI subsequences within the period of interest until all significant
bursts are found. Following this, the next sequence of spikes with ISIs below
$maxISI$ is examined in a similar fashion, and this process is continued until
the end of the spike train.

Robust Gaussian Surprise method [42]
In the Robust Gaussian Surprise (RGS) method, the distribution of log(ISI)s on
each spike train is found and centred around zero. The normalised log(ISI)s
from each spike train in the study are then pooled and the central
distribution of this joint data set is found using a procedure outlined in
[42]. A burst detection threshold for $maxISI$ is set at the 0.5 percentile of
this central distribution, which is estimated as 2.58 times the median
absolute deviation of the distribution.

The Gaussian burst surprise value in any interval on a spike train is defined
as

$GS_{B}=-\log(P)$

where $P$ is the probability that the sum of normalised log(ISI)s in the
interval is greater than or equal to the sum of an equal number of i.i.d.
Gaussian random variables with mean and variance equal to that of the central
distribution.

Any consecutive sequence of spikes separated by intervals less than $maxISI$
are classified as burst cores. These burst cores are then extended by adding
intervals to the beginning and end of the burst cores until the sequence with
the maximum value of $GS_{B}$ is found. In the case of overlapping bursts, the
burst with the largest $GS_{B}$ value is retained. Finally, any detected
bursts with $GS_{B}$ below a predefined threshold value are discarded. [42]
also propose a similar method for identifying pauses in spike trains.

Hidden Semi-Markov Model method [85]
This method is based on the assumption that neurons switch stochastically
between two states: ‘non-bursting’ (state 0) and ‘bursting’ (state 1), that
can be modelled using a Hidden semi-Markov model. The transition times between
the two states are modelled using two Gamma distributions, $f_{0}^{ITI}$ and
$f_{1}^{ITI}$. Within each of the states, the ISI times are modelled using two
additional gamma distributions, $f_{0}^{ISI}$ and $f_{1}^{ISI}$. The
parameters of these four distributions are learned from the data. A custom
Markov chain Monte-Carlo algorithm described in [85] is then used to compute
the posterior probability that a neuron is in a bursting state at any given
time. A fixed threshold value is chosen a priori, and any periods during which
the posterior probability exceeds this value are classified as bursts.

###  3.6 Evaluation of burst detection techniques

In [22], we performed a thorough evaluation of the burst detection methods
outlined above. This involved first assessing the methods against a list of
desirable properties that we deemed an ideal burst detector should possess
††margin:  (see Table 2). This was achieved by generating synthetic spike
trains with specific properties of interest to represent each desirable
property. The output of each burst detector when used to analyse each set of
spike trains was then compared to the ‘ground truth’ bursting activity.

Desirable properties
---
D1 |  Deterministic: The method should detect the same bursts over repeated runs on the same data, to ensure consistency and reproducibility of results
D2 |  No assumption of spike train distribution: The method should not assume ISIs follow a standard statistical distribution, to ensure wide applicability to a variety of spike trains
D3 |  Number of parameters: The method should have few parameters, to reduce the variability inherently introduced through parameter choice
D4 |  Computational time: The method should run in a reasonable amount of time using standard personal computers
D5 |  Non-bursting trains: The method should detect few spikes as being within bursts in spike trains containing no obvious bursting behaviour
D6 |  Non-stationary trains: The method should detect few spikes as being within bursts in spike trains with non-stationary firing rates that contain no obvious bursting behaviour
D7 |  Regular short bursts: The method should detect a high proportion of spikes in bursts in spike trains containing short well-separated bursts
D8 |  Non-stationary bursts: The method should detect a high proportion of spikes in bursts in spike trains containing bursts with variable durations and numbers of spikes per burst
D8 |  Regular long bursts: The method should detect a high proportion of spikes in bursts and accurate number of bursts in spike trains containing long bursts with low within-burst firing rates
D10 |  High frequency bursts: The method should detect a high proportion of spikes in bursts and accurate number of bursts in spike trains containing a large number of short bursts
D11 |  Noisy train: The method should classify a high number of within-burst spikes as bursting and a low number of interburst spikes as bursting in spike trains containing both bursts and noise spikes
Table 2:  ††margin:  Desirable properties for a burst detector. Table
reproduced from [22].  ††margin:

Figure 6 shows the performance of the chosen burst detectors on a sample of
these properties. Most burst detectors can accurately detect a small amount of
bursting activity in spike trains simulated to contain no bursting behaviour
(Figure 6A), with the exception of the HSMM and CMA methods, which detect a
significant amount of erroneous bursting. Conversely, most burst detectors
accurately identified most bursting activity in spike trains containing only
regular short bursts (Figure 6C). However, the RS, IRT and RGS methods
performed poorly here, only detecting a small proportion of the bursting
activity.

Figure 6: Fraction of spikes in bursts found by each burst detector in 100
synthetic trains with A no bursting (D5), B no bursting and non-stationary
firing rate (D6), C short regular bursts (D7), D bursts with non-stationary
burst lengths and durations (D8). Dotted line shows desired result from an
ideal burst detector; methods close to this line are deemed to work well. In
each ‘box-and-whisker’ plot, boxes show the median $\pm$ inter-quartile range
(IQR), and whiskers extend to median $\pm\ 1.5\times$ IQR. Outliers are
represented as points. Figure reproduced with permission from [22].

This approach ††margin:  of assessing the performance of each burst detection
method against desirable properties allowed us to determine a ranking for each
of the burst detectors, in which the Rank Surprise, Robust Gaussian Surprise
and ISI Rank Threshold methods ranked particularly poorly ††margin:  (see
Table 3) .

Desirable property | Burst detection method
---|---
| PS | MI | CMA | RS | IRT | RGS | logISI | HSMM
D1 Deterministic | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | $\times$
D2 Distribution assumption | $\times$ | ✓ | ✓ | ✓ | ✓ | $\times$ | ✓ | $\times$
D3 Number of parameters | ✓ | $\times$ | ✓ | ✓ | ✓ | ✓ | ✓ | $\times$
D4 Computational time | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | $\times$
D5 Non-bursting | 4 | 1 | 7 | 5 | 6 | 3 | 1 | 8
D6 Non-stationary | 6 | 2 | 7 | 4 | 5 | 3 | 1 | 8
D7 Regular bursting | 4 | 1 | 2 | 7 | 6 | 7 | 5 | 3
D8 Non-stat bursts | 4 | 3 | 5 | 7 | 6 | 8 | 2 | 1
D9 Long bursts | 2 | 4 | 3 | 8 | 5 | 7 | 6 | 1
D10 High frequency | 5 | 1 | 4 | 7 | 6 | 8 | 2 | 3
D11 Noisy bursts | 5 | 1 | 2 | 7 | 6 | 8 | 4 | 2
Total (Relative rank) | 30 (4) | 13 (1) | 30 (4) | 45 (8) | 40 (6) | 44 (7) | 21 (2) | 26 (3)
Table 3:  ††margin:  The performance of each method on the desirable
properties specified in Table 2. For binary properties, D1–-D4, each method
was judged to either possess the property or not, while for properties
D5–-D11, the performance of each method was ranked against the other methods
(1=best, 8=worst) and summed to produce an overall ranking. Table adapted from
[22].

Further assessment of the burst detectors was then achieved by examining the
coherence of the bursts detected by each method with visually annotated bursts
in experimental recordings of mouse retinal ganglion cells (RGCs). This
allowed us to analyse the specificity and sensitivity of the burst detectors
as their input parameters were varied. This analysis reinforced the low levels
of adaptability of the RS, RGS and IRT methods at analysing this type of data.
The HSMM method was also seen to have a consistently high false positive rate
compared to other burst detectors used to analyse this data.

Based on these assessments, four burst detectors, namely the MI, logISI, PS
and CMA methods, were chosen as the best performing burst detection methods,
and used to analyse bursting activity in novel recordings of networks of human
induced pluripotent stem cell (hiPSC)-derived neuronal networks over several
months of development. This analysis showed a slight increase in the
proportion of bursting activity observed in these networks as they mature,
although this increase was far lower than that which has been observed in
developing rodent neuronal networks [14, 19, 95].

From this analysis, we concluded that no existing burst detector possesses all
of the desirable properties required for ‘perfect’ identification of bursting
periods in highly variable networks. The CMA and PS methods possessed many of
the desirable properties, but had limitations such as their tendency to
overestimate bursting activity in spike trains containing sparse or no
bursting activity, particularly those with a non-stationary firing rate.

Overall, the MI and logISI methods showed the most promise for achieving
robust burst analysis in a range of contexts. These methods possessed most
properties we deemed desirable for a burst detection method and were generally
able to achieve high coherence with visually detected bursts in experimental
MEA recordings. These methods, however, still had limitations. The MI method
requires the choice of five parameters, the optimal values of which can be
challenging to determine, particularly when analysing recordings from a
variety of experimental conditions [22]. The logISI method had a tendency to
underestimate bursting in some spike trains, particularly those with non-
standard bursting activity.

The overall recommendation from this analysis was to choose a burst detector
from the several high performing methods outlined above based on the number of
freedom the user wishes to control. The MI method is a good first choice for
these purposes, and despite the large number of parameters this method
requires, these parameters are easy to interpret biologically and adjust to
achieve the desired burst detection results for the specific situations in
which it is utilised. If appropriate parameters cannot be found for the MI
method, a high performing alternative is the logISI method, which can be
implemented without choosing any input parameters. This method is most
effective when there is a clear distinction between the size of within-burst
and between-burst intervals on a spike train. In cases when this distinction
is not apparent, the PS and CMA methods are reasonably effective alternative
burst detection methods, however post hoc screening for outliers in terms of
burst duration is advisable when using either of these methods.

One robust approach to burst detection would be to use several burst detectors
to analyse the data of interest, and compare the results of each method. If
the burst detectors are largely in agreement, this provides confidence in the
nature of the bursting activity identified in the experimental data. Any major
discrepancies between the results from the methods can also be used to
identify areas where one or more burst detectors may be performing poorly,
which can be further investigated through inspection of the specific spike
trains of interest.

###  3.7 Network-wide burst detection

As well as single neuron bursts, synchronous bursting of networks of neurons,
termed ‘network bursts’, are a ubiquitous feature of various neuronal
networks. In rat cortical cultures, these network bursts have been observed to
arise from around one week in vitro, and comprise the dominant form of
spontaneous network activity at this age [19, 91]. Network bursts increase in
frequency and size before reaching a peak at around 3 weeks in vitro,
corresponding to the period in which synaptic density in the network reaches
is maximum [91, 93, 18].

As well as in rat cortical cultures, the presence of network bursting activity
has also been observed in a variety of other brain regions and species in
vitro [89, 7, 73, 31, 57] and in vivo [20, 47, 97]. Recently, synchronous
bursting resembling that in rat cortical cultures has also been observed in
networks produced from human embryonic or induced pluripotent stem-cell
derived neurons, generally arising 8 to 12 weeks after differentiation and
increasing in frequency over development [33, 62, 3].

####  3.7.1 Existing network burst detection techniques

A variety of techniques have been developed to detect these network-wide
bursts. Several of these methods identify bursts as increases in the network-
wide firing rate [55, 72]. These periods, however, do not necessarily consist
of single neurons bursts across multiple electrodes. Other methods define
network bursts only when single-neuron bursts occur simultaneously across
numerous recorded electrodes [95, 65]. ††margin:  For example, [4] combine the
spikes detected on all channels of an MEA into a single spike train and employ
the ISI histogram between every $n$th spike in this network-wide spike train
to determine an appropriate threshold for the maximum ISI within a network
burst. [94], on the other hand, detect ‘burstlets’ on each electrode
individually using an adaptive threshold based on the electrode’s average
firing rate. A network burst is then defined as any period in which burstlets
on multiple electrodes overlap.

Network-wide information can also be incorporated into single-neuron burst
detection techniques to improve their performance. [53] showed that the peaks
corresponding to intra and interburst spikes in an ISI histogram were better
separated when pooled ISIs from multiple electrodes of an MEA were included,
rather than simply those from a single spike train. They also proposed a pre-
processing technique designed to improve the detection of bursts, particularly
on noisy data. This involves creating a return map, which plots the ISI
immediately proceeding each spike ($ISI_{pre}$) against the ISI following the
spike ($ISI_{post}$). Background spikes lie in the region of this graph with
both high $ISI_{pre}$ and $ISI_{post}$, and are removed from consideration by
the burst detection method. The performance of various single channel burst
detection techniques were shown to be significantly improved when applied to
data pre-processed in this way, compared to the original data [53].

Additionally, [88] adapted the CMA method of [40] to incorporate information
from multiple MEA electrodes. In this multi-CMA method, instead of individual
histograms for each spike train, the ISI histogram from the combined ISIs from
multiple electrodes are used to calculate the threshold for burst detection in
an identical method to the original CMA method. This threshold is then used to
detect bursts on each electrode individually. The electrodes that are used for
combined analysis by this method can be chosen from a variety of options,
including analysing all electrodes in a single MEA simultaneously, analysing
the spike trains from a single electrode over several experimental time
points, or analysing all electrodes over all time points in the experiment.
This adaptation has been shown to reduce the number of excessively long sparse
bursts identified by the original CMA method, improving its performance at
analysing highly variable spike trains.

###  3.8 Summary and future directions

In this article we have summarised the main techniques of burst detection.
Moving from an informal definition (”bursts are groups of spikes that are
close to each other in time”) to a formal mathematical definition has proved
challenging. Our experience is that when the datasets are relatively clean,
there is good agreement between methods. However, when the data are noisy, not
only ††margin:  do different methods ††margin:  disagree , different human
observers will also disagree. Here we have outlined several of the methods
that we believe work relatively well, but are fallible when presented with
noisy data. Future work in this area might be centred around developing
methods that are more robust to noisy data. ††margin:  Possible steps towards
this may involve generating more realistic synthetic datasets to train and
assess burst detection techniques, or the incorporation of noise-reducing
preprocessing steps prior to burst detection, such as those developed by [53]

††margin:  Outside of neuroscience, the detection of ‘bursty’ events is also a
more general problem in time series analysis. For example, identifying bursts
of gamma rays can aid in the detection of black holes, and the detection of
periods of high trading volume of a stock is of relevance to regulators
looking for insider trading [102]. Various techniques have been developed for
detecting bursting periods in these and other data types, including sliding
window and infinite state automaton-based models [102, 101, 41, 9, 45]. Ideas
from these burst detectors developed in other domains may be useful for
informing future approaches to burst detection in a neuroscience context.

The increasing use of high density MEAs, which contain up to several thousand
electrodes [51, 50], to record in vitro neuronal activity as well as the
prevalence of multi-well MEAs in applications such as high-throughput
neurotoxicity screening [87, 61] and drug safety testing [28] also has
implications for burst detection. In particular, the computational complexity
of burst detection methods becomes increasingly relevant in these high-
throughput situations, as does the importance of minimising the manual
intervention required to run the burst detectors, such as through autonomous
parameter selection. The development of online burst detection techniques that
can detect bursting activity in real time is also necessary to facilitate
areas such as the study of real time learning in embodied cultured networks,
and applications involving bidirectional communication between biological
tissue and computer interfaces [94, 5]. This is another area in which ideas
adopted from burst detectors developed outside of neuroscience may benefit the
field.

In conclusion, years of study of bursting activity in cultured neuronal
networks has led to the development of many promising burst detection methods.
However, a ‘perfect’ method for analysing bursting activity remains elusive.
In the future, the development of improved burst detection methods will be
essential to keep up with advances in experimental techniques used to record
bursting activity, such as the use of higher density arrays and availability
of recordings from human stem cell-derived networks.

###  3.9 Acknowledgements

EC was supported by a Wellcome Trust PhD Studentship and a National Institute
for Health Research (NIHR) Cambridge Biomedical Research Centre Studentship.

###  3.10 Other resources

  * •

Open source R code for the burst detection methods outlined in this chapter
are available at https://github.com/ellesec/burstanalysis and archived at
https://doi.org/10.5281/zenodo.1284064.

## References

  * [1] M. Abeles, H. Bergman, I. Gat, I. Meilijson, E. Seidemann, N. Tishby, and E. Vaadia.  Cortical activity flips among quasi-stationary states.  Proc. Natl. Acad. Sci. U.S.A., 92:8616–8620, 1995.
  * [2] C. Allen and C. F. Stevens.  An evaluation of causes for unreliability of synaptic transmission.  Proc. Natl. Acad. Sci. U.S.A., 91:10380–10383, 1994.
  * [3] H. Amin, A. Maccione, F. Marinaro, S. Zordan, T. Nieus, and L. Berdondini.  Electrical responses and spontaneous activity of human iPS-derived neuronal networks characterized for 3-month culture with 4096-electrode arrays.  Front. Neurosci., 10:1–15, 2016.
  * [4] D. J. Bakkum, M. Radivojevic, U. Frey, F. Franke, A. Hierlemann, and H. Takahashi.  Parameters for burst detection.  Front. Comput. Neurosci., 7:193, 2013.
  * [5] D. J. Bakkum, A. C. Shkolnik, G. Ben-Ary, P. Gamblen, B. DeMarse, and S. M. Potter.  Removing some ’A’ from AI : Embodied Cultured Networks.  In F. Iida, R. Pfeifer, L. Steels, and Y. Kuniyoshi, editors, Embodied Artif. Intell., pages 130–146. Springer-Verlag, 2004.
  * [6] G. Barrionuevo, O. Benoit, and P. Tempier.  Evidence for two types of firing pattern during the sleep-waking cycle in the reticular thalamic nucleus of the cat.  Exp. Neurol., 72:486–501, 1981.
  * [7] Y. Ben-Ari.  Developing networks play a similar melody.  Trends Neurosci., 24:353–360, 2001.
  * [8] J. G. G. Borst.  The low synaptic release probability in vivo.  Trends Neurosci., 33:259–266, 2010.
  * [9] K. W. Boyack, K. Mane, and K. Börner.  Mapping Medline Papers, Genes and Proteins Related to Melanoma Research.  In Proc. Eigth Int. Conf. Inf. Vis., pages 965–971, 2004.
  * [10] T. Branco and K. Staras.  The probability of neurotransmitter release: variability and feedback control at single synapses.  Nat. Rev. Neurosci., 10:373–383, 2009.
  * [11] A. Burgos-Robles, I. Vidal-Gonzalez, E. Santini, and G. J. Quirk.  Consolidation of fear extinction requires NMDA receptor-dependent bursting in the ventromedial prefrontal cortex.  Neuron, 53:871–880, 2007.
  * [12] H. Câteau and A. D. Reyes.  Relation between single neuron and population spiking statistics and effects on network activity.  Phys. Rev. Lett., 96:058101, 2006.
  * [13] A. Cattaneo, L. Maffei, and C. Morrone.  Two firing patterns in the discharge of complex cells encoding different attributes of the visual stimulus.  Exp. Brain Res., 43:115–118, 1981.
  * [14] P. Charlesworth, E. Cotterill, A. Morton, S. G. Grant, and S. J. Eglen.  Quantitative differences in developmental profiles of spontaneous activity in cortical and hippocampal cultures.  Neural Dev., 10:1–10, 2015.
  * [15] P. Charlesworth, A. Morton, S. J. Eglen, N. H. Komiyama, and S. G. N. Grant.  Canalization of genetic and pharmacological perturbations in developing primary neuronal activity patterns.  Neuropharmacology, 100:47–55, 2016.
  * [16] L. Chen, Y. Deng, W. Luo, Z. Wang, and S. Zeng.  Detection of bursts in neuronal spike trains by the mean inter-spike interval method.  Prog. Nat. Sci., 19(2):229–235, 2009.
  * [17] Z. Chen and E. N. Brown.  Discrete- and Continuous-Time Probabilistic Models and Algorithms for Infferring Neuronal UP and DOWN States.  Neural Comput., 21(7):1797–1862, 2009.
  * [18] M. Chiappalone, M. Bove, A. Vato, M. Tedesco, and S. Martinoia.  Dissociated cortical networks show spontaneously correlated activity patterns during in vitro development.  Brain Res., 1093:41–53, 2006.
  * [19] M. Chiappalone, A. Novellino, I. Vajda, A. Vato, S. Martinoia, and J. van Pelt.  Burst detection algorithms for the analysis of spatio-temporal patterns in cortical networks of neurons.  Neurocomputing, 65-66:653–662, 2005.
  * [20] C. Chiu and M. Weliky.  Spontaneous activity in developing ferret visual cortex in vivo.  J. Neurosci., 21:8906–8914, 2001.
  * [21] J. H. Cocatre-Zilgien and F. Delcomyn.  Identification of bursts in spike trains.  J. Neurosci. Methods, 41(1):19–30, 1992.
  * [22] E. Cotterill, P. Charlesworth, C. W. Thomas, O. Paulsen, and S. J. Eglen.  A comparison of computational methods for detecting bursts in neuronal spike trains and their application to human stem cell-derived neuronal networks.  J. Neurophysiol., 116:306–321, 2016.
  * [23] L. N. Eisenman, C. M. Emnett, J. Mohan, C. F. Zorumski, and S. Mennerick.  Quantification of bursting and synchrony in cultured hippocampal neurons.  J. Neurophysiol., 114:1059–1071, 2015.
  * [24] J. Epsztein, M. Brecht, and A. K. Lee.  Intracellular determinants of hippocampal CA1 place and silent cell activity in a novel environment.  Neuron, 70:109–120, 2011.
  * [25] E. V. Evarts.  Temporal patterns of discharge of pyramidal tract neurons during sleep and waking in the monkey.  J. Neurophysiol., 27:152–171, 1964.
  * [26] R. C. Froemke, A. Tsay, I, M. Raad, J. D. Long, and Y. Dan.  Contribution of individual spikes in burst-induced long-term synaptic Modification.  J. Neurophysiol., 95:1620–1629, 2006.
  * [27] F. Gabbiani, W. Metzner, R. Wessel, and C. Koch.  From stimulus encoding to feature extraction in weakly electric fish.  Nature, 384:563–567, 1996.
  * [28] K. H. Gilchrist, G. F. Lewis, E. A. Gay, K. L. Sellgren, and S. Grego.  High-throughput cardiac safety evaluation and multi-parameter arrhythmia profiling of cardiomyocytes using microelectrode arrays.  Toxicol. Appl. Pharmacol., 288:249–257, 2015.
  * [29] A. Golbs, B. Nimmervoll, J.-J. Sun, I. E. Sava, and H. J. Luhmann.  Control of programmed cell death by distinct electrical activity patterns.  Cereb. Cortex, 21:1192–1202, 2011.
  * [30] B. Gourévitch and J. J. Eggermont.  A nonparametric approach for detection of bursts in spike trains.  J. Neurosci. Methods, 160(2):349–358, 2007.
  * [31] R. E. Harris, M. G. Coulombe, and M. B. Feller.  Dissociated retinal neurons form periodically active synaptic circuits.  J. Neurophysiol., 88:188–195, 2002.
  * [32] N. Heck, A. Golbs, T. Riedemann, J.-J. Sun, V. Lessmann, and H. J. Luhmann.  Activity-dependent regulation of neuronal apoptosis in neonatal mouse cerebral cortex.  Cereb. Cortex, 18:1335–1349, 2008.
  * [33] T. J. Heikkilä, L. Ylä-Outinen, J. M. A. Tanskanen, R. S. Lappalainen, H. Skottman, R. Suuronen, J. E. Mikkonen, J. A. K. Hyttinen, and S. Narkilahti.  Human embryonic stem cell-derived neuronal cells form spontaneously active neuronal networks in vitro.  Exp. Neurol., 218(1):109–116, 2009.
  * [34] M. H. Hennig, J. Grady, J. van Coppenhagen, and E. Sernagor.  Age-dependent homeostatic plasticity of GABAergic signaling in developing retinal networks.  J. Neurosci., 31(34):12159–12164, 2011.
  * [35] M. Ichikawa, K. Muramoto, K. Kobayashi, M. Kawahara, and Y. Kuroda.  Formation and maturation of synapses in primary cultures of rat cerebral cortical cells: an electron microscopic study.  Neurosci. Res., 16:95–103, 1993.
  * [36] S. Illes, W. Fleischer, M. Siebler, H.-P. Hartung, and M. Dihné.  Development and pharmacological modulation of embryonic stem cell-derived neuronal network activity.  Exp. Neurol., 207:171–176, 2007.
  * [37] M. E. Jackson, H. Homayoun, and B. Moghaddam.  NMDA receptor hypofunction produces concomitant firing rate potentiation and burst activity reduction in the prefrontal cortex.  Proc. Natl. Acad. Sci. U.S.A., 101:8467–8472, 2004.
  * [38] H. Kamioka, E. Maeda, Y. Jimbo, H. P. C. Robinson, and A. Kawana.  Spontaneous periodic synchronized bursting during formation of mature patterns of connections in cortical cultures.  Neurosci. Lett., 206:109–112, 1996.
  * [39] Y. Kaneoke and J. L. Vitek.  Burst and oscillation as disparate neuronal properties.  J. Neurosci. Methods, 68(2):211–223, 1996.
  * [40] F. E. Kapucu, J. M. A. Tanskanen, J. E. Mikkonen, L. Ylä-Outinen, S. Narkilahti, and J. A. K. Hyttinen.  Burst analysis tool for developing neuronal networks exhibiting highly varying action potential dynamics.  Front. Comput. Neurosci., 6:38, 2012.
  * [41] J. Kleinberg.  Bursty and Hierarchical Structure in Streams.  In Proc. 8th ACM SIGKDD Int. Conf. Knowl. Discov. Data Min., pages 91–101, 2002.
  * [42] D. Ko, C. J. Wilson, C. J. Lobb, and C. A. Paladini.  Detection of bursts and pauses in spike trains.  J. Neurosci. Methods, 211(1):145–158, 2012.
  * [43] R. Krahe and F. Gabbiani.  Burst firing in sensory systems.  Nat. Rev. Neurosci., 5:13–23, 2004.
  * [44] R. Krahe, G. Kreiman, F. Gabbiani, C. Koch, and W. Metzner.  Stimulus encoding and feature extraction by multiple sensory neurons.  J. Neurosci., 22:2374–2382, 2002.
  * [45] R. Kumar, H. Road, S. Jose, H. Road, S. Jose, R. Drive, A. Tomkins, and S. Jose.  On the Bursty Evolution of Blogspace.  In Int. World Wide Web Conf., pages 568–576, 2003.
  * [46] C. R. Legéndy and M. Salcman.  Bursts and recurrences of bursts in the spike trains of spontaneously active striate cortex neurons.  J. Neurophysiol., 53(4):926–939, 1985.
  * [47] X. Leinekugel, R. Khazipov, R. Cannon, H. Hirase, Y. Ben-Ari, and G. Buzsáki.  Correlated bursts of activity in the neonatal hippocampus in vivo.  Science, 296:2049–2052, 2002.
  * [48] J. E. Lisman.  Bursts as a unit of neural information: Making unreliable synapses reliable.  Trends Neurosci., 20(1):38–43, 1997.
  * [49] C. J. Lobb.  Abnormal bursting as a pathophysiological mechanism in Parkinson’s disease.  Basal Ganglia, 3:187–195, 2014.
  * [50] D. Lonardoni, S. Di Marco, H. Amin, A. Maccione, L. Berdondini, and T. Nieus.  High-density MEA recordings unveil the dynamics of bursting events in cell cultures.  Conf. Proc. IEEE Eng. Med. Biol. Soc., 2015:3763–3766, 2015.
  * [51] A. Maccione, M. H. Hennig, M. Gandolfo, O. Muthmann, J. van Coppenhagen, S. J. Eglen, L. Berdondini, and E. Sernagor.  Following the ontogeny of retinal waves: pan-retinal recordings of population dynamics in the neonatal mouse.  J. Physiol., 592(7):1545–1563, 2014.
  * [52] E. Maeda, H. P. Robinson, and A. Kawana.  The mechanisms of generation and propagation of synchronized bursting in developing networks of cortical neurons.  J. Neurosci., 15(10):6834–6845, 1995.
  * [53] M. B. Martens, M. Chiappalone, D. Schubert, and P. H. E. Tiesinga.  Separating burst from background spikes in multichannel neuronal recordings using return map analysis.  Int. J. Neural Syst., 24(04):1450012, 2014.
  * [54] J. Martinson, H. H. Webster, A. A. Myasnikov, and R. W. Dykes.  Recognition of temporally structured activity in spontaneously discharging neurons in the somatosensory cortex in waking cats.  Brain Res., 750:129–140, 1997.
  * [55] A. Mazzoni, F. D. Broccard, E. Garcia-Perez, P. Bonifazi, M. E. Ruaro, and V. Torre.  On the dynamics of the spontaneous activity in neuronal networks.  PLoS One, 2:e439, 2007.
  * [56] R. W. McCarley, O. Benoit, and G. Barrionuevo.  Lateral geniculate nucleus unitary discharge in sleep and waking: State- and rate-specific aspects.  J. Neurophysiol., 50:798–818, 1983.
  * [57] M. Meister, R. O. L. Wong, D. A. Baylor, and C. J. Shatz.  Synchronous bursts of action potentials in ganglion cells of the developing mammalian retina.  Science, 252:939–943, 1991.
  * [58] B. R. Miller, A. G. Walker, S. J. Barton, and G. V. Rebec.  Dysregulated Neuronal Activity Patterns Implicate Corticostriatal Circuit Dysfunction in Multiple Rodent Models of Huntington’s Disease.  Front. Syst. Neurosci., 5:26, 2011.
  * [59] Nex Technologies.  NeuroExplorer Manual.  Nex Technologies, 2014.
  * [60] Z. G. Ni, R. Bouali-Benazzouz, D. M. Gao, A. L. Benabid, and A. Benazzouz.  Time-course of changes in firing rates and firing patterns of subthalamic nucleus neuronal activity after 6-OHDA-induced dopamine depletion in rats.  Brain Res., 899:142–147, 2001.
  * [61] J. Nicolas, P. J. M. Hendriksen, R. G. D. M. van Kleef, A. de Groot, T. F. H. Bovee, I. M. C. M. Rietjens, and R. H. S. Westerink.  Detection of marine neurotoxins in food safety testing using a multielectrode array.  Mol. Nutr. Food Res., 58:2369–2378, 2014.
  * [62] A. Odawara, H. Katoh, N. Matsuda, and I. Suzuki.  Physiological maturation and drug responses of human induced pluripotent stem cell-derived cortical neuronal networks in long-term culture.  Sci. Rep., 6:1–14, 2016.
  * [63] J. O’Keefe and M. L. Recce.  Phase relationship between hippocampal place units and the EEG theta rhythm.  Hippocampus, 3:317–330, 1993.
  * [64] T. Otto, H. Eichenbaum, C. G. Wible, and S. I. Wiener.  Learning‐-related patterns of CA1 spike trains parallel stimulation parameters optimal for inducing hippocampal long‐term potentiation.  Hippocampus, 1:181–192, 1991.
  * [65] V. Pasquale, S. Martinoia, and M. Chiappalone.  A self-adapting approach for the detection of bursts and network bursts in neuronal cultures.  J. Comput. Neurosci., 29(1-2):213–229, 2010.
  * [66] V. Pasquale, P. Massobrio, L. L. Bologna, M. Chiappalone, and S. Martinoia.  Self-organization and neuronal avalanches in networks of dissociated cortical neurons.  Neuroscience, 153:1354–1369, 2008.
  * [67] O. Paulsen and T. J. Sejnowski.  Natural patterns of activity and long-term synaptic plasticity.  Curr. Opin. Neurobiol., 10:172–179, 2000.
  * [68] F. G. Pike, R. M. Meredith, A. W. A. Olding, and O. Paulsen.  Postsynaptic bursting is essential for ’Hebbian’ induction of associative long-term potentiation at excitatory synapses in rat hippocampus.  J. Physiol., 518:571–576, 2004.
  * [69] A. Pimashkin, I. Kastalskiy, A. Simonov, E. Koryagina, I. Mukhina, and V. Kazantsev.  Spiking signatures of spontaneous activity bursts in hippocampal cultures.  Front. Comput. Neurosci., 5:1–12, 2011.
  * [70] S. Pluta, A. Naka, J. Veit, G. Telian, L. Yao, R. Hakim, D. Taylor, and H. Adesnik.  A direct translaminar inhibitory circuit tunes cortical output.  Nat. Neurosci., 18:1631–1640, 2015.
  * [71] G. Radons, J. D. Becker, B. Dülfer, and J. Krüger.  Analysis, classification, and coding of multielectrode spike trains with hidden markov models.  Biol. Cybern., 71:359–373, 1994.
  * [72] N. Raichman and E. Ben-Jacob.  Identifying repeating motifs in the activation of synchronized bursts in cultured neuronal networks.  J. Neurosci. Methods, 170:96–110, 2008.
  * [73] B. K. Rhoades and G. W. Gross.  Potassium and calcium channel dependence of bursting in cultured neuronal networks.  Brain Res., 643:310–318, 1994.
  * [74] W. Schultz.  Predictive reward signal of dopamine neurons.  J. Neurophysiol., 80:1–27, 1998.
  * [75] W. Schultz, P. Dayan, and P. R. Montague.  A neural substrate of prediction and reward.  Science, 275:1593–1599, 1997.
  * [76] J. V. Selinger, N. V. Kulagina, T. J. O’Shaughnessy, W. Ma, and J. J. Pancrazio.  Methods for characterizing interspike intervals and identifying bursts in neuronal activity.  J. Neurosci. Methods, 162(1-2):64–71, 2007.
  * [77] V. Senn, S. B. E. Wolff, C. Herry, F. Grenier, I. Ehrlich, J. Gründemann, J. P. Fadok, C. Müller, J. J. Letzkus, and A. Lüthi.  Long-range connectivity defines behavioral specificity of amygdala neurons.  Neuron, 81:428–437, 2014.
  * [78] S. M. Sherman.  Tonic and burst firing: Dual modes of thalamocortical relay.  Trends Neurosci., 24:122–126, 2001.
  * [79] A. Singh, K. Mewes, R. E. Gross, M. R. DeLong, J. A. Obeso, and S. M. Papa.  Human striatal recordings reveal abnormal discharge of projection neurons in parkinson’s disease.  Proc. Natl. Acad. Sci. U.S.A., 113:9629–9634, 2016.
  * [80] M. Steriade, I. Timofeev, and F. Grenier.  Natural waking and sleep states: a view from inside neocortical neurons.  J. Neurophysiol., 85:1969–1985, 2001.
  * [81] D. Tam.  An alternate burst analysis for detecting intra-burst firings based on inter-burst periods.  Neurocomputing, 46:1155–1159, 2002.
  * [82] M. J. Thomas, A. M. Watabe, T. D. Moody, M. Makhinson, and T. J. O’Dell.  Postsynaptic complex spike bursting enables the induction of LTP by theta frequency synaptic stimulation.  J. Neurosci., 18:7118–7126, 1998.
  * [83] A. M. Thomson.  Activity-dependent properties of synaptic transmission at two classes of connections made by rat neocortical pyramidal axons in vitro.  J. Physiol., 502:131–147, 1997.
  * [84] P. N. Tobler, A. Dickinson, and W. Schultz.  Coding of predicted reward omission by dopamine neurons in a conditioned inhibition paradigm.  J. Neurosci., 23:10402–10410, 2003.
  * [85] S. Tokdar, P. Xi, R. C. Kelly, and R. E. Kass.  Detection of bursts in extracellular spike trains using hidden semi-Markov point process models.  J. Comput. Neurosci., 29(1-2):203–212, 2010.
  * [86] L. Turnbull, E. Dian, and G. Gross.  The string method of burst identification in neuronal spike trains.  J. Neurosci. Methods, 145(1-2):23–35, 2005.
  * [87] P. Valdivia, M. Martin, W. R. LeFew, J. Ross, K. A. Houck, and T. J. Shafer.  Multi-well microelectrode array recordings detect neuroactivity of ToxCast compounds.  Neurotoxicology, 44:204–217, 2014.
  * [88] I. A. Välkki, K. Lenk, J. E. Mikkonen, and F. E. Kapucu.  Network-wide adaptive burst detection depicts neuronal activity with improved accuracy.  Front. Comput. Neurosci., 11:40, 2017.
  * [89] A. N. Van Den Pol, K. Obrietan, and A. Belousov.  Glutamate hyperexcitability and seizure-like activity throughout the brain and spinal cord upon relief from chronic glutamate receptor blockage in culture.  Neuroscience, 74:653–674, 1996.
  * [90] F. Van Huizen, H. J. Romijn, and A. M. M. C. Habets.  Synaptogenesis in Rat Cerebral Cortex Cultures is Affected During Chronic Blockade of Spontaneous Bioelectric Activity by Tetrodotoxin.  Dev. Brain Res., 19:67–80, 1985.
  * [91] J. Van Pelt, M. A. Corner, P. S. Wolters, W. L. C. Rutten, and G. J. A. Ramakers.  Longterm stability and developmental changes in spontaneous network burst firing patterns in dissociated rat cerebral cortex cell cultures on multielectrode arrays.  Neurosci. Lett., 361:86–89, 2004.
  * [92] J. Van Pelt, I. Vajda, P. S. Wolters, M. a. Corner, and G. J. a. Ramakers.  Dynamics and plasticity in developing neuronal networks in vitro.  Prog. Brain Res., 147:173–188, 2005.
  * [93] J. Van Pelt, P. S. Wolters, M. A. Corner, W. L. C. Rutten, and G. J. A. Ramakers.  Long-term characterization of firing dynamics of spontaneous bursts in cultured neural networks.  IEEE Trans. Biomed. Eng., 51:2051–2062, 2004.
  * [94] D. Wagenaar, T. B. Demarse, and S. M. Potter.  MeaBench : A toolset for multi-electrode data acquisition and on-line analysis.  In Proc. 2nd Int. IEEE EMBS Conf. Neural Eng., pages 518–521, 2005\.
  * [95] D. A. Wagenaar, J. Pine, and S. M. Potter.  An extremely rich repertoire of bursting patterns during the development of cortical cultures.  BMC Neurosci., 7:11, 2006.
  * [96] A. G. Walker, B. R. Miller, J. N. Fritsch, S. J. Barton, and G. V. Rebec.  Altered information processing in the prefrontal cortex of huntington’s disease mouse models.  J. Neurosci., 28:8973–8982, 2008.
  * [97] M. Weliky and L. C. Katz.  Correlational structure of spontaneous neuronal activity in the developing lateral geniculate nucleus in vivo.  Science, 285:599–604, 1999.
  * [98] T. G. Weyand, M. Boudreaux, and W. Guido.  Burst and tonic response modes in thalamic neurons during sleep and wakefulness.  J. Neurophysiol., 85(3):1107–1118, 2001.
  * [99] W. Xu, W. Morishita, P. S. Buckmaster, Z. P. Pang, R. C. Malenka, and T. C. Südhof.  Distinct neuronal coding schemes in memory revealed by selectrive erasure of fast synchronous synaptic transmission.  Neuron, 73:990–1001, 2012.
  * [100] L. Ylä-Outinen, J. Heikkilä, H. Skottman, R. Suuronen, R. Aänismaa, and S. Narkilahti.  Human cell-based micro electrode array platform for studying neurotoxicity.  Front. Neuroeng., 3:1–9, 2010.
  * [101] X. Zhang and D. Shasha.  Better burst detection.  In Proc. 22nd Int. Conf. Data Eng., page 146, 2006.
  * [102] Y. Zhu and D. Shasha.  Efficient elastic burst detection in data streams.  In Proc. Ninth ACM SIGKDD Int. Conf. Knowl. Discov. data Min., pages 336–345, 2003.