# The Memory Remains (Undetected): Updates from the Second LIGO/Virgo
Gravitational-Wave Transient Catalog

Moritz Hübner<EMAIL_ADDRESS>Paul Lasky<EMAIL_ADDRESS>Eric
Thrane<EMAIL_ADDRESS>School of Physics and Astronomy, Monash
University, Clayton, VIC 3800, Australia  OzGrav: The ARC Centre of Excellence
for Gravitational Wave Discovery, Clayton, VIC 3800, Australia

###### Abstract

The LIGO and Virgo observatories have reported 39 new gravitational-wave
detections during the first part of the third observation run, bringing the
total to 50. Most of these new detections are consistent with binary black-
hole coalescences, making them suitable targets to search for gravitational-
wave memory, a non-linear effect of general relativity. We extend a method
developed in previous publications to analyse these events to determine a
Bayes factor comparing the memory hypothesis to the no-memory hypothesis.
Specifically, we calculate Bayes factors using two waveform models with
higher-order modes that allow us to analyse events with extreme mass ratios
and precessing spins, both of which have not been possible before. Depending
on the waveform model we find a combined $\ln\mathrm{BF}_{\mathrm{mem}}=0.024$
or $\ln\mathrm{BF}_{\mathrm{mem}}=0.049$ in favour of memory. This result is
consistent with recent predictions that indicate $\mathcal{O}(2000)$ binary
black-hole detections will be required to confidently establish the presence
or absence of memory.

##  I Introduction

The gravitational-wave memory effect is a non-oscillatory part of any
gravitational wave. It can be understood as the part of the gravitational wave
that is sourced by previously emitted waves. Memory causes a permanent
distortion in spacetime long after the wave has passed Zel’dovich, Ya. and
Polnarev (1974); Braginsky and Thorne (1987); Demetrios (1991); Thorne (1992);
Blanchet and Damour (1992). Memory effects are not included in most numerical
relativity models and are hence typically not incorporated in gravitational
waveforms of compact binary coalescences. This is because memory appears in in
the $m=0$ modes of the waveform which are difficult to resolve with numerical
relativity simulations Favata (2010); Pollney and Reisswig (2011); Mitman  _et
al._ (2020).

The slow build-up of memory during compact binary coalescences causes low-
frequency contributions to the gravitational-wave signal. Due to the
relatively low sensitivity at lower frequencies
($f\lesssim$20\text{\,}\mathrm{Hz}$$), measuring memory in any individual
compact binary detection with current generation gravitational-wave
observatories is extremely difficult Lasky  _et al._ (2016); Johnson  _et al._
(2019). However, recent advances in modelling and gravitational-wave signal
analysis have made it possible to coherently search for the presence of memory
in an ensemble of gravitational-wave signals Ashton  _et al._ (2019); Romero-
Shaw  _et al._ (2020a); Hübner  _et al._ (2020); Talbot  _et al._ (2018);
Lasky  _et al._ (2016). Previous studies have shown that thousands of
gravitational-wave detections with LIGO/Virgo Harry  _et al._ (2010); Aasi
_et al._ (2015); Abbott  _et al._ (2016); Tse  _et al._ (2019); Martynov  _et
al._ (2016); Buikema  _et al._ (2020); Abbott  _et al._ (2018); Acernese  _et
al._ (2015, 2019); Abbott  _et al._ (2019, 2020a) operating at their design
sensitivity may be required to confidently detect memory, a milestone which is
likely to occur during the LIGO A+/Virgo+ era Hübner  _et al._ (2020); Boersma
_et al._ (2020).

Since memory arises due to explicitly non-linear effects in general
relativity, any measurement of memory can be considered to be a test of the
theory in this regime. In this study, we focus on measuring the presence of
memory in the gravitational-wave signal itself. Beyond the detection of
memory, future studies will then be able to focus on tests of the exact
amplitude and shape of the memory part of the wave, which are motivated by
modified theories of gravity Garfinkle  _et al._ (2017); Hollands  _et al._
(2017); Satishchandran and Wald (2018); Yang and Martynov (2018), cosmology
Bieri  _et al._ (2016, 2017) and possibly cross-checks with waveform models
Khera  _et al._ (2021).

##  II Methods

We follow mostly the description laid out in Ref. Hübner  _et al._ (2020), but
make some adjustments to our waveform models. For a detailed description of
our methods see Refs. Lasky  _et al._ (2016); Hübner  _et al._ (2020). Recent
advances in waveform modelling have made it more feasible to perform lengthy
sampling processes with models that include higher-order modes (see e.g. Refs.
García-Quirós  _et al._ (2020); Pratten  _et al._ (2020); Colleoni  _et al._
(2021); Varma  _et al._ (2019a, b)). Including higher-order mode effects into
the analysis is not just necessary to avoid systematic errors in the waveform
models that could affect the analysis, they are also required to break a
degeneracy that leaves the sign of the memory ambiguous Lasky  _et al._
(2016). This degeneracy arises because memory changes its sign under a
simultaneous $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ rotation of the
polarisation angle $\psi$ and the phase at coalescence $\phi_{c}$. The same
transformation leaves the quadrupolar (leading order) modes unchanged.

For our analysis we use two waveform models. First, the IMRPhenomXHM waveform
model Pratten  _et al._ (2020); García-Quirós  _et al._ (2020); Colleoni  _et
al._ (2021), which is an aligned-spin model and includes several of the most
dominant higher-order modes. Aligned-spins in this context refer to the black-
hole spins being parallel (or antiparallel) to the orbital angular momentum of
the binary. IMRPhenomXHM covers mass ranges of the more extraordinary
gravitational-wave observations such as GW190814 Abbott  _et al._ (2020b),
which has not been possible for memory analyses until now Hübner  _et al._
(2020). Furthermore, IMRPhenomXHM natively provides the time domain
representation of the waveform in the spherical harmonic mode decomposition,
which is required to calculate the memory contribution using the gwmemory
package Talbot  _et al._ (2018). Finally, IMRPhenomXHM is among the most
computationally efficient waveform models that contain higher-order modes
García-Quirós  _et al._ (2020). Additionally, we perform our analysis with the
NRSur7dq4 waveform model Varma  _et al._ (2019b). While NRSur7dq4 only extends
to a mass ratio $q=m_{2}/m_{1}=1/4$, which mainly excludes GW190412, GW190814,
and parts of the posterior distributions of a few other events, it does
include spin precession effects, which have been shown to exist in the
population of binaries Abbott  _et al._ (2020c). NRSur7dq4 has been trained on
numerical-relativity simulations to a mass ratio of $q=1/4$, and has been
shown to work well down to $q=1/6$ Varma  _et al._ (2019b). We conservatively
restrict the minimum of the prior to $q=1/4$. Using both IMRPhenomXHM and
NRSur7dq4 provides us with a cross-check on our results and allows us to
deploy at least one model on all events 111There are more waveform models such
as IMRPhenomXPHM that both model extreme mass ratios and precession Pratten
_et al._ (2020), however it currently poses practical issues for memory
calculations. The IMPPhenomXPHM implementation in the LALSuite package LIGO
Scientific Collaboration (2018) only provides the spherical-harmonic mode
decomposition in the frequency domain. Frequency-domain waveforms within
LALSuite, unlike time-domain waveforms, are defined so that the merger is
exactly at the beginning/end of the data segment. Calculating memory, which
monotonically rises throughout the binary’s history, thus creates a
discontinuity at the merger time which can not fully be remedied by techniques
such as signal windowing. Furthermore, frequency-domain waveforms are cut-off
at $20\text{\,}\mathrm{Hz}$ which means that we can not perfectly restore the
time domain waveform when we apply an inverse Fourier transform. .

Although precession effects are present at the population level, we do not
expect these to cause substantially different memory estimates. None of the
individual binary black-hole mergers reported so far show significant signs of
precession, though, GW190412 and GW190521 exhibit a mild preference for
precession Abbott  _et al._ (2019, 2020a, 2020). Additionally, we find in
injection studies that precessing spins do not meaningfully change the signal-
to-noise ratio of the memory part of the waveform.

Recent studies of GW190521 suggest that it may be a highly eccentric ($e\geq
0.2$) binary Romero-Shaw  _et al._ (2020b); Gayathri  _et al._ (2020); Abbott
_et al._ (2020). Eccentric binaries may become an interesting target to
measure memory in future generation detectors Favata (2011), however, since
they are not yet firmly established, we assume all binaries to be circular for
this work.

