# Predicting discrete-time bifurcations with deep learning

Thomas M. Bury  Daniel Dylewsky  Department of Applied Mathematics, University
of Waterloo, Waterloo, Canada  Chris T. Bauch  Department of Applied
Mathematics, University of Waterloo, Waterloo, Canada  Madhur Anand  School of
Environmental Sciences, University of Guelph, Guelph, Canada  Leon Glass
Department of Physiology, McGill University, 3655 Promenade Sir William Osler,
Montreal, Quebec H3G 1Y6, Canada  Alvin Shrier  Gil Bub

###### Abstract

Many natural and man-made systems are prone to critical transitions—abrupt and
potentially devastating changes in dynamics. Deep learning classifiers can
provide an early warning signal (EWS) for critical transitions by learning
generic features of bifurcations (dynamical instabilities) from large
simulated training data sets. So far, classifiers have only been trained to
predict continuous-time bifurcations, ignoring rich dynamics unique to
discrete-time bifurcations. Here, we train a deep learning classifier to
provide an EWS for the five local discrete-time bifurcations of codimension-1.
We test the classifier on simulation data from discrete-time models used in
physiology, economics and ecology, as well as experimental data of
spontaneously beating chick-heart aggregates that undergo a period-doubling
bifurcation. The classifier outperforms commonly used EWS under a wide range
of noise intensities and rates of approach to the bifurcation. It also
predicts the correct bifurcation in most cases, with particularly high
accuracy for the period-doubling, Neimark-Sacker and fold bifurcations. Deep
learning as a tool for bifurcation prediction is still in its nascence and has
the potential to transform the way we monitor systems for critical
transitions.

## Introduction

Many systems in nature and society possess critical thresholds at which the
system undergoes an abrupt and significant change in dynamics [1, 2]. In
physiology, the heart can spontaneously transition from a healthy to a
dangerous rhythm [3]; in economics, financial markets can form a ‘bubble’ and
crash into a recession [4]; and in ecology, ecosystems can collapse as a
result of their interplay with human behaviour [5, 6]. These events,
characterised by a sudden switch to a different dynamical regime, are referred
to as critical transitions.

Critical transitions can be better understood with bifurcation theory [7, 8],
a branch of mathematics that studies how dynamical systems can undergo sudden
qualitative changes as a parameter crosses a threshold (a bifurcation). Many
bifurcations are accompanied by critical slowing down—a diminishing of the
local stability of the system—which results in systematic changes to
properties of a noisy time series, such as its variance, autocorrelation and
power spectrum [9, 10, 11]. These properties can be approximated analytically
in the presence of different bifurcations [12, 13, 14, 10], and a
corresponding observation in data can be used as an early warning signal (EWS)
for the bifurcation [11]. Systematic changes in variance and lag-1
autocorrelation have been observed prior to transitions in climate [15, 16,
17], geological [18], ecological [19, 20] and cardiac [21] systems, suggesting
the presence of a bifurcation. However, these EWS have limited ability in
predicting the type of bifurcation [22, 14] and can fail in systems with
nonsmooth potentials [23] or noise-induced transitions [24].

More recently, deep learning techniques have been employed to provide EWS for
bifurcations [25, 26, 27]. This involves training a neural network to classify
a time series based on the type of bifurcation it is approaching, as well as
appropriate controls [25, 26, 27]. Unlike many applications of deep learning,
this approach does not require abundant data from the study system, which, in
the context of critical transitions, is often unavailable. (Unfortunately we
do not have data from thousands or more ecosystems or climate systems that
went through bifurcation.) Instead, the approach generates a massive library
of simulation data from generic models that possess each type of bifurcation.
The neural network then learns generic features associated with each type of
bifurcation, that can be recognised in an unseen time series of the study
system. This is enabled by the existence of universal properties of
bifurcations that are manifested in time series as a dynamical system gets
close to a bifurcation [7, 9]. In our previous work, we trained a deep
learning classifier to provide an EWS for continuous-time bifurcations, and
found it was effective at predicting transitions for real thermoacoustic,
climate and geological transitions [25].

Bifurcations can be partitioned according to whether they occur in continuous
or discrete-time dynamical systems [8, 7]. This distinction is important,
since discrete-time dynamical systems (difference equations) can display very
different behaviour to their continuous-time counterparts (differential
equations). As an example, consider the logistic model for population growth.
When set up in continuous time (appropriate for populations with overlapping
generations e.g. humans), the population grows smoothly as the reproduction
rate increases. Whereas, when set up in discrete time (appropriate for
populations with non-overlapping generations, e.g. insects), the population
undergoes a sequence of period-doubling bifurcations to chaos [28]. This
involves transitioning from equilibrium dynamics to an oscillation of period
2, to an oscillation of period 4, 8, 16 etc. until the dynamics become
chaotic. It is therefore important to develop EWS suitable for both continuous
and discrete-time bifurcations. While indicators like variance and lag-1
autocorrelation can provide EWS for discrete-time bifurcations, the ability of
a deep learning classifier at this task has not been investigated.

As well as in ecology, discrete-time bifurcations arise naturally in
physiology [3], epidemiology [29], and economics [30], where events can take
place on a discrete timeline. To illustrate our approach, we will use model
simulations from ecology, physiology and economics, as well as experimental
data from spontaneously beating chick heart aggregates [31, 21]. Following
administration of a drug, in some aggregates the time interval between two
heart beats begins to alternate i.e. there is a period doubling bifurcation
(Figure 1). Such transitions can also occur for the human heart in the form of
T-wave alternans, which increases a patient’s risk for sudden cardiac death
[32]. The period-doubling bifurcation is accompanied by critical slowing down,
so systematic changes in variance and lag-1 autocorrelation are expected and
have been shown to provide an EWS in this system [21]. The chick heart
aggregates serve as a good study system to test the performance of EWS since
we have multiple recordings, not all of which underwent a transition, allowing
us to test for false positives.

Figure 1: Period-doubling bifurcation in a spontaneously beating aggregate of
embryonic chick heart cells following treatment with a potassium channel
blocker (E-4031, 1.5$\upmu$mol). (a) Interbeat intervals (IBI) for consecutive
beats. A period-doubling bifurcation occurs at approximately beat 230. Arrows
correspond to traces plotted in lower panels. (b-d) Normalised signal from
optical imaging of the aggregate’s motion. Traces are from a section well
before (b), just before (c), and after (d) the period-doubling bifurcation.

Among discrete-time bifurcations, there are many types, each with an
associated change in dynamics [7]. For this study, we focus on the five local
bifurcations of codimension one (Box 1). In being ‘local’, these bifurcations
are accompanied by critical slowing down, so systematic changes and variance
and autocorrelation are expected. However, not all of these bifurcations
result in a critical transition [22]. They can instead involve a smooth
transition to an intersecting steady state (transcritical) or to oscillations
with gradually increasing amplitude (supercritical Neimark-Sacker). Predicting
the type of bifurcation provides information on the nature of the dynamics
following the bifurcation, something variance and autocorrelation alone do not
provide.

