# Efficient Epileptic Seizure Detection Using CNN-Aided Factor Graphs

Bahareh Salafian1, Eyal Fishel Ben2, Nir Shlezinger2, Sandrine de
Ribaupierre1,3, and Nariman Farsad4 1B. Salafian is with the School of
Biomedical Engineering, University of Western Ontario, London, ON N6A 5B9
<EMAIL_ADDRESS>Fishel and N. Shlezinger are with the School of Electrical
and Computer Engineering, Ben-Gurion University of the Negev, Be’er Sheva,
Israel, 84105<EMAIL_ADDRESS><EMAIL_ADDRESS>de Ribaupierre is
with Department of Clinical Neurological Sciences and the School of Biomedical
Engineering, University of Western Ontario, ON N6A 5B9<EMAIL_ADDRESS>Farsad is with the Department of Computer Science at Ryerson
University,Toronto, ON M5B 2K3<EMAIL_ADDRESS>
###### Abstract

We propose a computationally efficient algorithm for seizure detection.
Instead of using a purely data-driven approach, we develop a hybrid model-
based/data-driven method, combining convolutional neural networks with factor
graph inference. On the CHB-MIT dataset, we demonstrate that the proposed
method can generalize well in a 6 fold leave-4-patient-out evaluation.
Moreover, it is shown that our algorithm can achieve as much as 5% absolute
improvement in performance compared to previous data-driven methods. This is
achieved while the computational complexity of the proposed technique is a
fraction of the complexity of prior work, making it suitable for real-time
seizure detection.

##  I Introduction

Epilepsy is one of the most common neurological disorders affecting about 50
million people worldwide. This disorder is associated with recurrent episodes
of abnormal neural activity in the central nervous system known as epileptic
seizures [1]. Based on the area where seizure starts and the intensity of
brain’s abnormal signals, patients with epilepsy may suffer from different
symptoms including auras, muscle contraction, and loss of consciousness [2].
Epilepsy affects the patients’ private and professional life; for instance,
some activities such as swimming, bathing, and climbing a ladder, become
dangerous as a seizure during that activity might result in unpredictable
injuries, and even death. Therefore, early detection of epilepsy can notably
improve the patient’s quality of life.

The most common tool used to diagnose seizures is electroencephalogram (EEG)
[3]. While other various techniques, such as magnetic resonance imaging (MRI)
[4], magnetoencephalography (MEG) [5], and positron emission tomography (PET)
[6] are sometimes used in conjunction to EEG, EEG is widely preferred as it is
economical, portable, non-invasive and shows clear rhythms in the frequency
domain [7]. However, the review of EEG signals is a time-consuming process, as
a neurologist needs to monitor the recording. Expertise is needed to diagnose
epilepsy as each case is quite variable, with different channels involved,
while the spectral content of the rhythmic activity varies across individuals
and signals, and is mostly contaminated by physiologic and non-physiologic
interference [8]. Hence, automatic seizure detection is a valuable clinical
tool to address this issue and reduce the dependency on human experts.

Many machine learning studies have been developed for automatic seizure
detection problems. One of the most common approaches applies a support vector
machine (SVM), which is mainly followed by additional preprocessing steps such
as discrete wavelet transform (DWT) and fast Fourier transform (FFT) to
extract more features of EEG signals [9, 10, 11, 12]. In the past decade, deep
learning (DL) techniques have become very popular in various applications,
including the analysis of time series EEG signals. Therefore, different DL
models have been investigated and tested in the area of seizure detection. For
instance, Khalilpour et al. [13] applied a 1D convolutional neural network
(CNN) to EEG signals of five patients to predict preictal and interictal
states of the brain. In [14], the spectrogram of EEG measurements was used as
input to a 1D CNN. The authors in [15] utilized DWT to represent the EEG
segments, which is used as input to a 1D CNN. Moreover, they combine the DWT
of the current, previous, and next block for predicting the label of the
current block, which helps exploit the temporal correlations. Another common
DL architecture is 2D CNN. Boonyakitanont et al. [16] applied 2D CNN to 24
epileptic recordings from 23 patients, where the signals are segmented into
4-second blocks. They showed state-of-the-art performance in terms of
detection accuracy.

Most of these prior techniques divide the EEG recording into blocks and treat
these blocks independently. This does not take advantage of the temporal
correlations that exists between consecutive blocks. While there are methods
based on CNN\- recurrent neural network (RNN) architectures [17, 18] that can
mitigate this issue, it is well known that RNNs have high computational
complexity for training. Another method used for exploiting temporal
correlations is hidden Markov model (HMM) [19]. HMMs belong to family of
factorizable joint distributions which admit low-complexity inference via
factor graph methods [20]. Despite proliferation of seizure detection
algorithms, having a computationally efficient algorithm that can generalize
to different patients and perform seizure detection reliably in a real-time
manner is lacking.

In this work, we propose a computationally efficient epileptic seizure
detection algorithm based on a hybrid model-based/data-driven approach using
CNN-aided factor graphs. First, we carefully design a 1D CNN for estimating
the probability that a 4-second block of EEG is a seizure block. Our goal is
to design a network that is applied to the EEG signals directly, without
feature engineering using transforms such as DWT or FFT. When using such
feature engineering, one must carefully select the parameters of the
transforms, and the performance of the trained model can vary considerably
based on these parameters. Moreover, while these transforms typically provide
information about the frequency component of the EEG signals, they are not
necessarily impactful for capturing the dependence of the EEG signals across
different channels. Such a dependence can be indicative of epileptic seizure.
Our proposed 1D CNN is designed to be able to capture the long term dependence
between EEG channels and operates on the signals directly with minimal
processing. To exploit the temporal correlation between consecutive blocks for
further improvement, we use factor graph inference, specifically using HMM
models to capture temporal correlations among the signals. Our proposed hybrid
method is highly efficient compared to previous works, where we reduce the
inference complexity by a factor of 2. Despite this decrease in computation,
our method achieves up to 5% absolute improvement in performance measures such
as precision, recall, and F1-score in a 6-fold leave-4-patients-out
evaluation.

The rest of this paper is organized as follows. In Section II we describe the
problem statement, the dataset that is used for model development and
evaluation, and the baseline models. Then, in Section III we describe our
proposed method, which is evaluated and compared to prior methods in Section
IV. Finally, Section V provides concluding remarks.

##  II Background, Dataset, and Baseline Models

In this section, we first discuss the seizure detection problem. We then
describe the dataset that is used for our hybrid model-based/data-driven
algorithm development and evaluation. We conclude this section by presenting
the baseline methods that are used for comparison in this work.

###  II-A Seizure Detection Using EEG Signals

EEG is the electrical recording of the brain activities which is the most
popular diagnostic and analytical tool for epileptic seizures. In seizure
detection problem, there are some basic terms as follows:

  * •

Ictal state which is the time when the seizure occurs (from start to end).

  * •

Preictal state is a period of time just before a seizure occurs.

  * •

Interictal state refers to the period between seizures.

  * •

Postictal state is a period of time just after the seizure ends.

During reading an EEG, experts don’t only look at seizures which are not
always occurring when the EEG is done, specifically with short EEG recordings,
but they also have to interpret interictal EEG signal and therefore, they are
trying to detect interictal epileptiform discharges (IED)[21]. IEDs are
generated by the synchronous discharges of a group of neurons in a region
referred to as the epileptic focus; however, the detection of these spikes is
difficult to accomplish due to their similarity to waves that are part of
normal EEG or artifacts and the wide variability in spike morphology and
background between patients. For instance, abnormalities like breach rhythm
(normal rhythm seen with skull defects) can have focal, sharply contoured
morphology and they might be inferred as epileptic seizures. As such, having
an automatic system to detect and predict seizures will resolve these issues.

###  II-B Data Description

In this section we detail the data used in our study of hybrid model-
based/data-driven epileptic seizure detection. We first describe the raw EEG
data, after which we describe the pre-processing carried out prior to its
usage for training and inference.