We measure memory by performing model comparison on each event by calculating
a Bayes factor ($\mathrm{BF}_{\mathrm{mem}}$) for the presence of memory. For
a more in-depth discussion of these methods we refer to our previous work
Lasky  _et al._ (2016); Hübner  _et al._ (2020). The memory Bayes factor can
be understood to be the fraction of evidences $Z_{\mathrm{osc+mem}}$ and
$Z_{\mathrm{osc}}$ for the combined oscillatory plus memory waveform and
oscillatory-only waveform, respectively

$\mathrm{BF}_{\mathrm{mem}}=Z_{\mathrm{osc+mem}}/Z_{\mathrm{osc}}\,.$ (1)

In order to calculate the Bayes factors in practise, we perform initial
sampling runs with the dynesty Speagle (2020) implementation within bilby
Ashton  _et al._ (2019); Romero-Shaw  _et al._ (2020a) using the IMRPhenomXHM
and NRSur7dq4 waveform models. We then use a modified version of those models
which contains the memory contribution and obtain the Bayes factor for each
event using importance sampling Payne  _et al._ (2019); Thrane and Talbot
(2019); Hübner  _et al._ (2020). In order to calculate the memory contribution
we adapt the gwmemory Python package Talbot  _et al._ (2018) to support
IMRPhenomXHM and NRSur7dq4 and make additional modifications produce waveforms
with a consistent length which is required for practical inference tasks.
Next, we obtain the memory Bayes factor by summing the memory weights
$w_{\mathrm{mem}}$ over all $n_{\mathrm{post}}$ posterior samples $\theta_{i}$

$\mathrm{BF}_{\mathrm{mem}}=\frac{1}{n_{\mathrm{post}}}\sum_{i=1}^{n_{\mathrm{post}}}w_{\mathrm{mem}}(\theta_{i})\,.$
(2)

The weights are defined to be the ratios between the likelihoods $L$ of the
two hypotheses

$w_{\mathrm{mem}}(\theta_{i})=\frac{L_{\mathrm{mem+osc}}(\theta_{i})}{L_{\mathrm{osc}}(\theta_{i})}\,.$
(3)

Using importance sampling effectively suppresses stochastic sampling noise in
the evidence calculation, making it more suitable than performing inference,
e.g. using nested sampling Skilling (2004), with both hypotheses separately
and comparing the resulting evidence values Hübner  _et al._ (2020).

Since including memory does not add any additional parameters to our problem,
we do not need to consider effects due to increased prior volume. We set our
prior odds on the presence of memory to be 1, i.e. we give equal weight to
memory either being present or not. This means we could also interpret the
memory Bayes factor as an odds.

##  III Results

Figure 1: Memory Bayes factors obtained for the first and second
gravitational-wave transient catalog using IMRPhenomXHM and NRSur7dq4 as
fiducial waveform models, as well as the results from Hübner  _et al._ (2020)
which used NRHybSur3dq8. We abbreviate the last six digits of the regular
event nomenclature with a simplified ‘A’ or ‘B’ to refer the first or second
event detected on that date for readability. Overall, there is no evidence for
or against the memory hypothesis.

As discussed in Ref Hübner  _et al._ (2020), $\ln\mathrm{BF}_{\mathrm{mem}}>8$
can be considered to be very strong evidence for the presence of
gravitational-wave memory. We present the new findings for the unambiguous 36
new binary-black hole observations plus GW190814 additionally to a re-analysis
of the first ten binary-black hole observations in Fig. 1. We calculate a
cumulative $\ln\mathrm{BF}_{\mathrm{mem}}=0.025$ using IMRPhenomXHM and
$\ln\mathrm{BF}_{\mathrm{mem}}=0.049$ using NRSur7dq4, which indicate that
there is no strong evidence favouring or disfavouring the presence of memory
in the signals.

We omit the analysis of GW190425 and GW190426_152155, which are most likely to
be a binary neutron star and a black hole-neutron star binary system,
respectively, for two reasons. First, low mass binaries produce far less
memory within the LIGO band than heavier binary black holes, which makes them
less useful for memory studies Johnson  _et al._ (2019). Second, neither
IMRPhenomXHM nor NRSur7dq4 model neutron star physics, and we would thus need
to implement and test another waveform model for very marginal benefit.
Additionally, GW190426_152155 has a relatively high false alarm rate of 1.4
per year and thus might not be of astrophysical origin. Since they are the
other most likely events to not be of astrophysical origin we also exclude
GW190719_215514 and GW190909_114149. We are otherwise liberal and include all
events that have been reported in the catalogs so far. For the NRSur7dq4 runs,
we exclude GW151012, GW190412, GW190814, GW190513_205428, GW190707_093326,
GW190728_064510, GW190924_021846, and GW190929_012149 since they show
substantial posterior support for $q<1/4$.

We visually verify that our parameter estimates are broadly consistent with
what has been reported in Abbott  _et al._ (2020a). While we occasionally find
minor differences, this is likely due to the fact that some of the runs in the
catalog Abbott  _et al._ (2020a) were performed using a sampling frequency of
$512\text{\,}\mathrm{Hz}$, which implies that no physics beyond
$256\text{\,}\mathrm{Hz}$ was included in the analysis. However, most events
have some contributions at higher frequencies due to higher-order modes, which
cause the posterior to shift.

While we do not find evidence for or against the presence of memory in any of
the observed systems, some of the obtained Bayes factors stand out and some
events deserve our attention.

###  III.1 GW190521

GW190521 is the highest mass event that has been reported so far Abbott  _et
al._ (2020). Even though higher mass systems should in principle create more
memory since they radiate off more energy in gravitational waves, their memory
is shifted outside of the observable LIGO/Virgo band Johnson  _et al._ (2019).
This is consistent with our finding that the measured memory Bayes factor
indicates GW190521 to be uninformative about memory.

###  III.2 GW190521_074359 and similar events

Figure 2: The distribution of Bayes factors for a simulated population of
events from Ref. Hübner  _et al._ (2020) exhibits wide tails. For
visualization, we do not show the two most extreme Bayes factors, which have
values of $\ln\mathrm{BF}_{\mathrm{mem}}=2.4$ and $4.0$.

GW190521_074359 is a near equal mass binary with a total mass of around
$75M_{\odot}$ and is also one of the loudest events observed so far
($\mathrm{SNR}\approx 26$) Abbott  _et al._ (2020a). These properties
generally point towards it being a favourable event with which to measure
memory. Despite this, we find $\ln\mathrm{BF}_{\mathrm{mem}}=-0.12$ using
IMRPhenomXHM and $\ln\mathrm{BF}_{\mathrm{mem}}=-0.05$ using NRSur7dq4, which
is the lowest memory Bayes factor but the highest by absolute value for both
waveform models. A negative $\ln\mathrm{BF}_{\mathrm{mem}}$ for any individual
event is not concerning since they are expected to arise from noise
fluctuations. To show this, we re-examine the population study in Ref. Hübner
_et al._ (2020). While the population study is not perfectly comparable with
the set of the actual measured events as the former is based on a point
estimate of the inferred population of GWTC-1 Abbott  _et al._ (2019), and
used different waveform models, it can still provide us with a cross-check to
see if the observed distribution of Bayes factors is sensible. We find that
out of 2000 simulated events, 28 have $\ln\mathrm{BF}_{\mathrm{mem}}<-0.12$
despite memory being present, which indicates that our measurement
GW190521_074359 is broadly consistent with our expectations. Furthermore, as
we show in Fig. 2, the distribution of Bayes factors has wide tails, which
means that single outlying values are to be expected.

How is it then that the event that looks most likely to contain measurable
memory returns the lowest memory Bayes factor in the catalog? Events that are
highly unfavourable to measure memory with will return
$\ln\mathrm{BF}_{\mathrm{mem}}\approx 0$ as they can only be uninformative. On
the other hand, events like GW190521_074359 are more informative about memory,
but memory is still weak relative to detector noise. Hence, their signals are
also prone to noise fluctuations that may randomly cancel out the memory
contributions, which results in a negative log Bayes factor. This is only true
for weak memory signals, though. If the signal-to-noise ratio of the memory is
greater than one, memory dominates over noise effects and it becomes much less
likely that we measure a negative log Bayes factor due to noise fluctuations.
As Fig. 2 demonstrates, the log Bayes factor of GW190521_074359 is still in
the regime in which we expect noise fluctuations to be able to change the
overall sign of the result.