Here, we train a deep learning classifier to provide a specific EWS for
bifurcations of discrete-time dynamical systems. We train the classifier using
simulation data of normal form equations appended with higher-order terms and
noise. We then test the classifier on simulation runs of five discrete-time
models used in cardiology, ecology and economics, and assess its performance
relative to variance and lag-1 autocorrelation. We vary the noise amplitude
and rate of forcing in model simulations to assess robustness of the EWS.
Finally, we test the classifier on experimental data of spontaneously beating
chick-heart aggregates that go through a period-doubling bifurcation. A
reproducible run of all analyses may be performed on Code Ocean
(https://codeocean.com/capsule/2209652/tree/v1) where the code is accompanied
by the necessary software environment.

Box 1. Local discrete-time bifurcations of codimension one.
A discrete-time dynamical system has the form $\vec{x}_{t+1}=f(\vec{x}_{t})$
(1) where $\vec{x}_{t}$ is a vector representing the state of the system and
$f$ is a function that maps to the next state. The system is stable about a
steady state $\vec{x}_{*}$ if all eigenvalues ($\lambda$) of the Jacobian
matrix (gradient of $f$ at $\vec{x}_{*}$) have a magnitude less than one i.e.
they lie within the unit circle on the complex plane (depicted in gray). A
local bifurcation occurs when one or a pair of these eigenvalues (the dominant
eigenvalue(s)) crosses the unit circle. The type of bifurcation depends on
where the unit circle is crossed. A bifurcation is said to be of codimension
one if it occurs upon changing only a single parameter of the system. There
are five types of local, discrete-time bifurcations of codimension one: (a)
period-doubling (b) Neimark-Sacker (c) fold (d) transcritical (e) pitchfork
Bifurcation diagrams are shown in the panels below [7]. The $x$-axis
represents the bifurcation parameter and the $y$-axis a state variable. Solid
lines show stable states or limit cycles, and dashed lines show unstable
states. Arrows indicate direction of motion at different points in state
space. The white points on the unit circle show the value of the dominant
eigenvalue(s) at the bifurcation. The Neimark-Sacker bifurcation has a pair of
complex conjugate eigenvalues, the angle ($\theta$) being the frequency of
oscillations at the bifurcation.

## Results

Performance of classifiers on withheld test data
As in previous work [25], we train two different types of classifiers and use
their ensemble average to make predictions. Classifier 1 is trained to
recognise bifurcation trajectories based on middle portions of the time
series, whereas Classifier 2 is trained on end portions (see Methods). In this
way, Classifier 1 provides an earlier signal of a bifurcation and Classifier 2
provides a more specific signal, as more information is revealed closer to the
bifurcation. To quantify performance of the classifiers, we use the F1 score,
which is a combined measure of sensitivity (how many of the true positives
were predicted correctly) and specificity (how many of the positive
predictions were actually true positives). On the withheld test data,
Classifier 1 and 2 achieved an F1 score of $0.66$ and $0.85$, respectively. On
the simpler, binary classification problem of predicting whether or not there
will be any bifurcation, the classifiers achieved an F1 score of $0.79$ and
$0.97$, respectively. Classifier 2 has a higher performance as it has the
easier task of classifying data closer to the bifurcation where fluctuations
are more pronounced. Performance on individual bifurcation classes is shown by
confusion matrices (Sup. Fig. 1). The period-doubling, Neimark-Sacker and fold
bifurcations are correctly classified with high sensitivity and specificity.
On the other hand, the transcritical and pitchfork bifurcations are often
mistaken for one another, likely due to having very similar normal forms
(identical linear terms). Despite this, Classifier 2 can distinguish them at
better than random, suggesting it is capable of recognising the different
higher order terms in the data. From here onward, we report results using the
ensemble prediction of the two classifiers, referred to collectively as the
deep learning classifier.

Performance of EWS on theoretical models
We monitor variance, lag-1 autocorrelation and the deep learning classifier as
progressively more of the time series is revealed. Variance and lag-1
autocorrelation are considered to provide an EWS if they display a strong
trend, which we quantify using the Kendall tau statistic. For each of the five
theoretical models (Fig. 2(a-e)), we observe an increasing trend in variance
(Fig. 2f-j), and an increasing or decreasing trend in lag-1 autocorrelation
(Fig. 2k-o). The direction of the trend in lag-1 autocorrelation prior to a
bifurcation depends on the frequency of oscillations ($\theta$) at the
bifurcation—equivalently the angle of the dominant eigenvalue from the
positive x-axis (Box 1). For $\theta\in[0,\pi/2)$ lag-1 autocorrelation
increases, whereas for $\theta\in(\pi/2,\pi]$ it decreases—insights that can
be obtained from analytical expressions of the autocorrelation function [14,
10]. The period-doubling bifurcation is characterised by $\theta=\pi$, and the
Neimark-Sacker bifurcation shown here has $\theta\approx\pi/4$. The trends in
variance and lag-1 autocorrelation therefore behave as expected and can be
used as an EWS.

Figure 2: Trends in indicators prior to five different bifurcations in the
theoretical models. (a-e) Trajectory (gray) and smoothing (black) of a
simulation going through a period-doubling, Neimark-Sacker, fold,
transcritical and pitchfork bifurcation, respectively. (f-j) Variance of
residual dynamics after smoothing, computed over a rolling window (arrow) of
size 0.5 times the length of the pre-transition data. (k-o) Lag-1
autocorrelation. (p-t) Probabilities assigned by the deep learning (DL)
classifier when given all preceding data. Orange line shows probability
assigned to the true bifurcation. Gray lines show probabilities assigned to
the other bifurcations. Blue line shows the sum of the five bifurcation
probabilities.

The deep learning classifier assigns a probability to each of the six possible
outcomes (null, period-doubling, fold, Neimark-Sacker, transcritical and
pitchfork). It is considered to provide an EWS when there is a heightening in
the sum of the bifurcation probabilities (blue line, Fig. 2p-t). The type of
bifurcation predicted is then taken as the highest individual bifurcation
probability. For each simulation, the classifier is becomes more confident of
an approaching bifurcation as time goes on, and its assigned bifurcation
probability for the true bifurcation increases. The period-doubling, Neimark-
Sacker and fold bifurcations are identified with high confidence well before
the transition. The transcritical and pitchfork bifurcations are assigned
roughly equal probability on their respective time series, suggesting they are
difficult to tell apart—an observation consistent with the classifier’s
performance on its within-sample test data.

To obtain a measure of performance for the EWS, we need to test their
predictions on both ‘forced’ time series (where a bifurcation is approached)
and ‘null’ time series (where no bifurcation is approached). For the
theoretical models, we generate null time series by keeping the bifurcation
parameter fixed. We also test the robustness of the EWS to the rate of forcing
and the noise amplitude of the simulations—two factors that have been shown to
influence the performance of variance and lag-1 autocorrelation as an EWS [33,
34]. To this end, we simulate 100 forced and null time series at five
different values of noise intensity and five different values of rate of
forcing, resulting in a total of 5,000 time series for each theoretical model.
Sample trajectories illustrating the different noise amplitudes and rates of
forcing are shown in Sup. Fig. 3. We compute the probabilities assigned by the
classifier and the Kendall tau values for variance and lag-1 autocorrelation
at $80\%$ of the way through the pretransition time series, and use these
values as discrimination thresholds to construct ROC curves (Fig. 4a-e). Using
the AUC score (area under the ROC curve) as a measure of performance, we find
that the classifier outperforms variance and lag-1 autocorrelation for each
theoretical model. When evaluated for each combination of noise amplitude and
rate of forcing separately, the classifier has the highest AUC score in
$100\%$ of cases for the Neimark-Sacker, fold, and pitchfork models, $84\%$ of
cases for the period-doubling model, and $80\%$ of cases for the transcritical
model (Sup. Fig. 4). Similar to variance and lag-1 autocorrelation, the
performance of the classifier is lower at higher rates of forcing. Noise
amplitude affects performance differently depending on the model. In terms of
predicting the correct bifurcation, the classifier typically performs better
at slower rates of forcing (Sup. Fig. 5) and was able to classify the period-
doubling and Neimark-Sacker bifurcations to high accuracy at all noise
amplitudes and rates of forcing considered.

Figure 3: Trends in indicators prior to a period-doubling bifurcation in five
chick heart aggregates treated with a potassium channel blocker. (a-e) Inter-
beat interval (gray) and smoothing. (f-j) Variance of residual dynamics after
smoothing, computed over a rolling window (arrow) of size 0.5 times the length
of the pre-transition data. (k-o) Lag-1 autocorrelation. (p-t) Probabilities
assigned by the deep learning (DL) classifier to the period-doubling
bifurcation (orange) and the other bifurcations (gray). Blue line shows the
sum of the five bifurcation probabilities.

Performance of EWS on chick heart data
In the chick heart data, we mostly observe an increasing trend in variance and
a decreasing trend in lag-1 autocorrelation prior to the period-doubling
bifurcation, as previously reported [21]. This is consistent with analytical
approximations for variance and lag-1 autocorrelation prior to a period-
doubling bifurcation in a noisy dynamical system [14, 10]. The 46 records and
their EWS are shown in Sup. Fig. 6-9, and a sample of five period-doubling
records are shown in Fig. 3. The classifier correctly predicts a period-
doubling bifurcation in 16 of the 23 period-doubling records. In other cases,
it incorrectly predicts a Neimark-Sacker bifurcation (e.g. Fig. 3e). This
seems to be linked to an early increase in lag-1 autocorrelation, perhaps
caused by a non-monotonic approach to the period-doubling bifurcation. For
predictions made at 60-100% of the way through the chick heart data, the
classifier obtains the highest AUC score (Fig. 4f), a slight improvement on
variance, with the advantage of also providing the bifurcation type.

Figure 4: ROC curves for predictions of an upcoming transition in model and
experimental data. ROC curves compare the performance of the deep learning
classifier (DL, blue), variance (Var, red) and lag-1 autocorrelation (AC,
green). For models (a-e), performance is assessed on 5,000 simulations with
different noise amplitudes and rates of forcing. For experimental data (f),
performance is assessed on 46 experimental runs. The area under the curve
(AUC), abbreviated to A, is a measure of performance. Insets show the
probabilities assigned by the classifier to each type of bifurcation (orange
being the true bifurcation) among the trajectories approaching a transition.
(a) Fox model going through a period-doubling bifurcation [35]. (b) Westerhoff
model going through a Neimark-Sacker bifurcation [30]. (c) Ricker model going
through a fold bifurcation [36]. (d) Lotka-Volterra model going through a
transcritical bifurcation [37]. (e) Lorenz model going through a pitchfork
bifurcation [38]. (f) Chick heart aggregates going through a period-doubling
bifurcation [31]. Predictions are made 80% of the way through the
pretransition data for the models, and 60-100% of the way for the experimental
data. PD: period-doubling. NS: Neimark-Sacker. TC: transcritical. PF:
pitchfork.

## Discussion

Many systems that evolve on a discrete timeline can undergo a sudden change in
dynamics via a discrete-time bifurcation. We have found that a deep learning
classifier is an effective tool for predicting discrete-time bifurcations in
systems with a range of noise levels and rates of approach to the bifurcation.
The classifier provides better sensitivity and specificity than variance and
lag-1 autocorrelation—two commonly used EWS for bifurcations. Moreover, the
classifier provides early indication of the type of bifurcation—an important
piece of information given the qualitatively different dynamics associated
with each bifurcation. A reliable early warning signal that specifies the type
of bifurcation will help us prevent harmful bifurcations (e.g. dangerous heart
rhythms [3]) and promote favourable transitions (e.g. ecosystem recovery
[39]).

It may be possible to design a deep learning classifier that achieves a higher
performance on our test data. First, there are many neural network
architectures that could be investigated. For example, transformers, which are
the current state-of-the-art for language models like GPT [40], may also be
useful for time series classification [41]. Second, the hyperparameters of the
classifier could be systematically tuned to optimise performance. Third, there
may be benefit to reframing bifurcation prediction as a hierarchical
classification problem [42]. One classifier could address the binary problem
of flagging an approaching bifurcation, and a second classifier could address
the multi-class problem of classifying the type of bifurcation _given_ that a
bifurcation is approaching. (The same way one might want to distinguish images
of dogs and cats before attempting to classify dog breeds.)

For a classifier to be effective, it must be trained on sufficiently diverse
training data. As such, the method by which training data is obtained needs
careful consideration. Our previous work on continuous-time bifurcations
obtained training data from randomly generated dynamical systems with
polynomial terms [25] and labelled the data using the bifurcation continuation
software AUTO. This approach is appealing as it imposes relatively few
restrictions on the models that are generated, and may include features
associated with higher-order terms. Here, we opted for a more restricted
approach that uses normal form models to generate the training data. This
method has the advantage of being faster computationally, since the
bifurcation of the model is known a priori. This study shows that even with
this more restricted training data, a classifier can generalise to detecting
bifurcations in more complex model and empirical systems.

We trained a classifier to provide EWS for a subset of bifurcations, namely
local, codimension-one discrete-time bifurcations. While these bifurcations
are present in many systems of interest, the real world presents many other
classes of bifurcation in both continuous and discrete-time, including global
bifurcations (e.g. homoclinic and heteroclinic), codimension-two bifurcations
(e.g. cusp and Bogdanov-Takens), and bifurcations of attractors. For systems
on attractors that explore a large portion of their phase space, empirical
dynamical modelling [43], reservoir computing [44, 45] and deep neural
networks [46] can be used to make forecasts that may help predict critical
transitions. In cases where spatial information is available, concepts from
statistical physics may be useful [47], particularly in combination with deep
learning [27].

One limitation of the present classifier is that it is only trained to predict
discrete-time bifurcations. Therefore, one needs to know ahead of time whether
continuous or discrete time is a better description for the system. In the
case of the chick heart cells, we had prior knowledge that they are well
described by a discrete-time dynamical system [48], and therefore appropriate
for the classifier. An interesting avenue for future research is to build a
classifier that works for both continuous and discrete-time bifurcations. This
may be achieved by generating a training library from models with a range of
discretised timesteps, from very large steps that generate discrete-time
bifurcations, down to the limit of a discrete timestep of zero, where
continuous-time bifurcations occur. With a large enough training set, one
would not need to assume ahead of time whether continuous or discrete time is
a better description for the system.

Our results demonstrate that combining dynamical system and deep learning
methodologies can provide EWS for critical transitions that are both more
reliable and more descriptive than non-hybrid approaches. This study has set a
baseline for prediction performance across a variety of popular discrete-time
models and an experimental data set. In providing a code capsule that
reproduces this study, we hope to facilitate the development, testing and
comparison of related methods. Building a universal predictor for critical
transitions is not a job for a single research team [46], and will benefit
from a variety of approaches and open source code. Depending on context,
critical transitions can be devastating or highly desirable. Improved EWS
would allow us to better prevent or promote such transitions.

## Methods

Generation of training data for the deep learning classifier
Training data consists of simulation data from a library of 50,000 models. The
models are generated at random from five different model frameworks, each
possessing one of the bifurcations studied (period-doubling, Neimark-Sacker,
fold, transcritical, pitchfork). The models are composed of the normal form of
the bifurcation [7], higher-order polynomial terms up to degree 10 with
coefficients drawn from a normal distribution, and additive Gaussian white
noise ($\epsilon_{t}$) with amplitude ($\sigma$) drawn from a uniform
distribution. In each case, the bifurcation occurs at $\mu=0$.

The model for the period-doubling bifurcation is

$x_{t+1}=-(1+\mu)x_{t}\pm
x_{t}^{3}+\sum_{i=4}^{10}\alpha_{i}x_{t}^{i}+\sigma\epsilon_{t},$ (2)

where $\alpha_{i}\sim\mathcal{N}(0,1)$. The positive (negative) cubic term
yields a supercritical (subcritical) bifurcation, and is chosen at random. The
model for the Neimark-Sacker bifurcation is

$\begin{pmatrix}x_{t+1}\\\
y_{t+1}\end{pmatrix}=(1+\mu)R(\theta)\begin{pmatrix}x_{t}\\\
y_{t}\end{pmatrix}\pm(x_{t}^{2}+y_{t}^{2})R(\theta)\begin{pmatrix}x_{t}\\\
y_{t}\end{pmatrix}+\sum_{i=4}^{10}\sum_{j=0}^{i}\begin{pmatrix}\alpha_{ij}\\\
\beta_{ij}\end{pmatrix}x_{t}^{i-j}y_{t}^{j}+\sigma\begin{pmatrix}\epsilon_{t}^{(1)}\\\
\epsilon_{t}^{(2)}\end{pmatrix}$ (3)

where $\alpha_{ij}$, $\beta_{ij}\sim\mathcal{N}(0,1)$, $R(\theta)$ is the
rotation matrix

$R(\theta)=\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix},$ (4)