####  II-B1 EEG Data

The dataset used in this study is the public CHB-MIT Scalp EEG Database
collected at the Children’s Hospital Boston and consists of EEG recordings
from pediatric subjects with intractable seizures [22]111This database is
available online at PhysioNet
(https://physionet.org/physiobank/database/chbmit/). Recordings were collected
from 23 subjects: 5 males aged 3-22 years, 17 females aged 1.5- 19 years, and
one anonymous subject. Each case contains 9 to 42 continuous EDF files from a
single subject. The duration of the recordings in each file varies between one
to four hours. All signals were sampled at frequency of 256 Hz with 16-bit
resolution. Since case 21 was obtained 1.5 years after case 1 from the same
female subject, we consider case 21 as a separate patient; therefore, our
experiment includes 24 subjects. Please note that since we are evaluating the
algorithms using a 6-fold leave-4-patients-out method, this might negligibly
effect only one of the folds.

####  II-B2 Data Pre-Processing

Since many of patient recordings do not contain any seizures, in order to have
a more balanced samples from seizures, for each patient, we only selected EDF
files that have at least one seizure. Moreover, since the length of the
seizures are very short (from 7 seconds to 753 seconds) compared to the
overall recording (from 959 seconds to 14427 seconds), we shorten the
recording to 3 times the seizure duration before and 3 times the seizure
duration after the seizure. Therefore, for every second of seizure data, there
are 6 seconds of non-seizure data.

From the EEG channels, we use the 18 bipolar montage: FP1-F7, F7-T7, T7-P7,
P7-O1, FP1-F3, F3-T3, T3-P3, P3-O1, FP2-F4, F4-C4, C4-P4, P4-O2, FP2-F8,
F8-T8, T8-P8, P8-O2, FZ-CZ, CZ-PZ. A notch filter is used to remove 60 Hz line
noise from each EEG signal. Then 4-second blocks, with 1024 sample points per
block, are used as moving window with step size of 1 second. The value of 4
seconds was chosen to provide a good trade off between the number of samples
in a block and the stationarity of the observed signals over a block[23]. We
have observed that when the width of window is increased, the seizure
detection procedure is not accurate enough. We now describe some of the prior
seizure detection algorithms developed using this dataset, which will be used
as baselines in this paper.

###  II-C Baseline Methods

We consider two recent works that presented state-of-the-art performance on
the CHB-MIT dataset as baselines in these papers [16, 24]. The method in [16]
is designed to take the EEG signals as input without applying any type of
transforms. Therefore, we employ the same structure used in [16] where the
input shape for the model is (18,1024,1), which implies considering each block
as an image with the size of (18,1024) and channel dimension of 1. The method
in [24] takes a graphical image of the EEG recording as input rather than the
EEG measurements. This type of image-based feature was shown to achieve better
detection performance compared to other features-based methods in seizure
detection such as spectrogram or periodogram in a recent study [25].

##  III Methodology

Our hybrid model-based/data-driven algorithm combines a carefully designed
CNN, which estimates the presence of seizure in a 4-second block, with factor
graph inference to exploit the temporal correlation between the blocks. We
begin this section by describing the the CNN architecture, followed by the
factor graph based inference step.

###  III-A 1D CNN Architecture

2D CNNs, as utilized in the baseline methods detailed in Subsection II-C,
exploit the notion of locality, exhibited by natural images, which implies
that the level of correlation between different elements typically grows the
closer they are in the image. However, in EEG segment matrices, each element
represents a single EEG measurement, which is likely to be correlated with all
the remaining measurements taken at that time instance, regardless of their
row index in the matrix representation. This structure makes 1D CNNs, which
combine all measurements taken from different channels at a given time
instance, more suitable compared to 2D CNNs used in [16]. That allows the
network to better learn the correlations that may exist between different
channels since during seizures, the channel measurements can become highly
correlated.

In order to have a comparable configuration with the baseline models, we use
the same number of layers where the input shape for this model is (1024,18).
The inputs to our 1D CNN are EEG signals with minimal prepossessing, which
only removes the 60 Hz component of the signals. The baseline models detailed
in Subsection II-C use a kernel size of 3 and 2. Given the number of layers,
this results in a receptive field of approximately 30 milliseconds. Since this
is not enough to capture low frequency components of the signal as well as the
long term correlations between the EEG channels, we design out kernel size to
be much larger, which results in a receptive field that covers approximately 1
second of the data. Fig. 1 shows the complete architecture of the proposed 1D
CNN.

As will be shown in Section IV, these simple changes, namely choosing the
kernel size carefully and using a 1D CNN, improve the results significantly
compared to the baseline CNN model in [16]. This is because the network can
capture a wider range of frequency components in the signal and also better
capture the correlations between the EEG channels. An additional advantage of
using 1D CNNs stems from their reduced complexity during inference compared to
their 2D counterparts. Although we have increased the receptive field of the
network by increasing the kernel size for our proposed 1D CNN architecture
compared to [16], the number of floating point operations (FLOPs) during
inference (i.e., the number of floating point multiplication and summation
operations) for the proposed method is almost down by a factor of 2 compared
to [16] and by a factor of 20 compared to [24] as summarized in Table I. For
training all networks, we use the ADAM optimizer with the learning rate of
0.001, batch size of 128 and 10 epochs of training.

TABLE I: Computational complexity in FLOPs for all models | Mega FLOPs
---|---
2D CNN [16] | 14.5
2D CNN [24] | 200
1D CNN | 9.81
1D CNN+FG | 9.81
1D CNN+GRU | 29.4
Figure 1: Proposed 1D CNN architecture.

###  III-B Factor Graph Based Inference

The 1D CNN model outputs an estimate of the probability that a given 4-seconds
block corresponds to a seizure. This probability is based solely on its input
EEG segment block, and does not account for the fact that the presence of a
seizure in a given block is likely to also reflect on its preceding and
subsequent blocks. To incorporate this temporal correlation, we combine the
probability estimates over multiple blocks by assuming that the underlying
temporal correlation can be represented using a factor graph, and utilize the
sum-product method for inference [26]. In the following we first describe how
the underlying dynamics of the seizure detection setup can be represented as a
factor graph, after which we discuss the sum-product algorithm and elaborate
on its combinations with the proposed 1-D CNN model.

####  III-B1 Factor Graph Representation of Underlying Dynamics

Factor graphs provide a visual representation of a multivariate function,
typically a joint distribution measure, which can be factorized into a
partition of local functions [20]. These partitions capture the inherent
statistical relationship among variables affecting each partition. To
represent a multivariate function as a factor graph, every partition and every
variable must be associated with a unique node. Edges connect function nodes
to variables nodes if and only if the function is explicitly dependent on the
corresponding variable. We adopt the Forney style representation of factor
graphs, where variable nodes are replaced by edges [27]. This graphical
representation enables desired quantities to be computed at reduced complexity
via message passing over the factor graph [20].

To implement factor graphs inference, we first fix the structure of the graph,
i.e., the interconnection between its nodes, which encapsulates our knowledge
on the underlying statistical relationships for the seizure detection problem.
A key feature preceding epileptic seizures is the de-synchronization of its
rhythmic activity [28]. To capture this behavior in our model study, we adopt
a first-order HMM. The Markovian architecture focuses on the temporal
characteristics associated with the epileptic episodes, highlights the effects
of temporal correlations in seizure detection, and facilitates efficient
classification at reduced complexity. The resulting factor graph under this
model is illustrated in Fig. 2.

Figure 2: Factor Graph of a HMM with function nodes computed using a 1-D CNN.

To formulate this mathematically, let
$\boldsymbol{y}=\\{y_{1},y_{2},\cdots,y_{N}\\}$, and
$\boldsymbol{s}=\\{s_{1},s_{2},\cdots,s_{N}\\}$ describe the observed EEG
measurements and latent seizure states, respectively, over $N$ consecutive
4-seconds blocks. The latent states takes binary values, i.e.,
$s_{i}\in\\{0,1\\}$ corresponding to the presence or absence of a seizure. We
assume that these states satisfy the Markovian property, and use
$P(s_{i}|s_{i-1})$ to denote the transition probability. The transition
probability is taken as a control parameter, whose entries are fixed
handcrafted values. In particular, we set $P(s_{i}=1|s_{i-1}=0)$ to correspond
to a $10.46\%$ of switching from non-seizure to seizure, and $17.9\%$ for
transitioning in the opposite direction. Ideally, one would like the
transition probability matrix to reflect the true transitions probabilities
between seizure and non-seizure state. However, the negligible seizure
occurrences throughout the recording of EEG episodes, produces a highly
unbalanced transition probabilities. Empirically, such transition
probabilities did not facilitate accurate inference, hence hyper-parameter
optimization of these probabilities was adopt, so that superior performance
can be achieved. To close, we assume each EEG measurement depends only on its
corresponding seizure state. Under this postulated HMM, the joint probability
density function of the measurements and the states can be written as:

$P(\boldsymbol{s},\boldsymbol{y})=\prod_{i=1}^{N}P(s_{i}|s_{i-1})P(y_{i}|s_{i}).$
(1)

####  III-B2 The Sum-Product Algorithm

Having established the mathematical foundation and suitable factor graph
representation, the objective is to distinguish between seizure and non-
seizures states. Classification of these states is achieved through accurate
inference of the marginal distribution $P(s_{i},\boldsymbol{y})$, which is the
metric used to compute the maximum a-posteriori probability detector. In
principle, evaluating $P(s_{i},\boldsymbol{y})$ from (1) involves
marginalization; a task whose computational burden scales exponentially with
$N$. However, by employing factor graph inference via the sum-product
algorithm, the same computation scales only linearly with $N$, making this
operation computationally feasible for the problem at hand.

The sum-product method relates the desired marginal probability to a product
of "messages", where for each $k\in\\{1,\ldots,N\\}$ we write

$P(s_{k},\boldsymbol{y})=\mu_{f_{j}\rightarrow
s_{k}}(s_{k})\cdot\mu_{f_{j+1}\rightarrow s_{k}}(s_{k}).$ (2)

The terms $\mu_{f_{j}\rightarrow s_{k}}(s_{k}),\mu_{f_{j+1}\rightarrow
s_{k}}(s_{k})$ in (2) are interpreted as the "forward message" and "backward
message" respectively, and are computed by [20]

$\mu_{f_{j}\rightarrow
s_{k}}(s_{k})=\sum_{\\{s_{1},\cdots,s_{k-1}\\}}\prod_{i=1}^{n}f_{i}(y_{i},s_{i},s_{i-1}),$
(3)

and

$\mu_{f_{j+1}\rightarrow
s_{k}}(s_{k})=\sum_{\\{s_{k+1},\cdots,s_{N}\\}}\prod_{i=n+1}^{N}f_{i}(y_{i},s_{i},s_{i-1}),$
(4)

where,

$f_{i}(y_{i},s_{i},s_{i-1})=P(s_{i}|s_{i-1})P(y_{i}|s_{i}).$ (5)

For a Markov chain factor graph structure as in Fig. 2, the computational
complexity is comprised of evaluating the forward (3) and the backward (4)
messages. The recursive nature of the computations implies that the number of
FLOPs grows linearly with the number of $4$-second blocks in a given EEG
episode, denoted by $N$. At the same time, the number of FLOPs also increases
with the cardinality of the seizure state, denoted by $|\mathcal{S}|$, as well
as the order of the Markov chain. Here, all seizure states are binary, hence
$|S|=2$, and the factor graph is that of a first-order Markov chain. As a
result, each message computation requires 4 multiplications and 2 addition
operations. The result of these computations over a complete EEG episode
requires merely $12N$ FLOPs, which is negligible compared to the complexity of
applying the neural network models, as reported in Table I.

####  III-B3 CNN-Aided Factor Graphs

Implementing the sum-product algorithm requires knowledge of the underlying
probability distribution $P(y_{i}|s_{i})$. In the case of seizure detection
from EEG measurements, acquiring such a distribution from first principles is
unattainable. Accurately characterizing epileptic seizures to high fidelity
requires a combination of dynamical features and spatial-temporal features,
for example the entropy, amplitude, synchronization, spectral power amongst
other features associated with epileptic episodes, making the underlying
distribution highly complex and intractable. To contend with this difficultly,
we follow the work [29], combining DL models and inference algorithm. In
particular, we utilize the output of the 1D CNN as an estimate of the
conditional distribution $P(y_{i}|s_{i})$ required in order to compute the
function nodes222In principle, a CNN classifier is trained with the cross-
entropy loss to detect $s_{i}$ from $y_{i}$, which is an estimate of the
conditional distribution $P(s_{i}|y_{i})$. Here we use the CNN output as an
estimate of the conditional $P(y_{i}|s_{i})$, instead of converting it via
Bayes rule as done in [30], to avoid instabilities due to dividing by the
estimated marginal of the seizure state $s_{i}$.. via (5). Combining these
approaches yields a hybrid model-based/data-driven detector, in which factor
graph inference constitutes a robust final stage incorporating temporal
correlation with the CNN outputs as illustrated in Fig. 2.

To complete the picture, the detection mechanism requires a comparison
measure, i.e., a threshold, which represents a decision boundary for one to
distinguish between different seizure states based on the estimated marginal
distribution produced by the sum-product method. Treating the threshold as
fixed, if the resultant probability of a seizure state yielded by the CNN-
aided factor graph exceeds that of the applied threshold, then the system is
said to be in a seizure state. Otherwise the system is said to be in a non-
seizure state. As shown in the sequel, this approach allows to achieve
accurate detection with controllable tradeoffs between detection and false
alarm.

##  IV Results and Discussion

We now discuss our evaluation results in this section. We begin by defining
the performance metrics we employed for evaluation followed by the results and
the discussions.

###  IV-A Evaluation Method and Performance Metrics

We apply 6-fold leave-4-patient-out cross validation to evaluate the
performance and the generalizability of our hybrid model-based/data-driven
approach. Specifically, in each fold, 4 patients are kept for testing and 20
patients for training. The dataset used for evaluation, similar to train data,
is segmented into 4-second blocks with moving window of 1 second where each
block is labeled. To evaluate the performance of the models trained in each
fold, we use the following metrics.

  * •

F1-score: is a harmonic mean of recall and precision where former indicates
the proportion of real positive cases that are correctly predicted positive
and latter denotes the proportion of predicted positive cases that are
correctly real positives [31].

  * •

AUC-PR: is area under the precision-recall curve where high area under the
curve means low false positive rate and low false negative rate.

  * •

AUC-ROC is area under receiver operating characteristics (ROC) curve which is
performance measurement for classification problems at various threshold
settings.

We choose these performance measures since the data is imbalanced (i.e., for
every block with seizure we have 6 blocks with no seizures, and AUC measures
are great indicator of performance over imbalanced data.

TABLE II: Summary of Results | AUC-ROC | AUC-PR | F1 score
---|---|---|---
2D CNN [16] | $87.27\pm 0.05$ | $72.05\pm 11.46$ | $85.58\pm 1.63$
2D CNN LOO [24] | N/A | N/A | $46.6\pm 31.0$
Our 1D CNN | $89.53\pm 0.04$ | $74.53\pm 10.56$ | $89.22\pm 2.36$
Our 1D CNN-FG | $90.23\pm 0.05$ | $76.73\pm 11.44$ | $90.55\pm 3.84$
Our 1D CNN-GRU | $90.56\pm 0.03$ | $76.77\pm 8.05$ | $90.42\pm 3.15$
Figure 3: Area under ROC curve for all architectures.

###  IV-B Numerical Results

We consider 5 models for comparison. The baselines models from [16, 24], the
1D CNN we proposed, the CNN-aided factor graph, and the 1D CNN with GRU for
capturing the temporal correlations. From the baseline models, we implement
the exact architecture in [16] and evaluate it using the 6 fold
leave-4-patient-out. For [24], since it uses images as features, we just
report the results of their leave-one-patient-out (LOO) evaluations from their
paper.

Fig. 3 shows the average AUC-ROC across the 6 folds, while Table II summarizes
all the results. As can be seen, changing the model from 2D CNN proposed in
[16] to our proposed 1D CNN architecture results in approximately 2%
improvement across all performance measure. Our hybrid model-based/data-driven
CNN-aided factor graph further improves the results by as much as 2%. Compared
to [24], we perform a leave-4-patient-out evaluation as opposed to leave-one-
patient-out. This reduces the number of patients we use for training compared
to [24]. Despite this reduction, we achieve a much higher F1 score compared to
[24] as shown in Table II. Moreover, while AUC-PR is not reported in [24], the
precision of $51.4\pm 34.1$ and recall of $53.1\pm 25.5$ is reported, which is
significantly lower than the performance achieved by our proposed approach.

The 1D CNN with GRU achieves the same performance as our proposed CNN-aided
factor graph inference method for all performance metrics. This demonstrates
that the hybrid model-based/data-driven approach proposed here can achieve
similar performance as a purely data-driven methods that employ deep and
highly parameterized neural networks. This performance is achieved by our CNN-
aided factor graph at a fraction of the computational complexity of deep
learning based approaches as summarized in Table I. This makes the proposed
method suitable for real-time seizure detection.

##  V Conclusions

In this paper, we proposed a computationally efficient hybrid model-
based/data-driven method using CNN-aided factor graphs for seizure detection.
First, we carefully designed a 1D CNN for estimating the probability that a
4-second block of EEG recording is a seizure block. We then used this neural
network in the factor node of the factor graph for inference. We demonstrated
that the proposed method generalizes well to other patients using a 6-fold
leave-4-patient-out cross validation. We also showed that our algorithm
achieves up to 5% improvement in performance compared to prior work, while
maintaining much lower computational complexity. This makes our approach ideal
for real-time seizure detection. For future work, we plan to expand our
approach to classifying focal and generalized seizures since this is a
challenging task in clinical procedures.

## References

  * [1] “Epilepsy.” https://www.who.int/news-room/fact-sheets/detail/epilepsy.
  * [2] C. Park, G. Choi, J. Kim, S. Kim, T.-J. Kim, K. Min, K.-Y. Jung, and J. Chong, “Epileptic seizure detection for multi-channel EEG with deep convolutional neural network,” in International Conference on Electronics, Information, and Communication (ICEIC), 2018.
  * [3] G. Zazzaro, S. Cuomo, A. Martone, R. V. Montaquila, G. Toraldo, and L. Pavone, “EEG signal analysis for epileptic seizures detection by applying data mining techniques,” Internet of Things, vol. 14, p. 100048, 2021.
  * [4] S. Kulaseharan, A. Aminpour, M. Ebrahimi, and E. Widjaja, “Identifying lesions in paediatric epilepsy using morphometric and textural analysis of magnetic resonance images,” NeuroImage: Clinical, vol. 21, 2019.
  * [5] N. van Klink, A. Mooij, G. Huiskamp, C. Ferrier, K. Braun, A. Hillebrand, and M. Zijlmans, “Simultaneous MEG and EEG to detect ripples in people with focal epilepsy,” Clinical Neurophysiology, vol. 130, no. 7, pp. 1175–1183, 2019.
  * [6] N. Pianou and S. Chatziioannou, “Imaging with PET/CT in Patients with Epilepsy,” in Epilepsy Surgery and Intrinsic Brain Tumor Surgery, pp. 45–50, Springer, 2019.
  * [7] A. Subasi, J. Kevric, and M. Abdullah Canbaz, “Epileptic seizure detection using hybrid machine learning methods,” Neural Computing and Applications, vol. 31, pp. 317–325, Jan. 2019.
  * [8] M. Golmohammadi, S. Ziyabari, V. Shah, S. L. de Diego, I. Obeid, and J. Picone, “Deep architectures for automated seizure detection in scalp EEGs,” arXiv preprint arXiv:1712.09776, 2017.
  * [9] I. B. Slimen, L. Boubchir, Z. Mbarki, and H. Seddik, “EEG epileptic seizure detection and classification based on dual-tree complex wavelet transform and machine learning algorithms,” Journal of biomedical research, vol. 34, no. 3, p. 151, 2020.
  * [10] S. R. R. Ahmad, S. M. Sayeed, Z. Ahmed, N. M. Siddique, and M. Z. Parvez, “Prediction of epileptic seizures using support vector machine and regularization,” in 2020 IEEE Region 10 Symposium (TENSYMP), pp. 1217–1220, IEEE, 2020.
  * [11] S. Raghu, N. Sriraam, S. V. Rao, A. S. Hegde, and P. L. Kubben, “Automated detection of epileptic seizures using successive decomposition index and support vector machine classifier in long-term EEG,” Neural Computing and Applications, vol. 32, no. 13, pp. 8965–8984, 2020.
  * [12] C. Li, W. Zhou, G. Liu, Y. Zhang, M. Geng, Z. Liu, S. Wang, and W. Shang, “Seizure Onset Detection Using Empirical Mode Decomposition and Common Spatial Pattern,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 29, pp. 458–467, 2021.
  * [13] S. Khalilpour, A. Ranjbar, M. B. Menhaj, and A. Sandooghdar, “Application of 1-D CNN to predict epileptic seizures using EEG records,” in 2020 6th International Conference on Web Research (ICWR), pp. 314–318, IEEE, 2020\.
  * [14] G. C. Jana, R. Sharma, and A. Agrawal, “A 1D-CNN-spectrogram based approach for seizure detection from EEG signal,” Procedia Computer Science, vol. 167, pp. 403–412, 2020.  Publisher: Elsevier.
  * [15] R. V. Sharan and S. Berkovsky, “Epileptic seizure detection using multi-channel EEG wavelet power spectra and 1-D convolutional neural networks,” in 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), pp. 545–548, IEEE, 2020.
  * [16] P. Boonyakitanont, A. Lek-uthai, K. Chomtho, and J. Songsiri, “A Comparison of Deep Neural Networks for Seizure Detection in EEG Signals,” bioRxiv, p. 702654, 2019.
  * [17] S. Roy, I. Kiral-Kornek, and S. Harrer, “Deep learning enabled automatic abnormal EEG identification,” in International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 2756–2759, 2018.
  * [18] W. Liang, H. Pei, Q. Cai, and Y. Wang, “Scalp EEG epileptogenic zone recognition and localization based on long-term recurrent convolutional network,” Neurocomputing, vol. 396, pp. 569–576, 2020.
  * [19] M. Lee, I. Youn, J. Ryu, and D.-H. Kim, “Classification of both seizure and non-seizure based on EEG signals using hidden Markov model,” in 2018 IEEE International Conference on Big Data and Smart Computing (BigComp), pp. 469–474, IEEE, 2018.
  * [20] H.-A. Loeliger, “An introduction to factor graphs,” IEEE Signal Process. Mag., vol. 21, no. 1, pp. 28–41, 2004.
  * [21] P. D. Emmady and A. C. Anilkumar, “EEG, Abnormal Waveforms,” StatPearls [Internet], 2020.
  * [22] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K. Peng, and H. E. Stanley, “PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals,” circulation, vol. 101, no. 23, pp. e215–e220, 2000.
  * [23] W. Klonowski, “Everything you wanted to ask about eeg but were afraid to get the right answer,” Nonlinear biomedical physics, vol. 3, no. 1, pp. 1–5, 2009.
  * [24] C. Gomez, P. Arbelaez, M. Navarrete, C. Alvarado-Rojas, M. Le Van Quyen, and M. Valderrama, “Automatic seizure detection based on imaged-EEG signals through fully convolutional networks,” Scientific Reports, vol. 10, p. 21833, Dec. 2020.
  * [25] K.-O. Cho and H.-J. Jang, “Comparison of different input modalities and network structures for deep learning-based seizure detection,” Scientific Reports, vol. 10, p. 122, Jan. 2020.
  * [26] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” IEEE Trans. Inf. Theory, vol. 47, no. 2, pp. 498–519, 2001.
  * [27] G. D. Forney, “Codes on graphs: Normal realizations,” IEEE Trans. Inf. Theory, vol. 47, pp. 520–548, 2001.
  * [28] F. Mormann, T. Kreuz, R. G. Andrzejak, P. David, K. Lehnertz, and C. E. Elger, “Epileptic seizures are preceded by a decrease in synchronization,” Epilepsy research, vol. 53, no. 3, pp. 173–185, 2003.
  * [29] N. Shlezinger, N. Farsad, Y. C. Eldar, and A. J. Goldsmith, “Learned factor graphs for inference from stationary time sequences,” arXiv preprint arXiv:2006.03258, 2020.
  * [30] N. Shlezinger, N. Farsad, Y. C. Eldar, and A. J. Goldsmith, “ViterbiNet: A deep learning based Viterbi algorithm for symbol detection,” IEEE Trans. Wireless Commun., vol. 19, no. 5, pp. 3319–3331, 2020.
  * [31] D. M. Powers, “Evaluation: from precision, recall and F-measure to ROC, informedness, markedness and correlation,” arXiv preprint arXiv:2010.16061, 2020.

  *[EEG]: electroencephalogram
  *[DWT]: discrete wavelet transform
  *[DL]: deep learning
  *[CNN]: convolutional neural network
  *[RNN]: recurrent neural network
  *[HMM]: hidden Markov model
  *[HMMs]: hidden Markov model
  *[CNNs]: convolutional neural network
  *[FLOPs]: floating point operations


# Granular Gases under Extreme Driving

W. Kang<EMAIL_ADDRESS>Department of Physics, University of Massachusetts,
Amherst, Massachusetts 01003 USA  J. Machta<EMAIL_ADDRESS>Department of Physics, University of Massachusetts, Amherst, Massachusetts
01003 USA  E. Ben-Naim<EMAIL_ADDRESS>Theoretical Division and Center for
Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico
87545 USA

###### Abstract

We study inelastic gases in two dimensions using event-driven molecular
dynamics simulations. Our focus is the nature of the stationary state attained
by rare injection of large amounts of energy to balance the dissipation due to
collisions. We find that under such extreme driving, with the injection rate
much smaller than the collision rate, the velocity distribution has a power-
law high energy tail. The numerically measured exponent characterizing this
tail is in excellent agreement with predictions of kinetic theory over a wide
range of system parameters.  We conclude that driving by rare but powerful
energy injection leads to a well-mixed gas and constitutes an alternative
mechanism for agitating granular matter. In this distinct nonequilibrium
steady-state, energy cascades from large to small scales. Our simulations also
show that when the injection rate is comparable with the collision rate, the
velocity distribution has a stretched exponential tail.

###### pacs:

45.70.Mg, 47.70.Nd, 05.40.-a, 81.05.Rm

Granular materials are ubiquitous in nature, but nevertheless, fundamental
understanding of the properties of granular materials presents many challenges
jaeger96 ; kadanoff99 ; gennes99 ; ig ; bp . Underlying these challenges are
structural inhomogeneities, macroscopic particle size, and energy dissipation,
all of which are defining features of granular matter.

Equilibrium gases have Maxwellian velocity distributions. Due to the
irreversible nature of the dissipative collisions, granular gases are out of
equilibrium. Indeed, non-Maxwellian velocity distributions are observed in a
wide range of experiments in driven granular matter including in particular
shaken grains olafsen98 ; gzb ; losert99 ; cafiero2000 ; rm ; blair01 ;
aranson02 ; prevost02 ; wp ; fm ; tmhs ; kohlstedt05 . In such experiments,
energy is injected over a wide range of scales and the measured velocity
distribution has a stretched exponential form. To a large extent, a kinetic
theory where energy injection through the system boundary is modeled by a
thermostat successfully describes these nonequilibrium steady-states ep ;
brey03 ; ve ; krb .

Furthermore, theoretical studies suggest that the steady-state is controlled
primarily by the ratio between the energy injection rate and the collision
rate mackintosh04 ; mackintosh05 . When the injection rate is much larger than
the collision rate, the velocity distribution is Maxwellian. However, when the
injection rate is smaller than the collision rate, the velocity distribution
is non-Maxwellian, and has a stretched-exponential tail.

In this study, we focus on the limiting case where the injection rate is
vanishingly small and energy is injected at extremely large velocity scales bm
; bmm . Under such extreme driving, the injected energy cascades down from
large velocity scales to small scales and thereby counters the dissipation by
collisions. Kinetic theory shows that in the stationary state, the velocity
distribution has a power-law high energy tail. These theoretical predictions
were supplemented by Monte-Carlo simulations of the homogeneous Boltzmann
equation where spatial correlations are ignored. However, such nonequilibrium
steady-states have yet to be observed using more realistic molecular dynamics
simulations. In this work we carry out extensive molecular dynamics
simulations to investigate the behavior of inelastic hard disks under extreme
driving.

The goal of this investigation is to establish whether extreme driving is a
feasible mechanism for driving granular matter. Our main result is that under
rare but powerful injection of energy, a granular gas indeed reaches a
stationary state that is characterized by a power-law velocity distribution.
Moreover, our simulations quantitatively confirm the predictions of the
kinetic theory as the exponent characterizing the tail of the distribution is
validated over a wide range of parameters. Our results show that extreme
driving is a feasible mechanism for agitating granular matter.

Kinetic Theory. Our starting point is the observation that the purely-
collisional homogeneous Boltzmann equation supports stationary solutions bm ;
bmm ; etb ; bz . The evolution equation for the velocity distribution $f({\bf
v})$ of inelastic hard disks takes the form

$\displaystyle\frac{\partial f({\bf v})}{\partial t}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\iiint d\hat{\bf n}\,d{\bf u}_{1}\,d{\bf u}_{2}\left|({\bf
u}_{1}-{\bf u}_{2})\cdot\hat{\bf n}\right|f({\bf u}_{1})f({\bf u}_{2})$ (1)
$\displaystyle\times$ $\displaystyle[\delta({\bf v}-{\bf v}_{1})-\delta({\bf
v}-{\bf u}_{1})],$

with $\hat{\bf n}$ the impact direction. This Boltzmann equation is
supplemented by the inelastic collision rule which specifies the post-
collision velocities ${\bf v}_{1,2}$ as a linear combination of the pre-
collision velocities ${\bf u}_{1,2}$,

${\bf v}_{1,2}={\bf u}_{1,2}-\frac{1+r}{2}({\bf u}_{1,2}-{\bf
u}_{2,1})\cdot\hat{\bf n}\,\hat{\bf n}.$ (2)

In an inelastic collision, the normal component of the relative velocity
reverses sign and is scaled down by the restitution coefficient $0\leq r\leq
1$, $({\bf v}_{1}-{\bf v}_{2})\cdot\hat{\bf n}=-r({\bf u}_{1}-{\bf
u}_{2})\cdot\hat{\bf n}$. The energy loss equals $\Delta
E=-\frac{1-r^{2}}{4}|({\bf u}_{1}-{\bf u}_{2})\cdot\hat{\bf n}|^{2}$.

The collision rule (2) simplifies to a “fragmentation” rule ${\bf u}\to({\bf
w}_{1},{\bf w}_{2})$ for collisions involving one extremely energetic particle
with velocity ${\bf u}$ and a second, implicit, particle with speed much less
than $|{\bf u}|$. The post-collision velocities ${\bf w}_{1}=\frac{1+r}{2}{\bf
u}\cdot\hat{\bf n}\,\hat{\bf n}$ and ${\bf w}_{2}={\bf u}-\frac{1+r}{2}{\bf
u}\cdot\hat{\bf n}\,\hat{\bf n}$ follow by substituting ${\bf u}_{1}=0$ and
${\bf u}_{2}={\bf u}$, respectively, into (2). In the limit $|{\bf
v}|\to\infty$, the nonlinear Boltzmann equation (1) becomes linear and its
stationary form is

$\displaystyle 0\\!=\\!\iint d\hat{\bf n}\,d{\bf u}\left|{\bf u}\cdot\hat{\bf
n}\right|f({\bf u})\left[\delta({\bf v}\\!-\\!{\bf w}_{1})\\!+\\!\delta({\bf
v}\\!-\\!{\bf w}_{2})\\!-\\!\delta({\bf v}\\!-\\!{\bf u})\right].$ (3)

For arbitrary dimension and for arbitrary collision parameters, this linear
and homogeneous equation admits the power-law solution

$f(v)\sim v^{-\sigma}.$ (4)

In two-dimensions, the subject of our investigation, the exponent $\sigma$
obeys the transcendental equation bm

$\frac{1-\\!{}_{2}F_{1}\left(\frac{3-\sigma}{2},1,\frac{3}{2},1-(\frac{1-r}{2})^{2}\right)}{\left(\frac{1+r}{2}\right)^{\sigma-3}}=\frac{\Gamma(\frac{\sigma-1}{2})\Gamma(\frac{3}{2})}{\Gamma(\frac{\sigma}{2})},$
(5)

where ${}_{2}F_{1}$ is the hypergeometric function. The exponent $\sigma$
grows monotonically with the restitution coefficient. The limiting values are
$\sigma=4.14922$ for completely inelastic collisions ($r=0$) and $\sigma=5$ in
the elastic limit ($r\rightarrow 1$).

The power-law distribution (4) is a stationary solution of the linear
Boltzmann equation (3). Yet, Monte Carlo methods show that the full nonlinear
Boltzmann equation (1) does admit a stationary solution with a tail given by
(4). These numerical solutions are computed by injecting energy at a rate that
is much smaller than the collision rate. In an individual injection event a
randomly chosen particle is given a velocity much larger than the typical
velocity. Such extreme driving maintains a steady-state in which energy
injection balances energy dissipation.

The physical mechanism underlying these driven steady-states is a cascade in
which a high energy particle collides with a particle of typical energy
yielding two high energy particles, each with energies less than that of the
original high energy particle. These two high energy particles produce two
more high energy particles, again by collisions with the much more numerous
particles with typical energies, and so on.

Scaling Analysis. Interestingly, there is a family of steady-states generated
by extreme driving. If $f(v)$ is a stationary solution of (1), then
$v_{0}^{-2}f(v/v_{0})$ with the arbitrary typical velocity $v_{0}$ is also a
stationary solution because the collision rule (2) is invariant under the
scale transformation ${\bf v}\to{\bf v}/v_{0}$. The energy injection rate
$\gamma$, the velocity injection scale $V$, and the typical velocity $v_{0}$
are related by the energy balance requirement.

We relate these three quantities by a heuristic argument and first note that
the energy injection rate is simply $\gamma V^{2}$. Also, we anticipate that
the velocity distribution is truncated at the injection scale $V$. The energy
dissipation is dominated by the tail of the distribution and is controlled by
the upper cut-off $V$. We estimate the dissipation rate $\Gamma$ as follows
norm ,

$\displaystyle\Gamma\sim\rho\\!\int^{V}\\!\\!\\!v\cdot
v^{2}\frac{1}{v_{0}^{2}}f\left(\frac{v}{v_{0}}\right)\,v\,dv\sim\rho
V^{3}(V/v_{0})^{2-\sigma},$ (6)

where $\rho$ is the particle density. In the integrand, the first term $v$
accounts for the collision rate, and the second term $v^{2}$ accounts for the
energy dissipation in an inelastic collision. The integration is performed
using the velocity distribution (4), and since $\sigma<5$ the dissipation is
indeed dominated by the high velocity tail of the distribution. Balancing
energy injection with dissipation we obtain a relationship between the
injection rate $\gamma$, the injection scale $V$, and the typical velocity
$v_{0}$,

$\gamma\sim\rho\,V(V/v_{0})^{2-\sigma}.$ (7)

Since the collision rate is proportional to $\rho\,v_{0}$, the dimensionless
ratio $\psi$ of the injection rate to the collision rate scales as a power of
the velocity ratio $V/v_{0}$,

$\psi\sim(V/v_{0})^{3-\sigma}.$ (8)

Since $\sigma>3$, we expect a wide power-law range, $V\gg v_{0}$, when the
injection rate is much smaller than the collision rate, $\psi\ll 1$.  There is
no lower cutoff on the injection rate $\psi$, below which power-law is not
observed; the smaller is $\psi$, the broader the power-law range. When $\psi$
is order one, the velocity distribution no longer has a power-law tail, and of
course, when $\psi\gg 1$, the velocity distribution should simply mirror the
distribution of injected velocities.

Molecular Dynamics Simulations. We used molecular dynamics alder59  to
simulate inelastic hard disks in a square box with elastic walls. In these
event driven simulations, upon impact, the velocities of the colliding
particles are updated according to the collision rule (2). Subsequent to each
collision, we identify the time and location of the next collision. The
particles undergo purely ballistic motion between two successive collisions.

We implemented the following velocity-dependent restitution coefficient
goldman98

$r(\delta_{n})=\begin{cases}1-(1-r)(\delta_{n}/v_{c})^{3/4}&\text{$\delta_{n}<v_{c}$},\\\
r&\text{$\delta_{n}\geq v_{c}$},\end{cases}$ (9)

where $\delta_{n}=(\mathbf{v_{1}}-\mathbf{v_{2}})\cdot\mathbf{\hat{n}}$ is the
normal component of the relative velocity. Here, $r$ is the nominal value of
the restitution coefficient, valid at large velocities, and $v_{c}$ is the
cutoff velocity, below which collisions become elastic. With this realistic
restitution coefficient luding96 ; bizon98 , we avoid inelastic collapse where
an infinite number of collisions can occur in a finite time my . Typically, we
set $v_{c}$ much smaller than the typical velocity $v_{0}$, but for small
restitution coefficients, we must set $v_{c}$ comparable to $v_{0}$ to avoid
inelastic collapse.

To maintain a steady-state, we periodically boost a single randomly-selected
particle to a large, random velocity. These injection events are rare and they
are governed by a Poisson process with rate $\gamma$, that is, with
probability $\gamma\,dt$ injection is implemented during the time interval
$[t,t+dt]$. The injection speed is selected from a Gaussian distribution with
zero mean and standard deviation $V$. By taking a long-time average, we
confirmed that the total energy approaches a constant, and hence, that the
system reaches a statistical steady-state where energy injection and energy
dissipation balance. Moreover, the velocity distributions were produced by
sampling particle velocities at a very large number ($10^{8}$) of equally
spaced time intervals.  We stress that the velocities are sampled at time
intervals that are completely uncorrelated with either collision events or
injection events. We tested that our sampling produces robust velocity
distributions, and that the velocity distribution, representing an average
over the entire system, is truly stationary. In particular, the system does
not enter the homogeneous cooling state pkh in between the rare injection
events.

Figure 1: The velocity distribution $f(v)$ versus the velocity $v$. Molecular
dynamics simulation results (bullets) are compared with the power-law tail
predicted by kinetic theory (solid line).

We performed numerical simulations using a system of $N=10^{3}$ identical
particles with diameter $2R=1$ in a square box of size $L=400$, corresponding
to the low area fraction $\phi=N\pi\,(R/L)^{2}=4.9\times 10^{-3}$. Unless
noted otherwise, these parameters are used throughout this study. Energy was
injected at rate $\gamma=5\times 10^{-7}$ and the injection scale was $V=850$.
First, we considered weakly-inelastic particles, $r=0.9$, with the cutoff
$v_{c}=0.1$. With these parameters, the injection rate is much smaller than
the measured collision rate and their ratio is $\psi=1.5\times 10^{-5}$.
Consequently, there is substantial scale separation between the typical
velocity $v_{0}$ and the injection velocity $V$. Over this range, the steady-
state velocity distribution obtained by the molecular dynamics simulations has
a power-law high energy tail as in (4) and the exponent $\sigma=4.74$ is in
very good agreement with the kinetic theory prediction given by (5),
$\sigma=4.74104$ (see figure 1). We also confirmed that the ratio
$V/v_{0}\approx 10^{2}$ is consistent with the scaling estimate (8).

Next, we varied the restitution coefficient and repeated the simulations. Over
the entire parameter range $0.1\leq r\leq 0.9$, we find stationary velocity
distributions with a power-law tail. In general, the ratio $V/v_{0}$ is
consistent with the scaling relation (8). Moreover, the exponent $\sigma$
obtained from the molecular dynamics simulations is in excellent agreement
with the kinetic theory predictions (5) for all restitution coefficients
(figure 2). We thus arrive at our main result that at least for dilute gases,
extreme driving in the form of rare but powerful energy injection generates a
steady-state with a broad distribution of velocities. The tail of the velocity
distribution is power-law and the characteristic exponent is nonuniversal as
it depends on the restitution coefficient.

Figure 2: The exponent $\sigma$ versus the restitution coefficient $r$. The
molecular dynamics results (bullets) are compared with the kinetic theory
predictions (solid line).

The fact that kinetic theory holds shows that, to good approximation, the gas
is well-mixed. We comment that it is remarkable that extreme driving results
in a well-mixed gas. On short time scales, energy injection clearly generates
spatial correlations because a just-energized particle transfers much of its
energy to nearby particles by inelastic collisions. Yet, on larger time scales
energetic particles break coherent structures which are known to be the
consequence of inelastic collisions gz . While these two mechanisms have
opposite effects, the simulations indicate that when a long time average is
taken, the latter effect dominates. Thus, energy injected at extreme velocity
scales in a tightly localized region of space, ends up evenly distributed
throughout the system.

Snapshots of the time evolution of the system following energy injection
demonstrate how the inelastic cascade works (figure 3). In the initial stages
of the cascade, the injection affects only a small region in space and,
moreover, there are strong spatial correlations between the velocities of the
particles (figure 3 a-c). However, after many inelastic collisions, the
injected energy ends up evenly distributed throughout the system (figure 3 d).
When a long time average is taken over many energy injection events at
different locations, the system is maintained in a homogeneous, well-mixed
state. Spatial correlations induced by inelasticity and the injection
mechanism do not affect the predicted power-law velocity distributions.

(a) (b)

(c) (d)

Figure 3:  The inelastic energy cascade. Shown are four time-ordered snapshots
of the gas shortly after an injection event. The top three figures show a
small window around the injection event in the early stages of the cascade:
(a) the initial energetic particle (red online) (b) two energetic particles
after one collision, (c) four energetic particles after three collisions.
Figure (d) shows the entire system, with energetic particles shown larger (red
online), in a late stage of the cascade. After many collisions, the injected
energy is evenly distributed throughout the system. The simulation parameters
are $r=0.8$, $N=1000$, $L=400$, $V=707$ and $\gamma=2\times 10^{-6}$.

Therefore, rare, powerful, and spatially localized energy injection is a
unique mechanism of agitating granular gases. This mechanism induces an
extended energy cascade which distributes the injected energy to the rest of
the system. This physical mechanism is different than energy injection by
walls gzb ; losert99  or by an effective thermostat ve  or by multiplicative
driving cafiero2000 where the injected energy directly affects only a small
region of space. In this sense, the energy cascade represents a novel
mechanism for agitating granular matter.

Figure 3 also illustrates that the velocity distribution is correlated with
injection times for our system so that the predicted power-law distribution
arises only after averaging over many measurements taken at times that are
uncorrelated with the injection times. However, for very large systems driven
with a fixed but small injection rate per particle there would be many
temporally overlapping but spatially well-separated injection events. At any
instant of time in a very large system cascades at all stages of development
would be present somewhere in the system and the power law tail would be time-
independent.

We performed additional simulations to test whether the results are robust
with respect to change of parameters. In particular, we varied the area
fraction by fixing the number of particles and varying the system size. The
results shown in figure 4 are for three different area fractions:
$\phi=7.7\times 10^{-5}$, $\phi=4.9\times 10^{-3}$, and $\phi=7.9\times
10^{-2}$.  The corresponding values of $\psi$ are $4.7\times 10^{-4}$,
$5.4\times 10^{-5}$, and $3.4\times 10^{-6}$, respectively. Note that in all
cases $\psi\ll 1$. We find the same power-law tail in all three cases, and the
exponent is in good quantitative agreement with the kinetic theory prediction
(figure 4).  Thus, the energy cascade mechanism can be realized even at area
fractions as high as $\phi\approx 10^{-1}$ and with $\psi$ as small as order
$10^{-6}$.

Figure 4: The velocity distribution at three different densities (solid
lines). The simulation parameters are: $r=0.8$, $\gamma=2\times 10^{-6}$, and
$V=707$. Eq. (4) with $\sigma=4.57246$ is also shown as a reference (dashed
line). A best fit to a power-law yields $\sigma=4.6$, $4.6$, and $4.7$ for
$\phi=7.7\times 10^{-5}$, $4.9\times 10^{-3}$, and $7.9\times 10^{-2}$,
respectively.

We also studied the dependence on the ratio $\psi$ between the injection rate
and the collision rate by varying the injection rate $\gamma$ and the
injection velocity $V$. In accord with (8), we find that the range $[v_{0},V]$
of power-law behavior shrinks as $\psi$ increases. As long as $\psi$ is
sufficiently small, the distribution has a power-law tail (figure 5). When
this ratio becomes sufficiently large, the tail is no longer algebraic and a
sharper decay occurs. For $\psi=5.8\times 10^{-2}$, we find the stretched
exponential tail $f(v_{x})\sim\exp\big{(}-{\rm
const.}\times|v_{x}|^{\zeta}\big{)}$ with $\zeta=1.52$ (figure 6). This value
is consistent with the theoretical value $\zeta=3/2$ for inelastic gases
driven by white noise ve ; krb  and the experimental value observed in
vigorously shaken beads rm .  Indeed, in this intermediate injection rate
regime, the energy cascade becomes localized, and frequent, small injections
are similar to white noise driving. On the other hand, these stretched
exponential tails do not relate to those observed in bmm  as the system cools
down after injection is turned-off. We find stretched exponential tails for
the range $10^{-1}\lessapprox\psi\lessapprox 1$. When the injection rate
exceeds the collision rate ($\psi\gg 1$) the entire distribution becomes
Maxwellian as the velocity distribution simply mirrors the distribution of
injected velocities.

Figure 5: The velocity distribution $f(w)$ versus the normalized velocity
variable $w=v/v_{0}$ with different ratios of injection rate to collision
rate: $\psi=1.4\times 10^{-3}$ (bullets), $2.7\times 10^{-3}$ (squares),
$1.1\times 10^{-2}$ (diamonds), $2.3\times 10^{-2}$ (up-triangles), and
$4.4\times 10^{-2}$ (down-triangles). The restitution coefficient is $r=0.8$.
Also shown is the reference theoretical curve (4) with $\sigma=4.57246$
(dashed line).

Finally, we mention that we even varied the energy injection mechanism itself.
In particular, to implement injection strictly at large energy scales, we used
an “energy loss counter” to keep track of the total energy dissipated by
collisions since the last injection. When the dissipated energy equals a fixed
large value, we inject this amount of energy into a single randomly chosen
particle bm ; bmm . Using this variant of extreme driving to maintain the
steady-state, we found similar power-law velocity distributions.  From these
studies, we conclude that the parameter $\psi$ controls the velocity
distribution, and that different energy injection mechanisms lead to a power-
law distribution with the very same exponent $\sigma$, as long as $\psi\ll 1$.

Conclusions. Extensive event-driven simulations show that under extreme
driving in the form of rare but powerful energy injections, an inelastic gas
reaches a steady-state with a broad distribution of energies. Such driven
steady-states are observed for a wide range of collision parameters,
densities, and energy injection rates, as long as the injection rate is much
smaller than the collision rate. The velocity distributions have a power-law
tail and the characteristic exponent is in good agreement with the kinetic
theory predictions. When the ratio between the energy injection rate and the
collision rate becomes sufficiently large, the velocity distribution has a
much sharper stretched exponential tail.

Figure 6: The distribution $f(v_{x})$ of the horizontal component of the
velocity $v_{x}$ at a moderate $\psi=5.8\times 10^{-2}$. These simulations are
performed with $r=0.8$, $\gamma=6\times 10^{-4}$, and $V=1.41$. The solid line
is a best-fit to the stretched exponential $f(v_{x})\sim\exp\big{(}-{\rm
const.}\times|v_{x}|^{\zeta}\big{)}$ with $\zeta=1.52$.

We conclude that extreme driving where energy is injected only at very large
scales presents an alternative mechanism for agitating granular matter and
that such driving leads to a fundamentally different steady-state compared
with traditional driving where energy is injected over all scales. Realizing
this driving in experiments is a challenge because the agitation must be
applied only at very large velocities. One possible mechanism is shooting very
fast particles into the system. While such a system would involve a growing
number of particles, the injection of energetic particles should lead to
transfer of energy from large scales to small scales by a cascade of
collisions.

Figure 7:  The time development, left to right, of a cascade in a denser
system ($\phi=7.9\times 10^{-2}$ and the same parameters as in Fig. 4.)
Particles with speeds greater than a fixed threshold are shown larger (red
online).

Interestingly, the power-law velocity distribution appears to hold in systems
that include dense clusters. We have observed that when the density is
increased, there is a tendency for clustering near the walls with no
measurable deviation from the predicted power law (figure 7). In this case
injection leads to explosive breakup of dense regions. It is intriguing that
the power-law tail is quite robust, and extends to situations where the
assumptions underlying the kinetic theory approach can no longer be justified.

In our simulations, inelastic collapse does not play a role because collisions
become elastic at small relative velocities. Yet, if the collisions are purely
inelastic, there should be a competition between the formation of high-density
regions by inelastic collisions and the destruction of such clusters by high
energy particles. Elucidating this competition is another possibility for
further investigation. Nonetheless, our simulations suggest that extreme
driving generates well-mixed steady-states despite the fact that this driving
is very inhomogeneous in both space and time.

Acknowledgments. We thank Narayanan Menon, Felix Werner for useful discussions
and Hong-Qiang Wang for sharing the molecular dynamics code. We acknowledge
Financial support from NSF grant DMR-0907235 (JM and WK) and DOE grant DE-
AC52-06NA25396 (EB).

## References

  * (1) H. M. Jaeger, S. R. Nagel and R. P. Behringer, Rev. Mod. Phys. 68, 1259 (1996).
  * (2) L. P. Kadanoff Rev. Mod. Phys. 71, 435 (1999).
  * (3) P. G. de. Gennes, Rev. Mod. Phys. 71, S374 (1999).
  * (4) I. Goldhirsch, Ann. Rev. Fluid. Mech. 35, 267 (2003).
  * (5) Kinetic theory of granular gases, N. Brilliantov and T. Pöschel, (Oxford, Oxford, 2003).
  * (6) J. S. Olafsen and J. S. Urbach, Phys. Rev. Lett. 81, 4369 (1998); Phys. Rev. E 60, R2468 (1999).
  * (7) E. L. Grossman, T. Zhou, and E. Ben-Naim, Phys. Rev. E 55, 4200 (1997).
  * (8) W. Losert, D. G. W. Cooper, J. Delour, A. Kudrolli and J. P. Gollub, Chaos 9, 682, 1999.
  * (9) R. Cafiero, S. Luding, and H. J. Herrmann, Phys. Rev. Lett. 84, 6014 (2000).
  * (10) F. Rouyer and N. Menon, Phys. Rev. Lett. 85, 3676 (2000).
  * (11) D. L. Blair and A. Kudrolli, Phys. Rev. E 64, R050301 (2001).
  * (12) I. S. Aranson, and J. S. Olafson, Phys. Rev. E 66, 061302 (2002).
  * (13) K. Kohlstedt, A. Snezhko, M. V. Sakozhnikov, I. S. Aranson, J. S. Olafson, E. Ben-Naim, Phys. Rev. Lett. 95, 068001 (2005).
  * (14) A. Prevost and et al, Phys. Rev. Lett. 89, 084301 (2002).
  * (15) R. D. Wildman and D. J. Parker, Phys. Rev. Lett. 88, 064301 (2002).
  * (16) K. Feitosa and N. Menon, Phys. Rev. Lett. 88, 198301 (2002).
  * (17) S. Tatsumi, Y. Murayama, H. Hayakawa, M. Sanoi, J. Fluid Mech. 641, 521 (2009).
  * (18) S. E. Esipov and T. Pöschel, J. Stat. Phys. 86, 1385 (1997).
  * (19) J. J. Brey and M. J. Ruiz-Montero, Phys. Rev. E 67, 021307 (2003).
  * (20) T. P. C. van Noije and M. H. Ernst, Gran. Matt. 1, 57 (1998).
  * (21) P. L. Krapivsky, S. Redner, and E. Ben-Naim, A Kinetic View of Statistical Physics (Cambridge University Press, Cambridge, 2010).
  * (22) J. S. van Zon, and F. C. MacKintosh, Phys. Rev. Lett. 93, 038001 (2004).
  * (23) J. S. van Zon, and F. C. MacKintosh, Phys. Rev. E 72, 051301 (2005).
  * (24) E. Ben-Naim and J. Machta, Phys. Rev. Lett. 94, 138001 (2005).
  * (25) E. Ben-Naim, B. Machta, and J. Machta, Phys. Rev. E 72, 021302 (2005).
  * (26) M. H. Ernst, E. Trizac, and A. Barrat, Europhys. Lett. 76, 56 (2006).
  * (27) E. Ben-Naim and A. Zippelius, J. Stat. Phys. 129 677 (2007).
  * (28) Throughout this paper we use the normalization $2\pi\int_{0}^{\infty}f(v)\,v\,dv=1$ where $v\equiv|{\bf v}|\equiv\sqrt{v_{x}^{2}+v_{y}^{2}}$.
  * (29) B. J. Alder and T. E. Wainwright, J. Chem. Phys. 31, 459 (1959).
  * (30) D. Goldman, M. D. Shattuck, C. Bizon, W. D. McCormick, J. B. Swift, and H. L. Swinney, Phys. Rev. E 57, 4831 (1998).
  * (31) C. Bizon, M. D. Shattuck, J. B. Swift, W. D. McCormick, and H. L. Swinney, Phys. Rev. Lett. 80, 57 (1998).
  * (32) S. Luding, E. Clement, J. Rajchenbach, and J. Duran, Europhys. Lett. 34, 247 (1996).
  * (33) S. McNamara and W. R. Young, Phys Fluids A 4, 496 (1992).
  * (34) P. K. Haff, J. Fluid. Mech. 134, 401 (1983).
  * (35) I. Goldhirsch, and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993).