Other events with high absolute memory log Bayes factors (e.g.
GW190630_185205, GW190828_065509, GW190910_112807) follow a similar pattern to
GW190521_074359 in that they are relatively high signal-to-noise ratio events,
close to equal mass, and have a total mass between $50-80M_{\odot}$.

While most of the differences between the Bayes factors from our two waveform
models are minor, they do appreciably diverge for GW190521_074359. The
difference is unlikely to be due to stochastic sampling noise as this is
strongly suppressed in the importance sampling step Hübner  _et al._ (2020).
In order to understand this difference, we examine the posteriors of both
IMRPhenomXHM and NRSur7dq4. In Fig. 3, we display the posterior as a contour
plot in terms of the obtained memory log weights $\ln w_{\mathrm{mem}}$, and
the inclination angle $\theta_{jn}$ as well as the luminosity distance
$d_{L}$. The observed memory strain $h_{\mathrm{mem}}$ is highly sensitive to
inclination angle $\theta_{jn}$ we are viewing the binary at

$h_{\mathrm{mem}}(\theta_{jn})\propto\sin^{2}\theta_{jn}(17+\cos^{2}\theta_{jn})\,.$
(4)

See e.g. Ref. Favata (2009) for a detailed derivation of this relation. Thus,
memory is most easily seen edge-on ($\theta_{jn}=\pi/2$) as opposed to the
oscillatory part which is preferably emitted face-on ($\theta_{jn}=0,\pi$). In
the posteriors, NRSur7dq4 has stronger support to be closer to face-on whereas
IMRPhenomXHM shows support for GW190521_074359 being an edge-on binary. This
leads to the IMRPhenomXHM weights obtaining larger absolute log weights. The
preference of NRSur7dq4 being closer to face-on conversely corresponds to a
higher inferred luminosity distance than IMRPhenomXHM. Overall, these inferred
differences in posteriors are expected due to systematic differences in the
waveform models, e.g. NRSur7dq4 contains precession effects and all modes up
to $(\ell,|m|)=(4,4)$, whereas IMRPhenomXHM only has aligned spins and modes
$(\ell,|m|)=(2,2),(2,1),(3,3),(3,2),(4,4)$.

###  III.3 Comparison with GWTC-1 analysis

Figure 3: Two dimensional contour plots of the inferred posterior in terms of
inclination and luminosity distance versus the calculated memory weights for
IMRPhenomXHM (blue) and NRSur7dq4 (orange). We show inclination
($\theta_{jn}$) in the top and luminosity distance ($d_{L}$) in the bottom
subfigure. The top subfigure demonstrates why the IMRPhenomXHM have samples
with on average larger weights. NRSur7dq4 samples are somewhat further
constrained away from $\theta_{jn}=\pi/2$ which corresponds to an edge-on
binary for which observed memory is maximal. The bottom subfigure demonstrates
that samples at closer distances have larger absolute log weights on average,
corresponding to the fact that closer events are more informative than ones
further away.

As part of our analysis we redo the analysis in Ref. Hübner  _et al._ (2020)
for the events of the first two observing runs in which we originally used the
hybridized surrogate model NRHybSur3dq8 Varma  _et al._ (2019a). We find that
the difference between NRHybSur3dq8 and either IMRPhenomXHM and NRSur7dq4 to
be in the same order of magnitude as differences between IMRPhenomXHM and
NRSur7dq4. Again, this is most likely to be due to systematic differences in
the waveform models, as well as, to some extent, sampling noise causing slight
deviations. We also note that stochastic sampling error scales with the square
root of the number of events so an error of $10^{-2}$ per event would only
scale up to an error of $\Delta(\ln\mathrm{BF})\approx 0.5$ for the
$\mathcal{O}(2000)$ events required to reach $\ln\mathrm{BF}=8$.

##  IV Conclusion and Outlook

We implement the memory waveforms associated with two waveform models,
IMRPhenomXHM and NRSur7dq4, using the memory calculation method laid out in
Refs. Favata (2010); Talbot  _et al._ (2018) We perform Bayesian model
comparison to search for the presence of memory in the data. Using the
IMRPhenomXHM (NRSur7dq4) model we find a combined
$\ln\mathrm{BF}_{\mathrm{mem}}=0.025\,(0.049)$ in the first and second
gravitational-wave transient catalog. This is consistent with our expectation
that $\mathcal{O}(2000)$ events are required to reach
$\ln\mathrm{BF}_{\mathrm{mem}}=8$, which we consider to be very strong
evidence Hübner  _et al._ (2020); Boersma  _et al._ (2020). We find that
differences in the Bayes factors for each event are likely due to systematic
differences in the waveforms and to a lesser extent due to stochastic sampling
noise.

We have shown that our approach outlined in our previous paper (Ref. Hübner
_et al._ (2020)) is scalable up to a large number of events, demonstrating the
possibility to coherently search for memory in the future. Given the rapid
developments in the waveform community and innovations such as massively
parallel Bayesian inference Smith  _et al._ (2020), we anticipate that more
advanced waveform models can be used for inference in the future. These
waveform models may allow us to calculate higher-order and precessing effects
at greater mass ratios and thus remove the need for using multiple waveform
models to cover all observed events.

##  V Acknowledgements

We would like to thank Ethan Payne, Cecilio Garcia-Quiros, Colm Talbot, and
Bernard Whiting for helpful discussions. This work is supported through
Australian Research Council (ARC) Centre of Excellence CE170100004. PDL is
supported through ARC Future Fellowship FT160100112 and ARC Discovery Project
DP180103155. This material is based upon work supported by NSF’s LIGO
Laboratory which is a major facility fully funded by the National Science
Foundation. This work was performed on the OzSTAR national facility at
Swinburne University of Technology. The OzSTAR program receives funding in
part from the Astronomy National Collaborative Research Infrastructure
Strategy (NCRIS) allocation provided by the Australian Government. This is
LIGO Document No. DCC P2100125.