and $\theta\sim\mathcal{U}[0,\pi]$ is the angular frequency of oscillations at
the bifurcation. The positive (negative) cubic term yields a subcritical
(supercritical) bifurcation, and is chosen at random. The model for the fold
bifurcation is

$x_{t+1}=-\mu+x_{t}-x_{t}^{2}+\sum_{i=3}^{10}\alpha_{i}(x_{t}-\sqrt{-\mu})^{i}+\sigma\epsilon_{t},$
(5)

where $\alpha_{i}\sim\mathcal{N}(0,1)$. The model for the transcritical
bifurcation is

$x_{t+1}=(1+\mu)x_{t}-x_{t}^{2}+\sum_{i=3}^{10}\alpha_{i}x_{t}^{i}+\sigma\epsilon_{t},$
(6)

where $\alpha_{i}\sim\mathcal{N}(0,1)$. Finally, the model for the pitchfork
bifurcation is

$x_{t+1}=(1+\mu)x_{t}\pm
x_{t}^{3}+\sum_{i=4}^{10}\alpha_{i}x_{t}^{i}+\sigma\epsilon_{t},$ (7)

where $\alpha_{i}\sim\mathcal{N}(0,1)$. The positive (negative) cubic term
yields a subcritical (supercritical) bifurcation, and is chosen at random.

The library is composed of 10,000 models from each framework. For each model,
we run a ‘forced’ simulation where the bifurcation parameter $\mu$ is
increased linearly across the interval $[\mu_{0},0]$, and a ‘null’ simulation
where $\mu$ is fixed at $\mu_{0}$. The initial value for the bifurcation
parameter $\mu_{0}$ is drawn from a uniform distribution across all values
that correspond to $|\lambda|<0.8$, where $\lambda$ is the eigenvalue of the
Jacobian matrix in the model. This ensures that the training data contains
simulations that start close to and far from a bifurcation. For the period-
doubling, Neimark-Sacker, transcritical and pitchfork models, this means
drawing $\mu_{0}$ from $\mathcal{U}[-1.8,-0.2]$. For the fold bifurcation,
this means drawing $\mu_{0}$ from $\mathcal{U}[-0.9,-0.1]$. After a burn-in
period of 100 iterations, we simulate each model for 600 iterations and keep
the last 500 data points, or the 500 data points immediately preceding a
transition if one occurs. We define a transition as a time when the deviation
from equilibrium exceeds ten times the noise amplitude $\sigma$. We simulate
one forced and one null simulation from each model, resulting in $50,000$
forced and $50,000$ null trajectories. To balance the number of entries for
each class, we take $10,000$ null simulations at random, resulting in a total
of $60,000$ entries in the training data set. Example trajectories for each
class are shown in Sup. Fig. 2.

Architecture and training of the deep learning classifier
We use a neural network with a CNN-LSTM architecture and hyperparameters as in
[25]. This consists of a single convolutional layer with max pooling followed
by two LSTM layers with dropout followed by a dense layer that maps to a
vector of probabilities over the six possible classes. For training, we use
Adam optimisation with a learning rate of $0.0005$, a batch size of $1024$,
and sparse categorical cross entropy as the loss function. We use a
training/validation/test split of $0.95/0.025/0.025$. We found 200 epochs was
sufficient to obtain optimal accuracy on the validation set.

To expose the classifier to time series of different lengths, we censor each
time series in the training data. We train two classifiers independently using
different censoring techniques. Classifier 1 is trained on time series
censored at the beginning and the end, forcing it to learn from data in the
middle of the time series. Classifier 2 is trained on time series only
censored at the beginning, allowing it to learn from data right up to the
bifurcation. The length for each censored time series $L$ is drawn from
$\mathcal{U}[50,500]$. Then, for Classifier 1, the start time of the censored
time series $t_{0}\sim U[0,500-L]$ and for Classifier 2, $t_{0}=500-L$. The
censored time series are then normalised by their mean absolute value and
prepended with zeros to make them $500$ points in length. We report results
using the average prediction of the two classifiers.

Theoretical models
To test the deep learning classifier on out-of-sample data, we simulate a
variety of nonlinear, discrete-time models, each containing one of the studied
bifurcations. To account for stochasticity, we include additive Gaussian white
noise. We run forced simulations, where the bifurcation parameter is increased
linearly up to the bifurcation point, and null simulations, where the
bifurcation parameter remains fixed. To create a diverse set of test data, we
vary the noise amplitude ($\sigma$) and rate of forcing (rate of change of the
bifurcation parameter). We run 100 forced and null simulations at five
different noise amplitudes and five different rates of forcing, resulting in
5000 simulations of each model. Values for the noise amplitude are on a
logarithmic scale and values for the rate of forcing result in time series of
length 100, 200, 300, 400, and 500. Sample simulations for each model at
different noise amplitude and rate of forcing are shown in Sup. Fig. 3. The
transition time for each forced simulation is taken as the moment when the
bifurcation parameter crosses the bifurcation, or the moment when the state
variable crosses a threshold, if specified.

(1) Fox model. To test detection of a period-doubling bifurcation, we use a
model of cardiac alternans [35] with additive Gaussian white noise. This is
given by

$\displaystyle D_{n+1}$ $\displaystyle=(1-\alpha
M_{n+1})\left(A+\frac{B}{1+e^{-(I_{n}-C)/D}}\right)+\sigma\epsilon_{n},$ (8)
$\displaystyle M_{n+1}$
$\displaystyle=e^{-I_{n}/\tau}[1+(M_{n}-1)e^{-D_{n}/\tau}],$ (9)
$\displaystyle I_{n}$ $\displaystyle=T-D_{n},$ (10)

where $D_{n}$ is the action potential duration of the $n$th beat, $M_{n}$ is a
memory variable, $I_{n}$ is the rest duration following the action potential,
$T$ is the stimulation period, $\tau$ is the time constant of accumulation and
dissipation of memory, $\alpha$ is the influence of memory on the action
potential duration, and $A$, $B$, $C$ and $D$ are parameters governing the
shape of the restitution curve. Following [35], we take $A=88$, $B=122$,
$C=40$, $D=28$, $\tau=180$, $\alpha=0.2$, which give dynamics in good
agreement with a complex ionic model. This yields a period-doubling
bifurcation at approximately $T=200$. Forced simulations are run with $T$
decreasing linearly on the interval $[300,150]$ and null simulations are run
with $T=300$. Values for noise amplitude are
$0.1\times\\{2^{0},2^{-1},2^{-2},2^{-3},2^{-4}\\}$.

(2) Westerhoff model. To test detection of a Neimark-Sacker bifurcation, we
use a simple model of business cycles based on consumer sentiment [30] with
additive Gaussian noise. This is given by

$Y_{t}=a+(b-d)Y_{t-1}+dY_{t-2}+\frac{cY_{t-1}}{1+\text{Exp}\left[-(Y_{t-1}-Y_{t-2})\right]}+\sigma\epsilon_{t},$
(11)

where $Y_{t}$ is the national income at time step $t$, $a$ is the level of
autonomous expenditures of agents, $b$ and $c$ govern a curve that determines
the fraction of income consumed by the agents, and $d$ is the policy-maker’s
control parameter to offset income trends. Following [30], we take $b=0.45$,
$c=0.1$, and $d=0$, which yields a Neimark-Sacker bifurcation at $a=20$
corresponding to the onset of business cycles. Forced simulations are run with
$a$ increasing linearly on the interval $[10,22.5]$ and null simulations are
run with $a=10$. Values for noise amplitude are
$0.1\times\\{2^{0},2^{-1},2^{-2},2^{-3},2^{-4}\\}$.

(3) Ricker model. To test detection of a fold bifurcation, we use the Ricker
model [36] with a sigmoidal harvesting term and additive Gaussian noise. This
is given by

$x_{t+1}=x_{t}e^{r(1-x_{t}/k)}-F\frac{x_{t}^{2}}{x_{t}^{2}+h^{2}}+\sigma\epsilon_{t},$
(12)

where $x_{t}$ is the population size at time step $t$, $r$ is the intrinsic
growth rate, $k$ is the carrying capacity, $F$ is the harvesting rate, and $h$
governs the steepness of the sigmoidal harvesting term. We take $r=0.75$,
$k=10$, $h=0.75$, which yields a fold bifurcation at $F=2.36$. Forced
simulations are run with $F$ increasing linearly on the interval $[0,3.54]$
and null simulations are run with $F=0$. We define a transitions as the time
when $x_{t}$ drops below 0.45. Values for noise amplitude are
$0.2\times\\{2^{0},2^{-1},2^{-2},2^{-3},2^{-4}\\}$.

(4) Lotka-Volterra model. To test detection of a transcritical bifurcation, we
use the discrete-time analogue of the Lotka-Volterra model, first studied by
Maynard Smith [37]. This system is especially relevant to arthropod predator-
prey and host-parasitoid interactions. The rescaled equations [49] are

$\displaystyle x_{t+1}$
$\displaystyle=(r+1)x_{t}-rx_{t}^{2}-cx_{t}y_{t}+\sigma\epsilon_{t}^{(1)},$
(13) $\displaystyle y_{t+1}$
$\displaystyle=cx_{t}y_{t}+\sigma\epsilon_{t}^{(2)},$ (14)

$r$ relates to the growth rate of the prey ($x_{t}$), and $c$ relates to the
foraging efficiency of the predator ($y_{t}$). We take $r=0.5$, which yields a
transcritical bifurcation at $c=1$, the critical foraging efficiency beyond
which the predator population can sustain themselves. Forced simulations are
run with $c$ increasing linearly on the interval $[0.5,1.25]$ and null
simulations are run with $c=0.5$. We look for early warning signals in the
prey population. Values for noise amplitude are
$0.01\times\\{2^{0},2^{-1},2^{-2},2^{-3},2^{-4}\\}$.

(5) Lorenz model. To test detection of a pitchfork bifurcation, we use the
reduced discrete Lorenz system, which was first introduced as a demonstration
of computational chaos [38]. This is given by

$\displaystyle x_{t+1}=(1+ah)x_{t}-hx_{t}y_{t}+\sigma\epsilon_{t}^{(1)}$ (15)
$\displaystyle y_{t+1}=(1-h)y_{t}+hx_{t}^{2}+\sigma\epsilon_{t}^{(2)}$ (16)

where state variables and parameters are derived from the full Lorenz
equations [38]. We take $h=0.5$, which yields a pitchfork bifurcation at
$a=0$. Forced simulations are run with $a$ increasing linearly over the
interval $[-1,0.25]$ and null simulations are run with $a=-1$. We look for
early warning signals in $x_{t}$. Values for noise amplitude are
$0.01\times\\{2^{0},2^{-1},2^{-2},2^{-3},2^{-4}\\}$.

Experimental data
We use data from spontaneously beating aggregates of chick heart cells, which
can undergo a transition to an alternating rhythm with the addition of the
potassium-channel blocker E-4031 [31]. To study the dynamics of these
aggregates, we focus on the interbeat intervals—the time between consecutive
beats. The transition to an alternating rhythm may be associated with a
period-doubling bifurcation. Trends in variance and autocorrelation of the
data have been shown to agree with theoretical expectations for a period-
doubling bifurcation [21]. We define the onset of the period-doubling
bifurcation as the first time when the slope of a linear regression of the
return map composed of a sliding window of interbeat intervals is below -0.95
for the next 10 beats. We have 23 time series that undergo a period-doubling
bifurcation (Sup. Fig. 5, 6), and 18 time series that do not, from which we
extract 23 segments at random with a random length between 100 and 500, to
serve as null time series (Sup. Fig. 7, 8). For details on the chick heart
experiments, see [31, 21].

Computing and assessing the performance of EWS
EWS are computed using the Python package _ewstools_ [50]. This involves first
detrending the (pretransition) time series. For the model simulations, we use
a Lowess filter with a span of 0.25 the length of the data. For the heart cell
data, we use a Gaussian filter with a bandwidth of 20 beats. Variance and
lag-1 autocorrelation are then computed over a rolling window of 0.5, which
had higher performance than a rolling window of 0.25. The deep learning
predictions at a given point in the time series are obtained by taking the
preceding data from the time series, normalising it, prepending it with zeroes
to make it 500 points in length, and feeding it into the classifier.

To compare performance of variance, lag-1 autocorrelation and the deep
learning classifier, we use the AUC (area under curve) score of the ROC curve.
The ROC curve plots the true positive rate vs. the false positive rate as a
discrimination threshold is varied. We use the Kendall $\tau$ value as the
discrimination threshold for variance and lag-1 autocorrelation, and the sum
of the bifurcation probabilities for the discrimination threshold of the deep
learning classifier.