## References

  * Zel’dovich, Ya. and Polnarev (1974) B. Zel’dovich, Ya. and A. G. Polnarev, Soviet Ast. 18, 17 (1974).
  * Braginsky and Thorne (1987) V. B. Braginsky and K. S. Thorne, Nature 327, 123 (1987).
  * Demetrios (1991) C. Demetrios, Phys. Rev. Lett. 67, 1486 (1991).
  * Thorne (1992) K. S. Thorne, Phys. Rev. D 45, 520 (1992).
  * Blanchet and Damour (1992) L. Blanchet and T. Damour, Phys. Rev. D 46, 4304 (1992).
  * Favata (2010) M. Favata, Class. Quantum Grav. 27 (2010).
  * Pollney and Reisswig (2011) D. Pollney and C. Reisswig, ApJ 732 (2011).
  * Mitman  _et al._ (2020) K. Mitman, J. Moxon, M. A. Scheel, S. A. Teukolsky, M. Boyle, N. Deppe, L. E. Kidder, and W. Throwe, Phys. Rev. D 102, 104007 (2020).
  * Lasky  _et al._ (2016) P. D. Lasky _et al._ , Phys. Rev. Lett. 117, 061102 (2016).
  * Johnson  _et al._ (2019) A. D. Johnson _et al._ , Phys. Rev. D 99, 044045 (2019).
  * Ashton  _et al._ (2019) G. Ashton _et al._ , ApJS 241, 27 (2019).
  * Romero-Shaw  _et al._ (2020a) I. M. Romero-Shaw _et al._ , MNRAS (2020a), 10.1093/mnras/staa2850.
  * Hübner  _et al._ (2020) M. Hübner _et al._ , Phys. Rev. D 101, 023011 (2020).
  * Talbot  _et al._ (2018) C. Talbot _et al._ , Phys. Rev. D 98, 064031 (2018).
  * Harry  _et al._ (2010) G. M. Harry _et al._ , Classical and Quantum Gravity 27, 084006 (2010).
  * Aasi  _et al._ (2015) J. Aasi _et al._ , Class. Quantum Grav 32, 074001 (2015).
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 131103 (2016).
  * Tse  _et al._ (2019) M. Tse _et al._ , Phys. Rev. Lett. 123, 231107 (2019).
  * Martynov  _et al._ (2016) D. V. Martynov _et al._ , Phys. Rev. D 93, 112004 (2016).
  * Buikema  _et al._ (2020) A. Buikema _et al._ , Phys. Rev. D 102, 062003 (2020).
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ , Living Reviews in Relativity 21, 3 (2018).
  * Acernese  _et al._ (2015) F. Acernese _et al._ , Classical and Quantum Gravity 32, 024001 (2015).
  * Acernese  _et al._ (2019) F. Acernese _et al._ , Phys. Rev. Lett. 123, 231108 (2019).
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ ,  9, 031040 (2019).
  * Abbott  _et al._ (2020a) R. Abbott _et al._ , (2020a), arXiv:2010.14527 .
  * Boersma  _et al._ (2020) O. M. Boersma _et al._ , Phys. Rev. D 101, 083026 (2020).
  * Garfinkle  _et al._ (2017) D. Garfinkle _et al._ , Classical and Quantum Gravity 34, 145015 (2017).
  * Hollands  _et al._ (2017) S. Hollands _et al._ , Classical and Quantum Gravity 34, 155005 (2017).
  * Satishchandran and Wald (2018) G. Satishchandran and R. M. Wald, Phys. Rev. D 97, 024036 (2018), arXiv:1712.00873 [gr-qc] .
  * Yang and Martynov (2018) H. Yang and D. Martynov, Phys. Rev. Lett. 121, 071102 (2018).
  * Bieri  _et al._ (2016) L. Bieri _et al._ , Phys. Rev. D 94, 064040 (2016).
  * Bieri  _et al._ (2017) L. Bieri _et al._ , Class. Quantum Grav. 34, 215002 (2017).
  * Khera  _et al._ (2021) N. Khera _et al._ , Phys. Rev. D 103, 044012 (2021).
  * García-Quirós  _et al._ (2020) C. García-Quirós _et al._ , Phys. Rev. D 102, 064002 (2020).
  * Pratten  _et al._ (2020) G. Pratten _et al._ , arXiv e-prints (2020), arXiv:2004.06503 [gr-qc] .
  * Colleoni  _et al._ (2021) M. Colleoni _et al._ , Phys. Rev. D 103, 024029 (2021), arXiv:2010.05830 [gr-qc] .
  * Varma  _et al._ (2019a) V. Varma _et al._ , Phys. Rev. D 99, 064045 (2019a).
  * Varma  _et al._ (2019b) V. Varma _et al._ , Physical Review Research 1, 033015 (2019b).
  * Pratten  _et al._ (2020) G. Pratten _et al._ , Phys. Rev. D 102, 064001 (2020).
  * Abbott  _et al._ (2020b) R. Abbott _et al._ , ApJ 896, L44 (2020b).
  * Abbott  _et al._ (2020c) R. Abbott _et al._ , arXiv e-prints (2020c), arXiv:2010.14533 .
  * Note (1) There are more waveform models such as IMRPhenomXPHM that both model extreme mass ratios and precession Pratten  _et al._ (2020), however it currently poses practical issues for memory calculations. The IMPPhenomXPHM implementation in the LALSuite package LIGO Scientific Collaboration (2018) only provides the spherical-harmonic mode decomposition in the frequency domain. Frequency-domain waveforms within LALSuite, unlike time-domain waveforms, are defined so that the merger is exactly at the beginning/end of the data segment. Calculating memory, which monotonically rises throughout the binary’s history, thus creates a discontinuity at the merger time which can not fully be remedied by techniques such as signal windowing. Furthermore, frequency-domain waveforms are cut-off at $20\text{\,}\mathrm{Hz}$ which means that we can not perfectly restore the time domain waveform when we apply an inverse Fourier transform.
  * Abbott  _et al._ (2020) R. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. D 102, 043015 (2020).
  * Romero-Shaw  _et al._ (2020b) I. Romero-Shaw, P. D. Lasky, E. Thrane, and J. Calderón Bustillo, ApJ 903, L5 (2020b).
  * Gayathri  _et al._ (2020) V. Gayathri _et al._ , (2020), arXiv:2009.05461 .
  * Abbott  _et al._ (2020) R. Abbott _et al._ , Phys. Rev. Lett. 125, 101102 (2020).
  * Favata (2011) M. Favata, Phys. Rev. D 84, 124013 (2011).
  * Speagle (2020) J. S. Speagle, MNRAS 493, 3132 (2020).
  * Payne  _et al._ (2019) E. Payne _et al._ , Phys. Rev. D 100, 123017 (2019).
  * Thrane and Talbot (2019) E. Thrane and C. Talbot, PASA 36, e010 (2019).
  * Skilling (2004) J. Skilling, in  _American Institute of Physics Conference Series_, Vol. 735, edited by R. Fischer, R. Preuss, and U. V. Toussaint (2004) pp. 395–405.
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ , ApJ 882, L24 (2019), arXiv:1811.12940 [astro-ph.HE] .
  * Favata (2009) M. Favata, ApJ 696, 20 (2009), 0902.3660 .
  * Smith  _et al._ (2020) R. J. E. Smith _et al._ , MNRAS 498, 4492 (2020).
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite,” free software (GPL) (2018).


# A quadratic linear-parabolic model-based classification to detect epileptic
EEG seizures

Antonio Quintero-Rincón, Carlos D’Giano
Fundación Lucha contra las Enfermedades Neurológicas Infantiles (FLENI)
Buenos Aires, Argentina
<EMAIL_ADDRESS>
&Hadj Batatia
University of Toulouse, IRIT - INPT
Toulouse, France

###### Abstract

The two-point central difference is a common algorithm in biological signal
processing and is particularly useful in analyzing physiological signals. In
this paper, we develop a model-based classification method to detect epileptic
seizures that relies on this algorithm to filter EEG signals. The underlying
idea is to design an EEG filter that enhances the waveform of epileptic
signals. The filtered signal is fitted to a quadratic linear-parabolic model
using the curve fitting technique. The model fitting is assessed using four
statistical parameters, which are used as classification features with a
random forest algorithm to discriminate seizure and non-seizure events. The
proposed method was applied to 66 epochs from the Children Hospital Boston
database. Results show that the method achieves fast and accurate detection of
epileptic seizures, with a 92% sensitivity, 96% specificity, and 94.1%
accuracy.

_K_ eywords Two-point central difference $\cdot$ curve fitting $\cdot$
parabolic curves $\cdot$ Epilepsy $\cdot$ EEG $\cdot$ Random forest

##  1 Introduction

Epilepsy is a neurological disease that affects people of all ages. It is
characterized by unpredictable seizures resulting from the hyperexcitability
of neurons. The electroencephalogram (EEG) is the predominant modality to
study cerebral activity. The analysis of EEG signals is nearly completely
dependent on visual inspection by the physician to quantify or qualify the
morphology of waves, their goal being the identification and classification of
abnormal patterns in order to provide aid for an epilepsy diagnosis.

This study proposes a simple, fast and adaptable method, implementable in
real-time, to help physicians visually inspect EEG signals for epileptic
seizure detection. The proposed method fits within the framework of model-
based classification. It is based on the statistical parameters obtained from
fitting a quadratic parabolic model to specifically filtered and transformed
EEG signals. Precisely, the two-point central difference algorithm is used to
build a filter for the EEG signals. The filtered signal is subsequently
represented using a quadratic linear-parabolic model, using the curve fitting
approach. Four statistical parameters are then calculated to characterize the
model fitting. These parameters are considered as a feature-vector in the
classification stage.

The literature abounds with research work dealing with the detection of
epileptic seizure onset. Various methods and techniques have been proposed for
this purpose. Many existing methods fit within the large framework of the
feature-based machine learning approach.

The two-point central difference algorithm has proven useful in analyzing
physiological signals, with various medical applications reported in the
literature. A typical work has been reported by Frei et al. to detect
epileptiform discharges in muscles in EEG signals [1]. This algorithm is
usually used to estimate the derivative of a function, giving an estimate
valid only over a limited frequency range. Here, our idea is to use this
algorithm to design a filter that enhances the frequency range of EEG signals
corresponding to brain activity characteristic of epileptic seizures. To
single out the related waveform, we transform the filtered signal into a
quadratic form. A quadratic linear-parabolic model is then fitted to the
transformed signal, using the curve fitting technique [2, 3]. The fitting is
assessed using four statistical parameters (weighted sum of squared residuals,
R-square, adjusted R-square, and root mean squared error). In order to show
that our quadratic linear-parabolic model is pertinent to characterize
epileptic seizures, we develop a classification method where the model-fitting
parameters are used as features. We show that these features are good bio-
markers of epileptic seizures.