## Data availability

Data from the chick heart aggregates, as well as synthetic data used to train
and test the deep learning classifier are available at the GitHub repository
https://github.com/ThomasMBury/dl_discrete_bifurcation.

## Code availability

Code and instructions to reproduce the analysis of this article are available
at the GitHub repository
https://github.com/ThomasMBury/dl_discrete_bifurcation. A reproducible run can
be performed on Code Ocean at https://codeocean.com/capsule/2209652/tree/v1
where the code is accompanied by the necessary software environment.

## Acknowledgements

This research was supported by a Fonds de Recherche du Québec – Nature et
technologies (FRQNT) postdoctoral fellowship to T.M.B, a Canadian Institutes
of Health Research (CIHR) grant (#PJT-169008) to A.S., and compute services
from the Digital Research Alliance of Canada (www.alliancecan.ca).

## References

  * [1] Marten Scheffer  “Critical transitions in nature and society”  Princeton University Press, 2020
  * [2] Simon A Levin  “Ecosystems and the biosphere as complex adaptive systems”  In _Ecosystems_ 1 Springer, 1998, pp. 431–436
  * [3] Leon Glass and Michael C Mackey  “From Clocks to Chaos”  Princeton University Press, 2020
  * [4] Didier Sornette  “Why Stock Markets Crash”  Princeton University Press, 2017
  * [5] Lee-Ann Barlow, Jacob Cecile, Chris T Bauch and Madhur Anand  “Modelling interactions between forest pest invasions and human decisions regarding firewood transport restrictions”  In _PloS one_ 9.4, 2014, pp. e90511
  * [6] Kirsten A Henderson, Chris T Bauch and Madhur Anand  “Alternative stable states and the sustainability of forests, grasslands, and agriculture”  In _Proceedings of the National Academy of Sciences_ 113.51, 2016, pp. 14552–14559
  * [7] Yuri Kuznetsov  “Elements of Applied Bifurcation Theory”  Springer, 1998
  * [8] Steven H Strogatz  “Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering”  CRC press, 2018
  * [9] C Wissel  “A universal law of the characteristic return time near thresholds”  In _Oecologia_ 65.1, 1984, pp. 101–107
  * [10] Kurt Wiesenfeld  “Noisy precursors of nonlinear instabilities”  In _Journal of Statistical Physics_ 38.5, 1985, pp. 1071–1097
  * [11] Marten Scheffer et al.  “Early-warning signals for critical transitions”  In _Nature_ 461.7260, 2009, pp. 53–59
  * [12] Christian Kuehn  “A mathematical framework for critical transitions: normal forms, variance and applications”  In _Journal of Nonlinear Science_ 23.3, 2013, pp. 457–510
  * [13] Suzanne M O’Regan and Danielle L Burton  “How stochasticity influences leading indicators of critical transitions”  In _Bulletin of Mathematical Biology_ 80.6, 2018, pp. 1630–1654
  * [14] Thomas M Bury, Chris T Bauch and Madhur Anand  “Detecting and distinguishing tipping points using spectral early warning signals”  In _Journal of the Royal Society Interface_ 17.170, 2020, pp. 20200482
  * [15] Vasilis Dakos et al.  “Slowing down as an early warning signal for abrupt climate change”  In _Proceedings of the National Academy of Sciences_ 105.38, 2008, pp. 14308–14312
  * [16] Niklas Boers  “Early-warning signals for Dansgaard-Oeschger events in a high-resolution ice core record”  In _Nature Communications_ 9.1, 2018, pp. 2556
  * [17] Niklas Boers  “Observation-based early-warning signals for a collapse of the Atlantic Meridional Overturning Circulation”  In _Nature Climate Change_ 11.8, 2021, pp. 680–688
  * [18] Rick Hennekam et al.  “Early-warning signals for marine anoxic events”  In _Geophysical Research Letters_ 47.20, 2020, pp. e2020GL089183
  * [19] Michael L Pace et al.  “Reversal of a cyanobacterial bloom in response to early warnings”  In _Proceedings of the National Academy of Sciences_ , 2016, pp. 201612424
  * [20] Rong Wang et al.  “Flickering gives early warning signals of a critical transition to a eutrophic lake state”  In _Nature_ 492.7429 Nature Publishing Group UK London, 2012, pp. 419–422
  * [21] Thomas Quail, Alvin Shrier and Leon Glass  “Predicting the onset of period-doubling bifurcations in noisy cardiac systems”  In _Proceedings of the National Academy of Sciences_ 112.30, 2015, pp. 9358–9363
  * [22] Sonia Kéfi et al.  “Early warning signals also precede non-catastrophic transitions”  In _Oikos_ 122.5, 2013, pp. 641–648
  * [23] Alan Hastings and Derin B Wysham  “Regime shifts in ecological systems can occur with no warning”  In _Ecology letters_ 13.4, 2010, pp. 464–472
  * [24] Peter D Ditlevsen and Sigfus J Johnsen  “Tipping points: Early warning and wishful thinking”  In _Geophysical Research Letters_ 37.19, 2010
  * [25] Thomas M Bury et al.  “Deep learning for early warning signals of tipping points”  In _Proceedings of the National Academy of Sciences_ 118.39, 2021, pp. e2106140118
  * [26] Smita Deb et al.  “Machine learning methods trained on simple models can predict critical transitions in complex natural systems”  In _Royal Society Open Science_ 9.2, 2022, pp. 211475
  * [27] Daniel Dylewsky et al.  “Universal Early Warning Signals of Phase Transitions in Climate Systems”  In _arXiv_ , 2022
  * [28] Robert M May  “Biological populations with nonoverlapping generations: stable points, stable cycles, and chaos”  In _Science_ 186.4164, 1974, pp. 645–647
  * [29] Linda JS Allen  “Some discrete-time SI, SIR, and SIS epidemic models”  In _Mathematical Biosciences_ 124.1, 1994, pp. 83–105
  * [30] Frank H Westerhoff  “Consumer sentiment and business cycles: a Neimark–Sacker bifurcation scenario”  In _Applied Economics Letters_ 15.15, 2008, pp. 1201–1205
  * [31] Min-Young Kim et al.  “Stochastic and spatial influences on drug-induced bifurcations in cardiac tissue culture”  In _Physical Review Letters_ 103.5, 2009, pp. 058101
  * [32] Richard L Verrier et al.  “Microvolt T-wave alternans: physiological basis, methods of measurement, and clinical utility—consensus guideline by International Society for Holter and Noninvasive Electrocardiology”  In _Journal of the American College of Cardiology_ 58.13, 2011, pp. 1309–1324
  * [33] Christopher F Clements and Arpat Ozgul  “Rate of forcing and the forecastability of critical transitions”  In _Ecology and Evolution_ 6.21, 2016, pp. 7787–7793
  * [34] Induja Pavithran and RI Sujith  “Effect of rate of change of parameter on early warning signals for critical transitions”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 31.1, 2021, pp. 013116
  * [35] Jeffrey J Fox, Eberhard Bodenschatz and Robert F Gilmour Jr  “Period-doubling instability and memory in cardiac tissue”  In _Physical Review Letters_ 89.13, 2002, pp. 138101
  * [36] William Edwin Ricker  “Stock and recruitment”  In _Journal of the Fisheries Board of Canada_ 11.5, 1954, pp. 559–623
  * [37] J Maynard Smith  “Mathematical Ideas in Biology”  CUP Archive, 1968
  * [38] Edward N Lorenz  “Computational chaos-a prelude to computational instability”  In _Physica D: Nonlinear Phenomena_ 35.3, 1989, pp. 299–317
  * [39] Christopher F Clements, Michael A McCarthy and Julia L Blanchard  “Early warning signals of recovery in complex systems”  In _Nature Communications_ 10.1, 2019, pp. 1681
  * [40] Ashish Vaswani et al.  “Attention is all you need”  In _Advances in Neural Information Processing Systems_ 30, 2017
  * [41] Qingsong Wen et al.  “Transformers in time series: A survey”  In _arXiv_ , 2022
  * [42] Carlos N Silla and Alex A Freitas  “A survey of hierarchical classification across different application domains”  In _Data Mining and Knowledge Discovery_ 22, 2011, pp. 31–72
  * [43] Hao Ye et al.  “Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling”  In _Proceedings of the National Academy of Sciences_ 112.13, 2015, pp. E1569–E1576
  * [44] Dhruvit Patel and Edward Ott  “Using machine learning to anticipate tipping points and extrapolate to post-tipping dynamics of non-stationary dynamical systems”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 33.2, 2023, pp. 023143
  * [45] Ling-Wei Kong, Hua-Wei Fan, Celso Grebogi and Ying-Cheng Lai  “Machine learning prediction of critical transition and system collapse”  In _Physical Review Research_ 3.1, 2021, pp. 013090
  * [46] Marcus Lapeyrolerie and Carl Boettiger  “Teaching machines to anticipate catastrophes”  In _Proceedings of the National Academy of Sciences_ 118.40, 2021, pp. e2115605118
  * [47] George I. Hagstrom and Simon A. Levin  “Phase Transitions and the Theory of Early Warning Indicators for Critical Transitions”  In _arXiv_ , 2021
  * [48] Thomas Quail et al.  “Chaotic dynamics in cardiac aggregates induced by potassium channel block”  In _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 22.3, 2012, pp. 033140
  * [49] Michael G Neubert and Mark Kot  “The subcritical collapse of predator populations in discrete-time predator-prey models”  In _Mathematical Biosciences_ 110.1, 1992, pp. 45–66
  * [50] Thomas M Bury  “ewstools: A Python package for early warning signals of bifurcations in time series data”  In _Journal of Open Source Software_ 8.82, 2023, pp. 5038


[pages=1-8,color=gray!20,scale=3,angle=-45,xpos=0,ypos=0]IEEE ICCWS 2022

# Sequential Embedding-based Attentive (SEA) classifier for malware
classification

1st Muhammad Ahmed  Department of Computer Science
National University of Computer and Emerging Sciences
Karachi, Pakistan
<EMAIL_ADDRESS>2nd Anam Qureshi  Department of Computer Science
National University of Computer and Emerging Sciences
Karachi, Pakistan
<EMAIL_ADDRESS>3rd Jawwad Ahmed Shamsi  Department of Computer
Science
National University of Computer and Emerging Sciences
Karachi, Pakistan
<EMAIL_ADDRESS>4th Murk Marvi  Department of Computer Science
NED University of Engineering and Technology
Karachi, Pakistan
<EMAIL_ADDRESS>
###### Abstract

The tremendous growth in smart devices has uplifted several security threats.
One of the most prominent threats is malicious software also known as malware.
Malware has the capability of corrupting a device and collapsing an entire
network. Therefore, it’s early detection and mitigation are extremely
important to avoid catastrophic effects. In this work, we came up with a
solution for malware detection using state-of-the-art natural language
processing (NLP) techniques. Our main focus is to provide a lightweight yet
effective classifier for malware detection which can be used for heterogeneous
devices, be it a resource constraint device or a resourceful machine. Our
proposed model is tested on the benchmark data set with an accuracy and log
loss score of 99.13% and 0.04 respectively.

###### Index Terms:

Security, Smart Devices, Malware Detection, Natural Language Processing

††publicationid: pubid: 978-1-6654-8939-3/22/$31.00 ©2022 IEEE

##  I Introduction

Malware is one of the biggest threats to the smart devices. According to a
recent Open Web Application Security Project (OWASP), malicious software is
among the top 5 threats to the smart devices [29]. A smart device holds
different types of information i.e., health, traffic, and agriculture data.
Some of the information is sensitive such as health data. Therefore, it is
important to protect the smart devices from malicious activities. Moreover, a
smart device has limited resources such as memory and processing. It is one of
the reasons these devices are an easy target for the attackers [30] as they
cannot run bulky antivirus programs. Therefore, there is a need of a light
weight malware detection technique which has capability to be integrated
within the embedded devices. In addition to this, the time taken for malware
detection should be minimal in order to protect the devices before it gets
exploited.

Figure 1: Abstract view of Sequential Embedding-based Attentive Model

There are three prominent techniques that are used to detect and classify a
malware [18]. Following are the details:

###  I-A Static Analysis (Signature Based)

In static analysis, program’s own structure is compared with the structure of
the known malware without considering the execution behavior and other
factors. This technique can be helpful where a malware does not have
capability to change it’s behavior. Nowadays, a malware has a polymorphic and
metamorphic characteristics which means that it will get emerge in a powerful
manner as it will propagate [31]. In such cases, it is impossible to detect a
malware by static analysis which is purely based on comparing the structural
signatures of malware [25].

###  I-B Dynamic Analysis (Behavioral Based)

The dynamic analysis is based on the execution behavior of the program that is
running in a sandbox environment. During it’s execution, some important
metrics are collected such as the system calls a program invokes and file
accessing information. These collected metrics are then compared with a
malware database and it is declared as malware or benign by using certain
threshold levels [26]. Comparatively, this technique works better than the
previous one. Yet, there exists the limitations such as if there is a zero-day
malware (due to which most of the systems are compromised) or a new malware
family for which the reference metrics are missing, the system would not be
able to identify the malware.

###  I-C Heuristic-based

The above two methods are used for the detection of known malware families.
However, the new threats are emerging rapidly and mostly fall into the
category of unknown malware families. The heuristic based methods has
capability to capture the unknown malicious codes. Recent advances in Machine
Learning have created opportunities to use heuristic techniques for the
detection and classification of polymorphic as well as metamorphic malware. In
this technique, when a binary arrives in a system, it is passed through a
Machine Learning (ML) or Deep Learning (DL) pipeline for classification into
malware or benign [27] buckets. This technique has the potential to learn the
internal representations and distributions of Malware data and thus enabling
the system to detect new malware and prevent zero-day attacks.

Inspired by the pros of heuristic based techniques, we have used the similar
in this research. Here, the machine learning pipeline uses a mix of pre-
processing techniques, and machine learning algorithms to predict the
probability of the incoming sequence belonging to a malware family as shown in
figure 1.

Furthermore, we use a benchmark data set, ”Microsoft BIG 15” which was
released by Microsoft on Kaggle[24]. The dataset comprises of nine different
malware families with highly imbalanced classes for the classification
problem. With that, following are the major contributions of this research
work:

  * •

To build up an efficient classifier (both computationally and time-wise) and
in doing so we also aim to analyze the impact of class imbalance problems on
the malware data.

  * •

To compare the results of our proposed method with the existing state of the
art methods which uses the same data set.

  * •

In addition to this, we also report the results on test data by submitting the
proposed model on the Kaggle [24] platform.

##  II Related Work

In the past years, different Machine learning, deep learning, and Natural
language processing based techniques have been used for pre-processing and
identifying the presence of a malware in a network of devices. The process
includes a thorough analysis of the files sent to the devices and identifying
whether the file is malware or benign.

The authors [19] have used sequential pattern mining technology to detect the
most frequent pattern of the opcode sequence for malware identification. In an
other similar research [20], to identify and categorize a malware, the
behavior sequence chain of some malware families has been constructed, and the
similarity between the behavior sequence chain and the sequence of the target
process has been calculated.

Signature based methods are being used in past for detecting malware in
various devices such as mobile [28]. The signature-based method is usually
involved in detecting already present malware from the past and it includes
matching the most similar structure with a known malware but it is
insufficient to detect malware that has remained undetected in the past. The
process includes first extracting the features from malware of a known family
and putting it in a database which is also known as a signature database, so
whenever a new program needs to be classified as malware or benign first the
features are extracted and then cross-checked and compared from the existing
signature database [21][22].

Antivirus vendors have been using signature-based detection methodology for
years, and it is very quick and effective in detecting known threats. This
method is commonly employed to detect malware from the same family. However,
it has limitations for detecting newer malware that employs obfuscation and
polymorphic methods [23].

Malware binary, usually with a file name extension of “.exe” or “.bin”, is a
malicious program that could harm computer operating systems. Sometimes, it
may have many variations with highly reused basic patterns. This implies that
malware binaries could be categorized into multiple families (classes), and
each variation inherits the characteristics of its own family. Therefore, it
is important to effectively detect malware binary and recognize possible
variations. A malware binary file can be visualized to a gray image.

##  III Data Exploration

In this section, we present the patterns and behaviors of a different malware
families using exploratory data analysis techniques.

###  III-A Exploratory Data Analysis

####  III-A1 Class-wise Data Samples

In this part we have performed an analysis on the dataset as a whole. Each
family has an opcode sequence and each malware file has an Id as a 20
character hash value that uniquely identifies the file. Fig. 2 represents the
number of files amongst different families.

Figure 2: Number of samples per class plot.

####  III-A2 Opcode-based relation

Here, we have calculated the frequency of different opcodes and tried to find
a relationship between them. These relations can be linear or non-linear. It
is observed that in any class if an opcode $push$ occurs so $pop$ also has the
same probability of occurring. It is represented using a scatter plot as shown
in Fig. 3 and Fig. 4.

Figure 3: Scatter plot of pop vs push occurrence in the samples.

Figure 4: Scatter plot of mov vs pop occurrence in the samples.

####  III-A3 Embedding Clusters

In this part, we explore the similarities and differences between clusters of
different families that are formed from the learned embedding of the malware
opcodes. We train our model on the training corpus, with the opcodes as
tokens. Each opcode is treated as a word. Then the document embedding vector
is formed using the average of tokens in that document. Fig. 5 and Fig. 6
shows the plots of Principal Component Analysis (PCA) [9] components and t-SNE
[8] components of those 100 dimension document embedding.

Figure 5: Pricipal Component Analysis components of document embedding.

Figure 6: t-SNE components of document embedding.

Although the plots does not show us the totally separable clusters among other
classes, One of the possible reasons is that these are the 2-dimensional
projections of 100-dimensional vectors therefore, there are high chances of
some information loss while projecting these vectors.

###  III-B Correlation of opcode among each file sample

Here, we found the pearson correlation [10] between different opcodes amongst
each other as shown in (Fig. 7). It is quite helpful as it let us know the
strength of opcodes’ relation with each other. the darker the color the
highest the positive correlation.

Figure 7: Correlation of opcodes between each other.

##  IV Methodology

Figure 8 depicts the proposed methodology, the goal is to make an efficient
(in terms of time and computing) yet effective classifier for malware
classification. Performance in terms of metrics and the time efficiency in
real time scenarios are often inverse to each other. Here, for each part, we
define a rule of thumb: 1) Use suitable method at each stage. 2) Use an
efficient alternate of that method to decrease the time consumption. Secondly,
a classifier which can easily identify the hidden patterns without being
biased towards the majority class.

Figure 8: Sequential Embedding-based Attentive Classifier: Our proposed Model

###  IV-A Window Sliding

Window sliding propagates from left to right in a sequential manner. However,
at a particular step at a time stamp $t$ it considers only the information
contained in a window.
For each of the window, the tokens at time stamp $t$ are passed to a network.
It learns embedding for each of the contained opcodes and the window passes
the sequence covering and learning all the opcodes embedding.

###  IV-B Learning Opcodes Embedding

We focus majorly on learning two things: 1) Context 2) Semantics. Both might
be the same in some contexts but here by “Context” we mean to capture the
information of opcode i.e., in what particular context it is being used. For
example, as we know opcode $push$ is being used for pushing to stack but the
same $push$, when used with some $register$, can lead to intentionally
generating an error. So by context, we target to capture this contextual
information.
While the semantics means to know the information or to learn a special
representation vector regarding that particular opcode, so based on our
previous example, we will be learning a vector representation that represents
the opcode $push$.
The sliding window technique is being used to consider $n$ continuous sequence
at a time $(n=5)$. At time $t$ we consider the sequence from $k$ to $k+n$,
where, $k$ is the starting point. The considered sequence in a particular
window is passed through a hidden layer to learn its hidden representations.
This process is continued till it reaches the end of the sequence.
During testing, same sliding window method is used and the learnt parameters
of the hidden layer are used to generate embedding for each opcode in the
sequence. The $d$ dimensional embedding is averaged and complete embedding of
whole sequence is also obtained.