The curve fitting technique adopted here is widely used in signal processing
in general and in biomedical applications in particular. It has also been
applied a few times in the domain of EEG and epilepsy. Orhan et al. used the
polynomial curve fitting method to fit a probability density function (PDF) of
EEG signals discretized using equal frequency binning. The method has been
applied to epileptic seizure detection using the comparison of PDFs of
seizures and non-seizures. In [4], second-order Fourier curve fitting has been
used to smooth EEG signals in a pre-processing stage of a model-based method
to predict epileptic seizures. A spline curve fitting to EEG data was proposed
in [5] to detect phase cone patterns to characterize epileptic activity. In
[6], linear curve fitting was applied to interictal heartbeat intervals in
order to combine EEG and ECG signals with the purpose of characterizing
modulation patterns in patients with drug-resistant epilepsy. More recently, a
first-order exponential function has been fit to voltage-gated sodium channel
type 2 to model the conductance-voltage curves in neonatal epileptic subjects
[7].

Supervised machine learning methods have widely been used to detect seizures
in EEG, with support vector machine (SVM) and K-Nearest Neighbor (KNN) being
the most popular. SVM classifiers have been used with various features
extracted from spectral and entropy analysis [8], matching pursuit algorithm
[9], tensor discriminant analysis [10], multifractal detrended fluctuation
analysis [11], and cross-bispectrum analysis [12]. While KNN classifiers have
been developed, among many others, with features from the fractal dimension
[13], and non-linear dimension reduction of frequency domain parameters [14].
In addition to these two popular approaches, Acharya et al. have recently been
the first to develop a seizure detection method using a convolutional neural
network [15].

In general, these methods have good performance thanks to the advanced signal
processing techniques they rely on. However, they consequently have high
computational costs [16]. In this work, we combine for the first time the two-
point difference algorithm and the quadratic linear-parabolic model to design
an original model-based statistical classification method to detect epileptic
seizures in EEG. The proposed feature-vector is simple and fast to calculate
from the quadratic model. Our classification method is based on bootstrap-
aggregated (bagged) decision trees. This approach combines results from
several decision trees to overcome the overfitting effect due to the
variability in the training data [17]. Using the different subset of features
(subspace sampling) to build the weak-tree classifiers reduces the training
time [18], see [19, 20, 21, 22] for some recent works in seizure
classification using EEG.

The remainder of the paper is organized as follows. Section 2 describes the
proposed methodology where the linear-parabolic curve fitting based on two-
point central difference is introduced. In Section 3 the proposed methodology
is applied to real EEG signals from patients suffering from epileptic seizures
using the Children Hospital Boston database. Results of classification using
the random forest to detect seizure events are presented. Conclusions,
limitations, advantages, and perspectives for future work are finally reported
in Section 4.

##  2 Methodology

Let $\boldsymbol{X}\in\mathbb{R}^{N\times M}$ denote the matrix of $M$ EEG
signals $\boldsymbol{x}_{m}\in\mathbb{R}^{N\times 1}$ measured simultaneously
on different channels and at $N$ discrete time instants.

The proposed methodology is composed of five stages. The first stage divides
the original signal $\boldsymbol{X}$ into a set of non-overlapping 1-second
segments using a rectangular sliding window
$\boldsymbol{\Omega}=\boldsymbol{\Omega}_{0}\left(w-\frac{W-1}{2}\right)$ with
$0\leq w\leq W-1$, so that
$\boldsymbol{X}[n]=\boldsymbol{\Omega}[n]\boldsymbol{X}$. In the second stage,
the two-point central difference algorithm is used to estimate the
coefficients $b[k]$ of the filter for each 1-second segment
$\boldsymbol{X}[n]$. The signal is then convolved by the resulting filter
$\widetilde{\boldsymbol{X}}[n]=\boldsymbol{X}[n]\ast b[k]$. The third stage
fits a quadratic linear-parabolic model to the filtered signal, and estimates
the associated parameters. The fourth stage calculates four statistical
parameters for each EEG segment, namely the weighted sum of squared residuals
($\zeta$), the R-square value ($\phi$), the adjusted R-square value ($\sigma$)
and the root mean squared error ($\psi$), to assess the quadratic model
fitting. Finally, in stage five, the feature-vector
$\rho=[\zeta,\phi,\sigma,\psi]$ associated with each segment is classified
using a random forest classifier in order to discriminate between seizure and
non-seizure.

The following sections present the two-point central difference algorithm, the
model fitting method based on curve fitting and the associated statistical
parameters, and the random forest classifier.

###  2.1 Two-point central difference algorithm

The principle of the two-point central difference algorithm consists in
subtracting non-adjacent, regularly spaced, pairs of points. The objective is
to extract a slope from $\boldsymbol{X}[n]$:

$\displaystyle
X^{\prime}[n]=\frac{\boldsymbol{X}[n+L]-\boldsymbol{X}[n-L]}{2LT_{s}},$ (1)

where $L$ is called the skip factor that defines the distance between points,
$T_{s}$ is the sample interval scaling. Taking the $Z$-transform of both
sides, we obtain

$\displaystyle\boldsymbol{X}^{\prime}(z)=\frac{z^{L}-z^{-L}}{2LT_{s}}X(z).$
(2)

The transfer function of (2) is given by

$\displaystyle\boldsymbol{Y}(z,L)=\frac{\boldsymbol{X}^{\prime}(z)}{\boldsymbol{X}(z)}=\frac{z^{L}-z^{-L}}{2LT_{s}}.$
(3)

Thus, the frequency response is obtained by replacing $z$ by $e^{j\Omega}$

$\displaystyle\boldsymbol{Y}(e^{j\Omega},L)=\frac{e^{jL\Omega}-e^{-jL\Omega}}{2LTs}=\frac{jsin(L\Omega)}{LT_{s}}.$
(4)

Using a symmetric FIR (finite impulse response) filter, the two-point central
difference algorithm is based on an impulse function containing two
coefficients of equal but opposite sign spaced $L$ points apart with the
following coefficients