###  IV-C Sequential Block

The embedding learned in the previous part is fed into sequential blocks. As
we already know, those learned embedding preserve both the contextual and the
semantic information which means that we can learn the information contained
in longer sequences. As our data is already massive so there will be many long
term dependencies between opcodes, maybe some starting code has some linkage
with the conclusive code. In order to capture such dependencies, we have some
state of the art architectures already published, one of them is Long Short
Term Memory (LSTM) [1]. Each sequence embedding is passed to multiple LSTM
cells.
On very basic level, the embedding are passed through the input gate, in order
to discard the irrelevant information contained in each of the embedding. We
are going to to obtain a probability in range of $0$ to $1$ indicating the
importance of each scalar value ($0$ being least important and $1$ being most
important) using this, the irrelevant information is discarded. It is done by
applying a $sigmoid$ function activation on the input embedding, this is
having a scalar level importance score. Now, this vector is $multiplied$ with
the output of the last cell to preserve only the important information. Next,
we are going to learn the information in longer term dependencies. The
obtained vector in the last part having only the useful information can be
further used to add new information to it. The new information can be added in
such a way that a single scalar value can either be very important or not
important at all, for that we are going to apply $tanh$ activation function to
obtain a value from range $-1$ to $+1$ indicating how much each scalar is
important for the sequence. This vector is multiplied with input vector and
the scalars contained in it are now representing new information and a new
vector $C_{t}$ is obtained. This new vector can be added to the useful
information preserved vector obtained through the forget gate. Finally the
output is passed through an output gate and the learnt vector is first passed
to the next cell in the form of $C_{t}$ and it is passed as a next state after
going through $tanh$ activation function.
Once the embedding are passed through the process, the dimensions will be
reduced while preserving and capturing long term dependencies.
As we know data is big and sequence contains many information, the model can’t
capture all of the information through single LSTM process.In order to add
more diversity and capture those dependencies which are ignored by LSTM, we
require another architecture in parallel. Here, we can add another LSTM but
that would be a compromise on time efficiency as each of the cell have 3
gates, if we reduce the number of gates, we can have an efficient alternative.
Luckily, we have GRU [2] already published having 2 gates only while
performing same as LSTM in our case.

###  IV-D Averaging Mean Values

From the previous part, two different feature vectors collected. The motive of
both of them is to represent the information regarding opcodes representation,
long-term relations and hidden contextual patterns contained in the sequences.
Now in order to get a single feature vector which benefits from the divergence
of two feature vectors, we cascade their outputs to the same dimensions, each
of the scalar in that dimension is supposed to represent the same part of the
data.
For example, if scalar $i$ in feature vector at position $j$ produced by GRU
represents the data at location $k$ to $k+m$. ($k$ is the starting position
and $m$ is the continuous sequential range) then it is expected that scalar
$h$ produced by LSTM at feature vector’s position $j$ will also represent the
same information of data at location $k$ to $k+m$.
Using this method we got the same and diverse representation of the same
space. In order to greatly benefit from this divergence, we combine them in
such a way that their means converge to the center while eliminating the
extreme values. Taking a simple mean of these will simply center their mean
and this will try to eliminate or at least smooth the extreme values.

###  IV-E Attention Block

This is one of the most important part of our model’s architecture where we
focus more on what’s important and forget the not-so-important information.
Suppose a sequence having 30,000 opcodes and assume that $98\%$ of the code is
normal and the rest of $2\%$ is suspicious or something which can produce a
mal-function in the code. Through our previous block (Sequential Block) we
have the patterns information for complete sequence which includes $98\%$ of
good code and $2\%$ of bad code. Now here is the question: is that $98\%$
useful for us in differentiating the particular sequence from other sequences?
Obviously that $2\%$ is more important for us but as that $2\%$ is in minority
it can easily get ignored and our model ends up focusing on the majority part
of the code. To tackle this issue, we have to assign very little attention to
$98\%$ and a high amount of attention to that $2\%$.
To achieve our motive, we are here using a self-attention mechanism which is
comparatively very light weight and captures the important details. The reason
we did not use the multi-head attention is that it can lead to very expensive
computations hence killing the motive of accurate but efficient classifiers.
Taking input from the sequential block, the attention block first prepares an
attention vector by passing the embedding through hidden layers and applying
$tanh$ activation function to it. The reason we applied the $tanh$ is that it
will give outputs in the range $-1$ and $+1$ which represents the importance
of a particular vector in $n-dimensional$ space. The dot product is taken
between the attention vector and the extracted embedding. It is passed to the
output layer to obtain the final outputs.
However, the proposed model architecture is finished here but it can also be
combined with other models easily using the average means method that we have
described above.

##  V Results

This section includes the details of experimental setup and the results of our
proposed approach. Moreover, it also adds the comparison of the proposed
approach with the state of the art approaches.

Figure 9: Time taken by each part of the network.

###  V-A Experimental Setup

In order to perform the experiments, a pipeline is created as following:

  * •

An executable (.exe) file is transferred as an input to the framework.

  * •

A reverse engineering process is used to convert the .exe file to an assembly
(.asm) file.

  * •

The opcodes are extracted from .asm file without disturbing their sequence of
opcodes.

  * •

At the end, the opcodes are given as an input to the proposed classifier which
will construct embeddings and later on use them for classification.

In this research, we have used tesla T4 GPU with 16GB RAM to train the model.

###  V-B Comparative results with other classifiers

For validation of the classifier, our primary metric is log-loss on 5-fold
validation data and log-loss on Kaggle’s testing data.
For the 5 fold validation, the model was trained on 50 epochs on each fold.
For all the folds, the loss was between 0.065 and 0.045 (can be seen in Fig
10) which is very similar to the existing methods while being the efficient
one.

Figure 10: Logloss on 5-fold validation.

Other than logloss, we also computed some other metrics such as ROC-AUC score,
accuracy and f1 score and compared them with existing methods as shown in I.

TABLE I: Comparison of SEA with state-of-the-art models Models | Metrics
---|---
| Accuracy | F1 Score | Logloss(Kaggle)
SEA | 0.9912 | 0.9908 | 0.0431
(Proposed Method) |  |  |
Hierarchical CNN [5] | 0.9913 | 0.9830 | 0.0419
CNN with | 0.9828 | 0.9830 | 0.0750
Structural Entropy [6] |  |  |
HYDRA [3] | 0.9975 | 0.9951 | -
Orthrus [4] | 0.9924 | 0.9872 | -
Bytes Files Classifier [7] | 0.9861 | 0.9719 | 0.3677

Figure 11 shows the confusion matrix based upon 5-folds cross validation. It
shows that, our model is successful in capturing the patterns of minority
classes as well. Even though only a few examples of minority classes were
given, still our classifier successfully identified the majority of them. It
indicates that our proposed method can give proper attention to hidden
patterns among different classes.

Figure 11: Confusion Matrix based on 5-folds

###  V-C Time efficiency

Another metric for our evaluation was the time efficiency of our classifier
which makes it the fastest classifier among other methods.We have divided our
model’s time into different parts which helps us to understand the time
consumption of each component.

####  V-C1 Training Time

Embedding Training

  * •

100 epochs took 16 minutes 46 secs with 1223 samples on Tesla GPU Neural
Network Training

  * •

250 Epochs in 25 minutes on Tesla GPU

  * •

250/25 = 10 Epochs per minute

####  V-C2 Time efficiency in real time

It took our model 686 seconds to process 11K massive data in 5 folds and each
of the fold took 137.2 seconds. It means that Each sample took 0.0124 seconds.

  * •

0.05 Seconds in embedding extraction

  * •

0.005 seconds in machine learning helping models

  * •

0.001 seconds in neural network part

The time taken in real time scenarios that is for test data on each part of
the network is shown in Fig 9.

TABLE II: Comparison of trainable parameters of different methods with efficient proposed method. Models | Trainable Parameters
---|---
SEA | 812,030
Hierarchical CNN [5] | 20,863,302
CNN with Structural Entropy [6] | 900,148
HYDRA [3] | 2,666,617
Orthrus [4] | 973,989

In table II, we compared trainable parameters of our proposed method with
other methods. The reported numbers show that our proposed method has less
number of trainable parameters. It makes SEA a suitable choice to be used in
light weight smart applications. Although, in terms of logloss, hierarchical
CNN was slightly better than our method and in terms of f1-score, HYDRA was
slightly better. However, if we compare the number of parameters we can see a
huge difference, we were able to almost equal them with very less number of
parameters. Moreover, in terms of f1-score and logloss our method has
outperformed most of them. Hence, the proposed method is a perfect balance
between performance and efficiency.

##  VI Conclusion and Future Work

In this work, we proposed a light weight malware detection technique, called
as SEA, for the detection of malicious software in smart applications. By
formulating the problem of detection of malicious executable as a highly
imbalanced multi-class classification, we developed an attention based model
to improve the results on the minority classes as well. Furthermore, after
benchmarking the developed model against state-of-the-art models, we found
that our proposed model is comparatively better not only in terms of accuracy
and log loss but the number of parameters and time consumption as well. The
other factors such as energy consumption, implementation on hardware devices
such as Raspberry Pi or ESP32 and delay in the network are the part of future
work.

## Acknowledgment

This work has been fully supported by Higher Education Commission (HEC) of
Pakistan under the Grant No. NRPU-5946

## References

  * [1] Greff, K., Srivastava, R.K., Koutník, J., Steunebrink, B.R. and Schmidhuber, J., 2016. LSTM: A search space odyssey. IEEE transactions on neural networks and learning systems, 28(10), pp.2222-2232.
  * [2] Dey, R. and Salem, F.M., 2017, August. Gate-variants of gated recurrent unit (GRU) neural networks. In 2017 IEEE 60th international midwest symposium on circuits and systems (MWSCAS) (pp. 1597-1600). IEEE.
  * [3] Daniel Gibert, Carles Mateu, Jordi Planes, HYDRA: A multimodal deep learning framework for malware classification, Computers & Security, Volume 95, 2020, 101873, ISSN 0167-4048, https://doi.org/10.1016/j.cose.2020.101873.
  * [4] D. Gibert, C. Mateu and J. Planes, ”Orthrus: A Bimodal Learning Architecture for Malware Classification,” 2020 International Joint Conference on Neural Networks (IJCNN), 2020, pp. 1-8, doi: 10.1109/IJCNN48605.2020.9206671.
  * [5] D. Gibert, C. Mateu and J. Planes, ”A Hierarchical Convolutional Neural Network for Malware Classification,” 2019 International Joint Conference on Neural Networks (IJCNN), 2019, pp. 1-8, doi: 10.1109/IJCNN.2019.8852469.
  * [6] Gibert, D., Mateu, C., Planes, J. and Vicens, R., 2018, April. Classification of malware by using structural entropy on convolutional neural networks. In Proceedings of the AAAI Conference on Artificial Intelligence (Vol. 32, No. 1).
  * [7] Gibert, D., Mateu, C., Planes, J. (2018). An End-to-End Deep Learning Architecture for Classification of Malware’s Binary Content. In: Kůrková, V., Manolopoulos, Y., Hammer, B., Iliadis, L., Maglogiannis, I. (eds) Artificial Neural Networks and Machine Learning – ICANN 2018. ICANN 2018. Lecture Notes in Computer Science(), vol 11141. Springer, Cham.
  * [8] Van der Maaten, L. and Hinton, G., 2008. Visualizing data using t-SNE. Journal of machine learning research, 9(11).
  * [9] Daffertshofer, A., Lamoth, C.J., Meijer, O.G. and Beek, P.J., 2004. PCA in studying coordination and variability: a tutorial. Clinical biomechanics, 19(4), pp.415-428.
  * [10] Benesty, J., Chen, J., Huang, Y. and Cohen, I., 2009. Pearson correlation coefficient. In Noise reduction in speech processing (pp. 1-4). Springer, Berlin, Heidelberg.
  * [11] Ling, C.X. and Li, C., 1998, August. Data mining for direct marketing: Problems and solutions. In Kdd (Vol. 98, pp. 73-79).
  * [12] Japkowicz, N., 2000, June. The class imbalance problem: Significance and strategies. In Proc. of the Int’l Conf. on Artificial Intelligence (Vol. 56, pp. 111-117).
  * [13] Alrubayyi, H., Goteng, G., Jaber, M. and Kelly, J., 2021. Challenges of Malware Detection in the IoT and a Review of Artificial Immune System Approaches. Journal of Sensor and Actuator Networks, 10(4), p.61
  * [14] Swain, P.H. and Hauska, H., 1977. The decision tree classifier: Design and potential. IEEE Transactions on Geoscience Electronics, 15(3), pp.142-147.
  * [15] Vishwanathan, S.V.M. and Murty, M.N., 2002, May. SSVM: a simple SVM algorithm. In Proceedings of the 2002 International Joint Conference on Neural Networks. IJCNN’02 (Cat. No. 02CH37290) (Vol. 3, pp. 2393-2398). IEEE.
  * [16] Hopfield, J.J., 1988. Artificial neural networks. IEEE Circuits and Devices Magazine, 4(5), pp.3-10.
  * [17] Ronen, R., Radu, M., Feuerstein, C., Yom-Tov, E. and Ahmadi, M., 2018. Microsoft malware classification challenge. arXiv preprint arXiv:1802.10135.
  * [18] Alrubayyi, H., Goteng, G., Jaber, M. and Kelly, J., 2021. Challenges of Malware Detection in the IoT and a Review of Artificial Immune System Approaches. Journal of Sensor and Actuator Networks, 10(4), p.61
  * [19] Nikolopoulos, S.D. and Polenakis, I., 2017. A graph-based model for malware detection and classification using system-call groups. Journal of Computer Virology and Hacking Techniques, 13(1), pp.29-46.
  * [20] Zhang, Y., Huang, Q., Ma, X., Yang, Z., and Jiang, J., 2016, August. Using multi-features and ensemble learning method for imbalanced malware classification. In 2016 IEEE Trustcom/BigDataSE/ISPA (pp. 965-973). IEEE.
  * [21] Yue, S., 2017. Imbalanced malware images classification: a CNN based approach. arXiv preprint arXiv:1708.08042.
  * [22] K. Alzarooni ., 2012, “Malware variant detection,” Ph.D. dissertation, Dept.Comput. Sci., Univ. College London, London, U.K
  * [23] Saxe, J. and Berlin, K., 2015, October. Deep neural network-based malware detection using two-dimensional binary program features. In 2015 10th international conference on malicious and unwanted software (MALWARE) (pp. 11-20). IEEE.
  * [24] Ronen, R., Radu, M., Feuerstein, C., Yom-Tov, E. and Ahmadi, M., 2018. Microsoft malware classification challenge. arXiv preprint arXiv:1802.10135.
  * [25] Iwamoto, K. and Wasaki, K., 2012, November. Malware classification based on extracted api sequences using static analysis. In Proceedings of the Asian Internet Engineeering Conference (pp. 31-38).
  * [26] Carlin, D., O’Kane, P. and Sezer, S., 2017. Dynamic analysis of malware using run-time opcodes. In Data analytics and decision support for cybersecurity (pp. 99-125). Springer, Cham.
  * [27] Liu, L., Wang, B.S., Yu, B. and Zhong, Q.X., 2017. Automatic malware classification and new malware detection using machine learning. Frontiers of Information Technology & Electronic Engineering, 18(9), pp.1336-1347.
  * [28] Venugopal, D. and Hu, G., 2008. Efficient signature based malware detection on mobile devices. Mobile Information Systems, 4(1), pp.33-49.
  * [29] OWASP https://owasp.org/
  * [30] Tommasini, Martino, Martin Rosso, Emmanuele Zambon, Luca Allodi, and Jerry den Hartog. ”Characterizing Building Automation System Attacks and Attackers.” In 2022 IEEE European Symposium on Security and Privacy Workshops (EuroS&PW), pp. 139-149. IEEE, 2022.
  * [31] Popli, Navneet Kaur, and Anup Girdhar. ”Behavioural analysis of recent ransomwares and prediction of future attacks by polymorphic and metamorphic ransomware.” In Computational Intelligence: Theories, Applications and Future Directions-Volume II, pp. 65-80. Springer, Singapore, 2019.