$\displaystyle b[k]=\left\\{\begin{matrix}\frac{1}{2LT_{s}}&k=-L\\\
-\frac{1}{2LT_{s}}&k=+L\\\ 0&k\neq\pm L.\\\ \end{matrix}\right.$ (5)

Note that, FIR filters are free from stability problems and they cause no time
delay and no phase distortion within the pass-band.

A convolution operation is estimated between each 1-second segment
$\boldsymbol{X}[n]$ and the coefficients $b[k]$ in order to obtain a new
filtered signal in a bandwidth frequency until 50 Hz:

$\displaystyle\widetilde{\boldsymbol{X}}[n]=\boldsymbol{X}[n]\ast b[k].$ (6)

Note that, this is the effective bandwidth of the filter where the EEG
activity has an important clinical relevance [23]. This allows to
automatically reject the high line-noise artifacts greater than 50 Hz as
$L=F_{s}/5$, as the sampling rate $F_{s}$ is the 256 Hz. Note that $L$ is
rounded. We refer the reader to [24, 25, 26] for a comprehensive treatment of
the mathematical properties of two-point central difference algorithm.

###  2.2 Quadratic linear-parabolic model

We propose to fit a quadratic linear-parabolic model to the filtered signal.
Precisely, the model is fitted to the couple
$\left(\widetilde{\boldsymbol{X}},\widetilde{\boldsymbol{X}}^{2}\right)$,
where $\widetilde{\boldsymbol{X}}$ is the complete filtered signal, obtained
by concatenating the segments given in equation (6). The idea is to fit the
signal with a curve of the form:

$\displaystyle y=a\;sin(x-\pi)+b\;(x-10)^{2}+c.$ (7)

The fitting is performed using the least squares method [27, 28, 29] to
estimate the three parameters $a,b,c$.

###  2.3 Model-fitting statistical parameters

In order to assess the fitting of the curve given in equation (7), the
following four statistical parameters are estimated. Note that $y$ is the
observed data, $\hat{y}$ is the predicted value using the quadratic linear-
parabolic model, and $\bar{y}$ is the mean of the observed data.

_Weighted sum of squared residuals_ ($\zeta$): It is used to measure the total
deviation of the response values from the predicted values, and is defined as

$\displaystyle\zeta=\sum_{i=1}^{n}w_{i}\left(y_{i}-\hat{y}_{i}\right)^{2}.$
(8)

The weights allow taking into consideration the different uncertainties of the
measurements, and are calculated as follows

$\displaystyle w_{i}=\frac{1}{\sigma_{i}^{2}},$ (9)

where $\sigma^{2}$ is the variance.
_R-square_ ($\phi$): It is the square of the correlation between the response
values and the predicted response values, and is defined as

$\displaystyle\phi=1-\sum_{i=1}^{n}\frac{\left(y_{i}-\hat{y}_{i}\right)^{2}}{\left(y_{i}-\bar{y}_{i}\right)^{2}}.$
(10)

_Adjusted R-square_ ($\sigma$): It adjusts the R-square residual degrees of
freedom, and is defined as

$\displaystyle\sigma$ $\displaystyle=1-\frac{(1-\phi)(n-1)}{n-m-1},$ (11)

where $n$ is the number of response values and $m$ is the estimated fitted
coefficients from the response values. A value of $\sigma$ closer to 1
indicates a better fit.
_Root mean squared error_ ($\psi$): It is an estimate of the standard
deviation of the random component in the data, and is defined as

$\displaystyle\psi=\sqrt{\frac{\zeta}{n-m}}.$ (12)

###  2.4 Random Forest classification

In this section, we present a classification method to identify seizures. The
four statistical parameters presented above are used as classification
features. A random forest classification technique is adopted. Random forest
is an ensemble learning technique that combines the Bagging algorithm and the
random subspace method using decision trees as the base classifier. The goal
is to discriminate between seizure and non-seizure events, specifically the
seizure onset through the feature predictor vector
$\rho=[\zeta,\phi,\sigma,\psi]$ associated with each EEG segment
$\widetilde{\boldsymbol{X}}$. The idea is to design a prediction function
$f(\rho)$ to estimate the binary response for two classes,
$\Gamma=\\{\gamma_{1},\gamma_{0}\\}$, where $\gamma_{1}=1$ for a seizure event
and $\gamma_{0}=0$ for a non-seizure event. A loss function
$L(\Gamma,f(\rho))$ determines the prediction function as follows

$\displaystyle E_{\rho,\Gamma}\left(L(\Gamma,f(\rho))\right),$ (13)

where $E_{\rho,\Gamma}$ denotes the expectation with respect to an unknown
joint distribution $P_{\rho,\Gamma}(\rho,\Gamma)$ of $\rho$ and $\Gamma$. For
the classification, the next zero-one loss function is used

$\displaystyle L(\Gamma,f(\rho))=I(\Gamma\neq f(\rho))=\begin{cases}1\text{ if
}&\Gamma=f(\rho)\approx\text{ for }\gamma_{1}\text{ or seizure},\\\ 0\text{ if
}&\Gamma\neq f(\rho)\approx\text{ for }\gamma_{o}\text{ or non-seizure}.\\\
\end{cases}$ (14)

The ensemble learning approach constructs $f$ by using $J$ decision-tree base-
classifiers

$\displaystyle\\{h_{1}(\rho,\Theta_{1}),\cdots,h_{J}(\rho,\Theta_{J})\\}$ (15)

where $\Theta_{j}$, $j=1,\cdots,J$, is an independent collection of random
variables. These base classifiers are combined by majority vote, so that

$\displaystyle
f(\rho)=\operatorname*{arg\,max}_{\gamma\in\Lambda}\sum_{j=1}^{J}I(\gamma=h_{j}(\rho)).$
(16)

Consider a dataset $\mathcal{D}=\\{d_{1},d_{2},..,d_{N}\\}$ with
$d_{i}=(\rho_{i},c_{i})$, where $c_{i}$ is a class label, $\rho$ denotes the
vector predictor and $\gamma_{i}$ the response. For a particular realization
$\theta_{j}$ of $\Theta_{j}$, the fitted tree is denoted
$\hat{h}_{j}(\rho,\theta_{j},\mathcal{D})$. The random forest classifier uses
random subspace in two ways. First, each tree is fitted to an independent
bootstrap sample from the original data. The randomization involved in
bootstrap sampling gives one part of $\Theta_{j}$. Second, the best split is
retained for all $\rho$ predictors according to the randomly selected subset
of $m$ predictor variables.

The randomization used to sample the predictors gives the remaining part of
$\Theta_{j}$. The resulting class predicted is by the majority combination of
unweighted voting of the trees, see algorithm 1. We refer the reader to [17,
30, 18] for a comprehensive treatment of the properties of random forest
classifier.

Data: data set $\mathcal{D}$; predictor vector $\rho$

Result: ensemble of tree models whose predictions are to be combined by voting
or averaging.

for  _$j=1$ to $J$_ do

1:Take a bootstrap sample $\mathcal{D}_{j}$ of size $N$ from $\mathcal{D}$.

2:Using the bootstrap sample $\mathcal{D}_{j}$ as the training data, fit a
tree using binary recursive partitioning:

$a$. Start with all observations in a single node ;

$b$. Repeat the following steps recursively for each unsplit node until the
stopping criterion is met:

$i$. Select $m$ predictors at random from the $\rho$ available predictors. ;

$ii$. Find the best binary split among all binary splits on the $m$ predictors
from Step i.;

$iii$. Split the node into two descendant nodes using the split from Step ii.;

end for

To make a prediction at a new point, use equation (16).

Algorithm 1 Random forest classifier algorithm

##  3 Experimentation and discussion

In this section, we evaluate the proposed methodology using the Children
Hospital Boston database, which presented in the following section.

###  3.1 Data

For the experimentation purpose, we considered a dataset from the Children’s
Hospital Boston database [31, 32], which consists of 22 EEG recordings from
pediatric subjects with intractable seizures. All signals were sampled at 256
Hz with 16-bit resolution by using the International 10-20 system. The set of
recordings lasted on average 35 minutes for 30 subjects in total; 2 hours for
4 subjects; and 12 hours for 2 other subjects. Taken together the recordings
account for 60 hours of EEG recordings and 139 seizures. No distinctions
regarding the types of seizure onsets were considered; the data contains
focal, lateral, and generalized seizure onsets. Furthermore, the recordings
were made in a routine clinical environment, so non-seizure activity and
artifacts such as head/body movement, chewing, blinking, early stages of
sleep, and electrode pops/movements are present in the data. In this work, we
have used 66 epochs from 9 different subjects, see Table LABEL:tab:data. Each
recording contains a seizure event, whose onset time has been labeled by an
expert neurologist. Moreover, for each seizure segment, the neurologist also
selected one adjacent non-seizure signal segment before and after the seizure,
of the same length to represent healthy brain activity. In total, we
considered 33 seizures and 33 non-seizure events. Non-seizure control events
have been chosen just before seizures. This choice is justified for two
reasons. First, we are interested in identifying seizure onset, therefore it
is most important to distinguish its signal from preceding normal activity.
Second, the EEG activity following a seizure remains chaotic for a long time
and comparing it to seizure onset is not pertinent. The selected signals had
the same montage by using 23 scalp EEG channels.

Table 1: Length of the 33 seizures used in this study. The duration has an average of 79,4 sec. Epoch | Seizure | Duration (sec) | Genre | Age (years)
---|---|---|---|---
01 | $01\\_03$ | 40 | F | 11
02 | $01\\_04$ | 27 |  |
03 | $01\\_15$ | 40 |  |
04 | $01\\_18$ | 90 |  |
05 | $01\\_21$ | 115 |  |
06 | $01\\_26$ | 101 |  |
07 | $02\\_16$ | 82 | M | 11
08 | $02\\_16$ | 21 |  |
09 | $02\\_19$ | 9 |  |
10 | $03\\_01$ | 52 | F | 14
11 | $03\\_02$ | 55 |  |
12 | $03\\_03$ | 69 |  |
13 | $03\\_04$ | 52 |  |
14 | $03\\_34$ | 47 |  |
15 | $03\\_35$ | 64 |  |
16 | $03\\_36$ | 67 |  |
17 | $05\\_17$ | 120 | F | 7
18 | $05\\_22$ | 117 |  |
19 | $06\\_18$ | 12 | F | 1.5
20 | $06\\_24$ | 16 |  |
21 | $07\\_12$ | 86 | F | 14.5
22 | $07\\_13$ | 144 |  |
23 | $07\\_19$ | 83 |  |
24 | $08\\_02$ | 189 | M | 3.5
25 | $08\\_05$ | 190 |  |
26 | $08\\_21$ | 338 |  |
27 | $09\\_06$ | 64 | F | 10
28 | $09\\_19$ | 64 |  |
29 | $10\\_20$ | 10 | M | 3
30 | $10\\_27$ | 65 |  |
31 | $10\\_30$ | 62 |  |
32 | $10\\_31$ | 76 |  |
33 | $10\\_89$ | 54 |  |

###  3.2 Model fitting

This section presents the results of the quadratic linear-parabolic model
fitting. Table 2 presents the three coefficients $(a,b,c)$ estimated according
to equation (7). Coefficients for seizures and non-seizures are shown
separately, for illustration purpose only. This allows us to note the
significant difference between values corresponding to the two classes. The
actual quadratic linear-parabolic model is fitted to the entire signal,
without distinction of seizure and non-seizure.

| Non-Seizure | Seizure
---|---|---
Coefficients | Value | CB95% | Value | CB95%
$a$ | -0.6855 | (-4.095, 2.725) | 1.368 | (-6.717, 9.453)
$b$ | 0.9989 | (0.9989, 0.9989) | 0.9998 | (0.9998, 0.9998)
$c$ | -36.05 | (-38.53, -33.57) | -35.45 | (-41.47, -29.43)
Table 2: The quadratic linear-parabolic model coefficients and the associated
95% confidence bounds (CB).

Figure 1.1 and 1.2 show the quadratic linear-parabolic curves obtained for
seizure and non-seizure segments, respectively. One notices that the proposed
model (7) provides a visually precise fit to the EEG waveform. In addition,
the prediction bounds have a small uncertainty throughout the entire data
range, hence new observations should be predicted with high accuracy.

(a) Non-Seizures (b) Seizures

Figure 1: Curves of the quadratic linear-parabolic model for seizure and non-
seizure events. Note that the scale of the seizure curve is much larger than
the one of the non-seizure curve. The x-axis represents
$\widetilde{\boldsymbol{X}}$ while the y-axis represents
$\widetilde{\boldsymbol{X}}^{2}$ (the unit on the X and Y axes correspond to
the amplitude of the signal). The prediction bounds have a small uncertainty
throughout the entire data range, therefore new observations can be predicted
with high accuracy.

###  3.3 Classification

The dataset described in section 3.1 has been used to train the classifier
according to the proposed method. The capacity of the proposed classification
scheme to discriminate between seizure and non-seizure events in order to
detect the seizure onset in EEG signals has been assessed. In Figure 2, one
can notice the great difficulty in discriminating between seizure and non-
seizure in EEG raw data. The start and end of the seizure in this EEG signal
were labeled by the neurologist using two lines. The first line divides the
EEG signal at 81 sec (onset) and the second at 162 sec (offset).

Figure 2: EEG raw example. We can observe the great difficulty in
discriminating between seizure and non-seizure. The onset detection begins in
81 sec according to the medical annotation but by non-expert visual inspection
cannot reach the same conclusion. The y-axis scale for each channel is between
$\pm 200$ mv.

Table 3 reports the four statistical parameters, described in section 2.3,
calculated for all seizure and non-seizure segments. The weighted sum of
squared residuals ($\zeta$) and the root mean squared error ($\psi$) show much
larger values for seizure events with respect to non-seizure events, which
suggests that these features can be used to discriminate between seizure and
non-seizures in EEG signals. While R-square ($\phi$) and adjusted R-square
($\sigma$) have values close to one for both types of events, which suggests
that the model has high accuracy in the fit. This can be visually corroborated
in Figure 1.

Events | $\zeta$ | $\phi$ | $\sigma$ | $\psi$
---|---|---|---|---
Non-Seizure | 3.624e+14 | 0.9996 | 0.9996 | 4839
Seizure | 2.03e+15 | 0.9999 | 0.9999 | 1.145e+04
Table 3: Goodness-of-fit statistics. The weighted sum of squared residuals
($\zeta$) and root mean squared error ($\psi$) show much larger values for
seizure events with respect to non-seizure events. While R-square ($\phi$) and
adjusted R-squared ($\sigma$) have values near one for both, seizure and non-
seizure events. These four parameters show good potential to discriminate
between seizure and non-seizures in EEG signals.

Estimating the four statistical parameters at each 1-second segment allows us
to detect changes between epileptic EEG seizures. When we apply equation (7)
to each 1-second segments, similar results are obtained as those applied to
the signal as a whole. Figure 3 shows the similarity with the previous figure
1. Both seizure and non-seizure events have an excellent fitting with the
linear-parabolic model, suggesting that the scale changes and the different
values of the estimated parameters can be used to differentiate seizure and
non-seizure events in EEG signals.

(a) Non-Seizures (b) Seizures

Figure 3: Quadratic parabolic curve fitting example for one segment
$\widetilde{\boldsymbol{X}}[n]$. We can see the waveform similarity with the
previous figure 1

In order to cope with the high variability of EEG signals, the features were
normalized to be in the range $[0,1]$. Table 4 reports the mean values for the
four statistical parameters $\rho=[\zeta,\phi,\sigma,\psi]$ calculated for
each 1-second segment from $\boldsymbol{X}[n]$. The weighted sum of squared
residuals ($\zeta$) has values closer to zero while R-square ($\phi$) and
adjusted R-square ($\sigma$) have values closer to one. This shows that the
model has a high fit accuracy and can be used for seizure detection, as
corroborated by Figure 3.

Events | $\zeta$ | $\phi$ | $\sigma$ | $\psi$
---|---|---|---|---
Non-Seizure | 0.14 | 0.99 | 0.99 | 0.30
Seizure | 0.11 | 1.00 | 1.00 | 0.25
Table 4: Means of the statistical parameters with values normalized between 0
and 1 for all each 1-second time segments. A value closer to zero in the
weighted sum of squared residuals ($\zeta$) and Root mean squared error
($\psi$) and values closer to one in R-square ($\phi$) and adjusted R-squared
($\sigma$) suggest a good fit quality and therefore can be used as a
predictor.

Figure 4 shows the main scatter plots of couples of the normalized statistical
parameters $(\zeta,\phi,\sigma,\psi)$. It is interesting to note that in
$\zeta$ versus $\psi$ (fig.4.a) and $\psi$ versus $\phi$ (fig.4.c) the seizure
events (red dots) are concentrated closer to one with respect to non-seizure
events (blue dots). In (fig.4.b) $\zeta$ versus $\psi$ has a big tendency to
go to the right in the seizure events (red dots) with respect to non-seizure
events (blue dots). In (fig.4.d) $\phi$ versus $\sigma$ shows a clear
concentration of seizure events (red dots) with respect to non-seizure events
(blue dots) on the upper right side. These scatter plots show that our
predictor-vector $\rho=[\zeta,\phi,\sigma,\psi]$ is potentially useful to
discriminate between seizures and non-seizures events in EEG signals.

(a) $\zeta$ vs. $\phi$ (b) $\zeta$ vs. $\psi$ (c) $\psi$ vs. $\phi$ (d) $\phi$
vs. $\sigma$

Figure 4: Scatter plots for the four statistical parameters
$\rho=[\zeta,\phi,\sigma,\psi]$ observed through $1518$ events from 66 epochs,
33 seizures (red dots) and 33 non-seizures (blue dots) before the seizure. The
weighted sum of squared residuals ($\zeta$), R-square ($\phi$), adjusted
R-squared ($\sigma$) and root mean squared error ($\psi$) shows a clear
concentration of seizure events (red dots) with respect to non-seizure events
(blue dots), which suggest that our model is potentially useful to detect
between seizures and non-seizures events in EEG signals.

To assess the performance of the proposed method, we adopted a supervised
testing approach and used $1518$ events (1-second signals) to train and test
the method with the 20-fold cross-validation technique of the predictor vector
$\rho=[\zeta,\phi,\sigma,\psi]$. As explained in section 3.1, the events are
extracted from 33 seizure epochs and 33 non-seizure epochs, giving a total of
66 epochs from 9 different subjects. The method gives good classification
performance, with 92% true positives rate (TPR) or sensitivity, 96% true
negative rate (TNR) or specificity, 4% false positive rate and 94.1% accuracy.
Please note that segments for training and testing are drawn randomly without
consideration of their epochs. This process is repeated 1000 times to ensure
not-bias in the partitioning.

##  4 Conclusion

This work presented a new method to detect epileptic seizures in EEG signals.
The proposed model-based classification method relies on the design of a
specific filter using the two-point central difference algorithm. A linear-
parabolic model is fitted to the filtered signal using mean squares. Four
statistical parameters associated with the model fitting are used as
classification features to discriminate seizure and non-seizure events. A
random forest classifier has been adopted to evaluate the ability of these
parameters to detect seizures. The proposed methodology was applied on $66$
balanced events from the Children Hospital Boston database. Results suggest
that the proposed algorithm is a powerful tool for detecting epileptic seizure
events in EEG signals achieving a 94.1% accuracy.

The proposed method has two main advantages compared to existing methods.
First, its low computational complexity makes it implementable in quasi-real-
time. Second, the method can be easily extended to work separately with
different brain rhythms, by adapting the skip factor $L$. Despite these
advantages, the method suffers some limitations, mainly noise and robustness.

The signal is affected by noise and artifacts due to the acquisition and pre-
processing. This affects the precision of the results, especially when the
method is generalized to complex epileptic forms. Outliers are not directly
considered in the estimation procedure that needs more robustness.

Future work will focus on robust approaches to consider noise and artifacts.
Optimization techniques will be investigated in order to remove outliers and
improve the accuracy of the detection.

## References

  * [1] M. G. Frei, R. L. Davidchack, I. Osorio, Least squares acceleration filtering for the estimation of signal derivatives and sharpness at extrema and application to biological signals, IEEE Transactions on Biomedical Engineering 46 (8) (1999) 971–977.
  * [2] N. Chernov, C. Lesort, Statistical efficiency of curve fitting algorithms, Computational Statistics and Data Analysis 47 (4) (2004) 713–728.
  * [3] E. López-Rubio, K. Thurnhofer-Hemsi, E. B. Bláquez-Parra, O. D. de Cózar-Macías, M. C. L. de Guevara-Muñoz, A fast robust geometric fitting method for parabolic curves, Pattern Recognition (2018) 301–316.
  * [4] H. Zhang, J. Su, Q. Wang, Y. Liu, L. Good, J. Pascual, Predicting seizure by modeling synaptic plasticity based on EEG signals - a case study of inherited epilepsy, Communications in Nonlinear Science and Numerical Simulation 56 (2018) 330–343.
  * [5] C. Ramon, M. D. Holmes, M. V. Wise, D. Tucker, K. Jenson, S. R. Kinn, Oscillatory patterns of phase cone formations near to epileptic spikes derived from 256-channel scalp EEG data, Computational and Mathematical Methods in Medicine 2018 (2018) 1–15.
  * [6] H. Liu, Z. Yang, F. Meng, L. Huang, W. Qu, H. Hao, J. Zhang, L. Li, Chronic vagus nerve stimulation reverses heart rhythm complexity in patients with drug-resistant epilepsy: An assessment with multiscale entropy analysis, Epilepsy & Behavior 83 (2018) 168–174.
  * [7] A. Begemann, M. A. A. na, M. Zweier, M. V. et. al, Further corroboration of distinct functional features in SCN2A variants causing intellectual disability or epileptic phenotypes, Molecular Medicine 25 (6) (2019) 1–15.
  * [8] S.-F. Liang, H.-C. Wang, W.-L. Chang, Combination of EEG complexity and spectral analysis for epilepsy diagnosis and seizure detection, EURASIP Journal on Advances in Signal Processing 2010 (1) (2010) 853434.
  * [9] T. L. Sorensen, U. L. Olsen, I. Conradsen, J. Henriksen, T. W. Kjaer, C. E. Thomsen, H. B. D. Sorensen, Automatic epileptic seizure onset detection using matching pursuit: a case study, 32nd Annual International Conference of the IEEE EMB 2010 (2010) 3277–3280.
  * [10] S. Nasehi, H. Pourghassem, A novel fast epileptic seizure onset detection algorithm using general tensor discriminant analysis., Journal of Clinical Neurophysiology 30 (4) (2013) 362–370.
  * [11] Z. Zhang, T. Wen, W. Huang, M. Wang, C. Li, Automatic epileptic seizure detection in EEGs using MF-DFA, SVM based on cloud computing, Journal of X-Ray Science and Technology 25 (2) (2017) 261–272.
  * [12] N, Mahmoodian, A. Boese, M. Friebe, J. Haddadnia, Epileptic seizure detection using cross-bispectrum of electroencephalogram signal, Seizure: European Journal of Epilepsy 66 (2019) 4–11.
  * [13] G. Polychronaki, P. Ktonas, S. Gatzonis, A. Siatouni, P. Asvestas, H. Tsekou, D. S. an KS Nikita, Comparison of fractal dimension estimation algorithms for epileptic seizure onset detection., Journal of Neural Engineering 7 (4) (2014) 046007.
  * [14] J. Birjandtalab, M. B. Pouyan, D. Cogan, M. Nourani, J. Harvey, Automated seizure detection using limited-channel EEG and non-linear dimension reduction, Computers in Biology and Medicine 82 (2013) 49–58.
  * [15] U. Acharya, S. L. Oh, Y. Hagiwara, J. Tan, H. Adeli, Deep convolutional neural network for the automated detection and diagnosis of seizure using EEG signals, Computers in Biology and Medicine 17 (2017) 1–9.
  * [16] A. Quintero-Rincón, M. Pereyra, C. D’Giano, M. Risk, H. Batatia, Fast statistical model-based classification of epileptic EEG signals, Biocybernetics and Biomedical Engineering 38 (4) (2018) 877–889.
  * [17] L. Breiman, Random forests, Machine Learning 45 (2001) 5–32.
  * [18] P. Flach, Machine Learning, Cambridge University Press, 2012.
  * [19] F. Manzouri, S. Heller, M. Dumpelmann, P. Woias, A. Schulze-Bonhage, A comparison of machine learning classifiers for energy-efficient implementation of seizure detection, Frontiers in Systems Neuroscience 12 (43) (2018).
  * [20] JE. Le Douget, A. Fouad, M. Maskani Filali, J. Pyrzowski, M. Le Van Quyen, Surface and intracranial EEG spike detection based on discrete wavelet decomposition and random forest classification, Annual International Conference of the IEEE Engineering in Medicine and Biology Society (2017) 475–478.
  * [21] C. Donos, M. Dumpelmann, A. Schulze-Bonhage, Early seizure detection algorithm based on intracranial EEG and random forest classification, International Journal of Neural Systems 25 (5) (2015) 1550023.
  * [22] S. López, G. Suarez, D. Jungreis, I. Obeid, J. Picone, Automated identification of abnormal adult EEGs, IEEE Signal Processing in Medicine and Biology Symposium.
  * [23] S. Sanei, Adaptive Processing of Brain, Wiley, 2013.
  * [24] A. T. Bahill, J. D. McDonald, Frequency limitations and optimal step size for the two-point central difference derivative algorithm with applications to human eye movement data, IEEE Transactions on Biomedical Engineering 30 (3) (1983) 191–194.
  * [25] C.-C. Tseng, S.-L. Lee, Design of digital differentiator using difference formula and formula and richardson extrapolation, IET Signal Processing 2 (2) (2008) 177–188.
  * [26] J. L. Semmlow, B. Griffel, Biosignal and Medical Image Processing, CRC Press, 2014\.
  * [27] P. C. Hansen, V. Pereyra, G. Scherer, Least Squares Data Fitting with Applications, Johns Hopkins University Press, 2013.
  * [28] C. Neuhauser, Calculus For Biology and Medicine, Pearson, 2014.
  * [29] G. Samarasinghe, A. Sowmya, D. A. Moses, A semi-quantitative analysis model with parabolic modelling for DCE-MRI sequences of prostate, International Conference on Digital Image Computing: Techniques and Applications (DICTA).
  * [30] C. Zhang, Y. Ma, Ensemble Machine Learning. Methods and Applications, Springer-Verlag New York, 2012.
  * [31] A. Shoeb, H. Edwards, J. Connolly, B. Bourgeois, S. T. Treves, J. Guttagf, Patient-specific seizure onset detection, Epilepsy and Behavior 5 (2004) 483–498.
  * [32] A. Goldberger, L. Amaral, L. Glass, J. Hausdorff, P. Ivanov, R. Mark, J. Mietus, G. Moody, C.-K. Peng, H. Stanley, Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals, Circulation 101(23): 101 (23) (2000) e215–e220.