Decoding Epileptogenesis in a Reduced State Space

1.3 MB·

by  François G. Meyer

Decoding Epileptogenesis in a Reduced State Space
Franc¸ois G. Meyer1, Alexander M. Benison2, Zachariah Smith2, and Daniel S. Barth2
1Electrical Engineering & Applied Mathematics, 2Psychology & Neuroscience
University of Colorado at Boulder, Boulder, CO 80309
Abstract—We describe here the recent results of a multidis- neuronal firing. The development of electrophysiological
ciplinary effort to design a biomarker that can actively and biomarkers has focused on the analysis of both epileptiform
7 continuouslydecodetheprogressivechangesinneuronalorga-
spikes [5], and high frequency oscillations [6]. Interictal
1 nizationleadingtoepilepsy,aprocessknownasepileptogenesis.
spikes are sharp electrical impulses. Their morphology has
0 Using an animal model of acquired epilepsy, we chronically
2 record hippocampal evoked potentials elicited by an auditory been shown to be correlated with the progression of epilep-
stimulus.Usingasetofreducedcoordinates,ouralgorithmcan togenesis [5]. Recent experiments combined with computa-
n
identifyuniversalsmoothlow-dimensionalconfigurationsofthe tionalmodels[5]suggestthatepileptogenesissystematically
a
auditory evoked potentials that correspond to distinct stages
J modifies the morphology of the spikes in a predictable
of epileptogenesis. We use a hidden Markov model to learn
5 manner. High frequency oscillations detected in local field
the dynamics of the evoked potential, as it evolves along these
2 smoothlow-dimensionalsubsets.Weprovideexperimentalevi- potentials are created by synchronously bursting pyramidal
dencethatthebiomarkerisabletoexploitsubtlechangesinthe cells, and have been observed in healthy patients when the
]
C evokedpotentialtoreliablydecodethestageofepileptogenesis frequency of the ripples is lower than 250 Hz. Conversely,
and predict whether an animal will eventually recover from
frequencies in the range of 250Hz to 800Hz are considered
N
the injury, or develop spontaneous seizures.
to be pathological ripples and are reliable biomarkers for
.
o theepileptogeniczone[7].Whilerecordingsofspontaneous
i I. INTRODUCTION
b neuronal spiking can be indicative of neuronal excitability,
- A. The Challenge: Decoding Epileptogenesis andthereforecorrelatewiththepropensityforseizures,such
q
Epilepsyisaneurologicaldiseasethatischaracterizedby methodscannotactivelyprobethehippocampalcircuitinliv-
[
theoccurrenceofseveralunprovokedseizures.Despitevari- inganimalsduringepileptogenesis.Consequently,onecould

1 ouscausesofepilepsyandvaryingdegreesofdiseasesever- argue that the passive electrophysiological recordings may
v
ity in the human population, hippocampal sclerosis is the not provide enough information to observe early changes in
3
4 most consistent neuropathological feature of temporal lobe neuronal excitability associated with epileptogenesis.
2 epilepsy [1]. Animal models have been developed to study The present work addresses this limitation and proposes
7 the neuronal changes underlying the clinical manifestations for the first time a “computational biomarker” that relies
0 of epilepsy (chronic-spontaneous seizures). One popular on actively probing the excitability of the hippocampus
.
1 model relies on controlled administration of a convulsant using an auditory stimulus. We advocate an active approach
0 drug (e.g., pilocarpine) to induce status-epilepticus, a life- whereby we probe the excitability of the hippocampus, a
7 threatening condition in humans. This condition is followed brain region known to be epileptogenic. Because the hip-
1
by a latent seizure-free period of weeks to months, where pocampus also receives several sensory inputs (through the
:
v progressive neuronal damage and network reorganization entorhinalcortex),weproposetorecordinthehippocampus
Xi eventually lead to the development of spontaneous seizures. the evoked potential elicited by an auditory stimulus. Since
Most of our understanding of the progression of epilepsy, epilepsy does not modify the primary auditory cortex, any
r
a or epileptogenesis, is derived from such animal models. alterations in the evoked potential should be indicative of
It is therefore critical to define a biomarker to monitor neuronal changes in the hippocampus underlying epilepto-
epileptogenesis.Anaccuratebiomarkerwouldbeinvaluable genesis. To quantify the property of this biomarker, we use
for the design of novel anti-epileptogenic drugs, and could a pilocarpine animal model of temporal lobe epilepsy, and
eventually be translated into a diagnostic tool for humans. chronically record hippocampal auditory evoked potentials
during epileptogenesis.
B. From Passive Recording to Active Sensing
We design a decoding algorithm to demonstrate that
Current efforts toward the development of a reliable and changes in the morphology of the hippocampal auditory
predictive biomarker of epileptogenesis fall in three main evoked potential have universal predictive value and can be
classes: molecular and cellular biomarkers [2], imaging used to accurately quantify the progression of epilepsy. The
biomarkers[3],andelectrophysiologicalbiomarkers[4].The authorsarenotawareofanyworkthatusesmachinelearning
present work focuses on the last class of biomarkers that methods to construct an active biomarker of epileptogenesis
rely on recordings of the electrical activity associated with (but see [5] for approaches based on passive recordings).
condition baseline silent latent chronic
e
deuvroaktieodn h1( 0 )(wt)eek h( 1 )(t) 72 hho( k )u(tr)s weeks − months 2−6 months elin silent
inteprovteennttiioanl pilocarpine paraldehyde als bas es latent
Fevigenutre1. TimstaetluIisIn e.epialTenpdHticEnuos Am(SEeN)ncIMlatAurLeoMftOheDdEiLfferefinrstt scpoonndtaintieoonuss .seizure evoked potentirof animalc waveletcoefficients 1 delay coordinat baseline 2 IRdζ(rk) 3ζ(rk)
ni chronic
All procedures were performed in accordance with the o
hr spectral embedding
University of Colorado Institutional Animal Care and Use c h(k) w(k) z(k)
r r r
Committee guidelines for the humane use of laboratory rats
in biological research. Twenty-four male Sprague-Dawley Figure 2. Overview of the decoding algorithm. 1: a vector of wavelet
rats (200-250 gm) were implanted with a hippocampal wire coefficients,wr(k),iscomputedfromtheevokedp(cid:44)otentialh(rk).Avector
electrode, a ground screw, and a reference screw. The 24 of time-delay wavelets coordinates, zr(k), is formed by concatenating τ
rats were divided into two groups: a group of 17 rats that dcoisntsaenccuetibveetwwere(kn).ζr((cid:44)k2):asnpdectthrealloewm-bdeidmdeinnsgiomnaalpsstrwucr(tku)retofoζrmr(ke)d.b(cid:44)3y:eatchhe
received lithium-pilocarpine, and a control group of 7 rats. conditioniscomputed.
The control group was composed of 2 rats that received
all drug injections associated with the lithium-pilocarpine wavelet and scaling coefficients, wr(k), (see 1 in Fig. 2).
The second stage involves characterizing th(cid:44)e association
model except for pilocarpine, which was substituted with
between the condition of the disease (baseline, silent, la-
saline; and 5 rats that received no drugs. After full recovery
tent, or chronic) and the vector of wavelet coefficients
from the electrode implantation (2 weeks) and at least one
w(k). We tackle this question by lifting w(k) into Rτ×s
additionalweekofchronicrecordingofbaselinevideo/EEG, r r
using time-delay embedding: we concatenate the consecu-
17 rats were given an injection of lithium chloride followed
tive vectors w(k),...,w(k+τ−1) to form a τ × s vector,
by an injection of pilocarpine hydrochloride 24 hours later r r
(seeFig.1).Afteronehourofstatusepilepticus,theanimals zr(k), of time-delay wavelet coordinates (see 1 in Fig. 2).
were administered a dose of paraldehyde to terminate con- Low-dimensionalstructures,whichuniquelych(cid:44)aracterizethe
vulsions.Throughouttheexperiment,every∆t=30min,an stage of epileptogenesis, emerge in the high-dimensional
auditory stimulus, composed of a sequence of 120 square- space. We use spectral embedding to parameterize these
wave clicks (0.1 ms duration, 2 sec ISI, 45dB SPL), was low-dimensional structures, and map zr(k) to ζr(k) (see 2
played in a top-mounted speaker. The 300 ms hippocampal in Fig. 2). The first decoding stage involves geometrical(cid:44)ly
responsestoeachclickwerefilteredandsampledat10kHz, computing the likelihood that a given vector ζr(k) corre-
and the average of the 120 responses was computed. In spondstooneofthefourconditions.Tothisend,wequantify
the remainder of the paper, we denote by h(k) the average the distance of ζr(k) to the low-dimensional cluster formed
evokedpotential,measuredattimek∆t.Tofurthersimplify by each condition (see 3 in Fig. 2). In the final decoding
the exposition, h(k) is simply referred to as the evoked stage,weuseahiddenM(cid:44)arkovmodeltocapturetheintrinsic
potential measured at time k. dynamics of epileptogenesis.
Figure 1 provides a detailed timeline, along with the Toalleviatethenotation,andunlessweexplicitlycompare
nomenclature of the different periods associated with the or combine several animals, we dispense with the subscript
progress and eventual onset of epilepsy. The period before r, denoting the dependency on rat r, when we discuss the
the injection of pilocarpine is called baseline. Conversely, analysis of the evoked potentials for a fixed animal.
the period following the first spontaneous seizure is called
chronic.Wefurtherdefinethesilentperiodtobethe72hour IV. DENOISINGTHEINPUTDATA
periodofrecoveryimmediatelyfollowingtheterminationof To compensate for the alteration of the tissue impedance
status epilepticus, and the latent period to be the remaining around the electrodes, the evoked potentials h(k),k =
r
period leading to the eventual onset of the first spontaneous 0,1,...werenormalizedforeachanimal.Foragivenanimal
seizure. r, all evoked potentials for that animal were rescaled, in
such a way that the average energy computed during the
III. OVERVIEWOFTHEDECODINGFRAMEWORK baseline condition, (cid:104)h2(cid:105), was one. In order to use the noisy
r
We give here a brief overview of the decoding approach. evoked potentials to predict the state of epileptogenesis,
Given an animal r = 1,...,24, we consider the sequence we extract a denoised representation of h(k). We use a
of evoked potentials h(r0),h(r1),.... The first stage involves discrete stationary wavelet transform (CDF 9-7) to compute
theconstructionofadenoisedrepresentationofh(rk) (see 1 a redundant representation of h(k). We used a “leave-one-
inFig.2).Wedecomposeh(k)(t)usingadiscretestationa(cid:44)ry animal-out” cross-validation to determine the time intervals
r
wavelettransformandretainavectorofs(carefullychosen) andthescalesofthewaveletcoefficientsthatresultedinthe
baseline silent latent chronic
70 45
930 930 0.6 930 0.6 930 0.6 0.6 60 baseline 3450 silent
Frequency (Hz)4216315217948894 Frequency (Hz)4215217631894 948 ---000000..24...642 Frequency (Hz)4215217631894 948 ---000000..24...642 Frequency (Hz)4215217631894 948 ---000000..24...642 ---000000..24...642 23450000 mσ e=a 0n. 0=7 4.0 1122305050 mσ e=a 0n. 2=0 3.53
Frequency (Hz)942152136318940948 50 100Tim1e5 (0ms2)00 25Frequency (Hz)094215213631389409480 0 -000050...2460.2 100Tim1e5 (0ms2)00 25Frequency (Hz)09421521363189430948 0 0 -000050...2460.2 100Tim1e5 (0ms2)00 25Frequency (Hz)094215213631389409480 0 -000050...246.02 100Tim1e5 (0ms2)00 250 300 -00000...246.2 1111460803000.5 lat3e.6nt 3.7 D3im.8ensio3n.9 4 4.1 4.2 110520300.2 c3h.3roni3c.4 3.5Dime3.n6sion3.7 3.8 3.9 4
7 50 1T00im1e5 0(m2s00)2507 50 10T0im15e0 (2m00s)2507 50 10T0im15e0 (2m00s)2507 50 100Tim15e0 (2m00s)250 300 110200 mσ e=a 0n. 1=9 3.68 80 mσ e=a 0n. 2=1 3.69
60
80
Figure 3. Top: average wavelet coefficients with the average evoked 60 40
potential(acrossallanimals).Bottom:averageapproximationcoefficients. 40
Onlythescalesj=3(top)to10(bottom)aredisplayed.Whiterectangles 20 20
delineatethetime-frequencyblocksusedtoconstructw(k). 03 3.2 3.4 3.6 3.8 4 4.2 03.13.2 3.3 3.4 3.D5ime3n.6sion3.7 3.8 3.9 4 4.1
the most accurate prediction of the condition of the animal Figure4. Distributionofthelocaldimensionofthestatespaceaccording
tothecondition.Top:baselineandsilent;bottom:latentandchronic.
that was not part of the training data. For most of the ani-
mals, the wavelet coefficients in the time-frequency region
embedding, yielded poor prediction of the animal’s condi-
[0,100] ms×[5,10] Hz, and the approximation coefficients
tion. While τ was optimized for each animal, the average
from the time-frequency regions [70,120] ms×[10,20] Hz
of the optimal values was τ =12, which corresponds to six
and [100,150] ms×[5,10] Hz (see Fig. 3) resulted in either
hours.
the optimal, or very close to optimal, decoding perfor-
To learn the geometry of the set formed by the different
mances. To simplify the decoding, we therefore kept the
trajectories z(k),k =0,1..., we consider the union – over
same time-frequency regions for all animals. In summary, r
each evoked potential was represented with a feature vector all the animals in the training set – of the vectors zr(k), and
w(k)of2,000entriescomposedof1,000waveletcoefficients define the set
and 1,000 approximation coefficients. This representation is Z= (cid:91) (cid:110)z(k),k =0,1,...(cid:111). (2)
consistent with reports of disruption in the θ rhythm (4-12 r
r ∈training animals
Hz) during the latent period preceding the onset of epilepsy
[8]. B. What is the local dimension of the state-space?
Because we have only a limited number of training data,
V. UNIVERSALCONFIGURATIONSOFEPILEPTOGENESIS
weneedtodrasticallyreducethedimensionalityofthestate
We now describe the construction of universal (stable space (on average τ × s = 24,000) in order to reliably
across all animals) low-dimensional smooth sets that are decodetheconditionofeachanimal.Weusedakernel-based
formed by all the hippocampal auditory evoked potentials versionofthecorrelationdimension[9]tocomputethelocal
collected during the same stage of epileptogenesis. These dimension of the point-cloud formed by the state space Z.
sets are created by lifting the wavelet coefficients w(k) Wediscoveredthatthepointsz(k) wereorganizedalonglow
r
of each h(k) into high-dimension. This lifting effectively
dimensional subsets that correspond to well defined stages
createssmoothlow-dimensionalcoherentstructuresthatcan
of epileptogenesis and exhibit different local dimensions.
then be parameterized with a drastically smaller number of
Figure 4 displays histograms (as well as the mean, and
coordinates.
standard deviation) of the estimates of the local dimensions
for the four different conditions. These experiments suggest
A. Time-Delay Embedding of the Wavelet Coordinates
thatweshouldbeabletorepresentthestate-spaceassociated
Given the time series (cid:8)w(k),k =0,1,...(cid:9) for a given
with the dynamic of the disease using d = 5 dimensions.
animal, we analyze the dynamics of this time series by
However,ourresults(notshown)indicatethatthetraditional
considering the time-delay wavelet coordinates formed by
linearapproach,PCA,providesaverypoorparameterization
concatenating τ consecutive vectors of wavelet coefficients,
of the set Z.
z(k) =(cid:2)w(k) w(k+1) ··· w(k+τ−1)(cid:3). (1)
C. Nonlinear Parameterization of the State Space
We characterize the dynamical changes in the evoked po- In order to identify the separate regions of Rτ×s that
tentials by studying the geometric structures formed by the correspond to the four conditions, we seek a smooth low-
trajectory of z(k) in Rτ×s, as k evolves. In practice, the dimensional parameterization of Z. A nonlinear approach
number of time-delay vectors, τ, is determined using cross- spectral embedding [10], yields a low-dimensional parame-
validation. We confirmed that τ =1, to wit no time-delay terizationofZthatnaturallyclustersthedifferentconditions
chronic are shown along this trajectory, confirming changes in the
morphology of h(k) during epileptogenesis.
baseline
0.08 D. Learning Epileptogenesis from the Reduced Coordinates
0.06
ThenaturaldivisionofZintocoherentsubsets,whichcor-
0.04 100μV
0.02 respond to well defined stages of epileptogenesis, suggests
0 100 ms that a purely geometric algorithm could be used to quantify
−0.02 −0.0−20.03 the development of epilepsy. Indeed, given an unclassified
latent −−00.0.064 0.010 −0.01 evoked potential, ζ(k), the distance from ζ(k) to each of the
−0.05−0.04−0.03−0.02−0.01 0 0.010.020.030.040.050.060.050.040.030.02 silent faonuersctilmusatteersof(btahseelliinkee,lishioleondt,olfatbeenitn,ganindtchherocnoircr)espproonvdidinegs
condition at time k.
Figure 5. The training set of evoked potentials, Z, displayed using the
reduced coordinates ζ(k). Each condition, baseline (blue), silent (cyan), We denote by Zc, c = 0,1,2,3, the four clusters formed
latent(green),andchronic(red),formsacoherentsub-cloud. by the evoked potentials in Z that were collected during the
baseline, silent, latent, and chronic conditions, respectively.
into coherent disconnected smooth subsets. Briefly, we de-
Wefoundthatusingamixtureofprobabilisticprincipalcom-
fineasimilaritymatrixKthatquantifieshowanytwoevoked ponents [11] to represent each Zc resulted in a remarkably
potentials h(k) and h(l) extracted from the same or from
low dimensional parameterization of that condition. Indeed,
different conditions, and from the same or from a different
the mixture model is able to describe with a small number
animal, at the respective times k and l co-vary, of components the geometric structure created by each Zc.
K(k,l)=exp(cid:16)−(cid:107)z(k)−z(l)(cid:107)2/σ2(cid:17). (3) Furthermore, unlike other models, the mixture model does
not require the stringent assumption that Zc be a smooth
The scaling constant σ was chosen to be a multiple of the sub-manifold.
median distance (cid:107)z(k)−z(l)(cid:107), k,l =0,1,.... We used the Formally, for each c = 0,1,2,3, we use a different mix-
d eigenvectors of K that optimally separated the dataset ture of probabilistic principal components to parameterize
Z into four clusters related to the four conditions. The Zc. To reduce the number of indices, we do not make
ith eigenvector of K provided the ith reduced coordinate, the dependency on c explicit in the following discussion.
ζ(k)(i), of z(k), For a given cluster Zc, the mixture is composed of M
local Gaussian densities, each describing the local principal
ζ(k) =(cid:2)ζ(k)(1) ··· ζ(k)(d)(cid:3)T . (4)
directions with a d×q matrix, U , around a collection of
i
AsexplainedinsectionV-B,weused=5.Figure5displays centers ζi, i = 1,...,M. The local spread of the ζr(k)
the nonlinear parameterization of Z. Each dot represents around ζi is given by the covariance matrix UiUiT +ε2iI.
an evoked potential h(k) measured at a time k from an The mixture model associated with the cluster Zc can be
r
animal r, and parameterized by the reduced coordinates written as
ζ(k). The color indicates the condition during which h(k) M
r r ζ(k) ∼(cid:88)π N (cid:0)U n +ζ , ε2I(cid:1), (5)
was recorded. When collapsed across animals and time of r i i i i i
recordings, four distinct clusters can be visually discerned, i=1
which can be interpreted in terms of the corresponding where the ni,i=1,...,M are q-dimensional latent Gaus-
conditions. The evoked potentials in the silent condition are sian variables, and the mixture proportions πi are non
grouped together around a low-dimensional structure. The negative and sum to one. The parameters of the model are
latentconditiondisplaysthelargestspread:someevokedpo- estimated using the EM algorithm [11]. For a given animal
tentialsaremorphologicallyclosetosilentevokedpotentials; r, we compute the posterior probability that the model (5)
while others are close to chronic evoked potentials. These associated with condition c would generate the reduced
visual observations are consistent with the computation of coordinates ζr(k),
thelocaldimensionalitydescribedinsectionV-B.Tofurther γ(k)(c)=Prob(c|ζ(k)). (6)
help with the interpretation of this nonlinear representation r r
of Z, we trace the trajectory of ζr(k), k = 0,1,... for a Weusetheposteriorprobabilityofclassmembershipγr(k)(c)
single rat (H37) over several weeks of recordings (see thick asarawmeasurementoftheprogressionofepileptogenesis.
black line in Fig. 5). The animal first spends several days in Because the instantaneous posterior probability γr(k)(c) is
the baseline cluster (blue), then briskly traverses the entire oblivious to the past trajectory {ζ(k−1),ζ(k−2),...}, it can
r r
spacetoreachthesilentcluster(cyan)afterstatusepilepticus benoisy.Toaddressthislimitation,andenforcethetemporal
has been induced. Eventually, as the animal recovers, it coherence that defines the progression of the disease, we
joins the latent condition (green), and advances toward the introduce a Markovian description of the dynamics of the
chronic cluster (red). Four representative evoked potentials changes in the evoked potential h(k) as a function of time.
SE baseline
1 silent
latent
0.5 chronic
0
0 500 1000 1500 2000 2500 3000 3500 4000
Time (0.5 hr)
SE SRS baseline
1
silent
0.5 latent
chronic
0
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
Time (0.5 hr)
Figure6. Biomarkerprobabilitiesπ(k)(c)ofbeinginconditioncforanon-epileptic(top)andanepileptic(bottom)animal.Seizurecount(clippedat1)
isdisplayedinmagenta.TheblackbarSEcorrespondstoStatusEpilepticus;thesecondbarSRScorrespondstothefirstSpontaneousRecurringSeizure.
VI. THEDYNAMICSOFEPILEPTOGENESIS seizure (SRS), and shifts from latent (green) to chronic
(red) before SRS. The probability of being in the chronic
We use a hidden Markov Model to capture both the
state, π(k)(3), remains at one after SRS, indicating that the
dynamics of epileptogenesis, and the resulting changes in
the shape of h(k) triggered by this hidden condition. We animal is irreversibly in the chronic state, as confirmed by
r the uninterrupted sequence of seizures.
define x(k) to be a discrete random variable that encodes
r
the state of animal r at time k, x(rk) ∈ {0,1,2,3}, where VII. RESULTS
the states 0, 1, 2, and 3 encode the baseline, silent, latent,
We evaluated our approach using a “leave-one-animal-
andchronicconditions,respectively.Wefurtherassumethat
out” cross-validation: the algorithm was trained on 11
x(k) is a Markov chain with probability transition matrix
r epileptic rats, and evaluated with one animal that was not
P = Prob(x(k+1) = i|x(k) = j). We define y(k) to
i,j r r r partofthetrainingdata.Foreachtestanimalr,thedecoding
be a discrete random variable, which takes its values in algorithm computed from h(k) the posterior probability
r
{0,1,2,3},andencodesthemostlikelyconditionofanimal π(k)(c), given by (8). Figure 7 displays the temporal profile
r at time k, given the measurement h(k), r
r of the median probability π(k)(c) computed over all the
y(k) =argmax γ(k)(c), (7) epilepticanimals(left,N=12),andthenonepilepticanimals
r r
c=0,1,2,3 (right, N=5). The computation of the median probability
π(k)(c) required care: because each condition varied in
where γ(k)(c) is the posterior probability that condition c
r length across the different animals, we could not use the
generates the reduced coordinates ζr(k), and is given by beginning of each recording as the origin, and bluntly
(6). We denote by Qi,j = Prob(yr(k) = i|x(rk) = j) the average across the animals.
measurementprobabilitydistributionmatrix.Theprobability As an alternative, we used two physiologically mean-
transition matrix P and the measurement probability matrix ingful temporal landmarks to align the time series
(cid:110) (cid:111)
Qwereestimatedfromthetrainingdata.Thetrainedhidden π(k)(c),k =0,1,... of the different animals before av-
r
Markov model was then used in a “reverse mode” [12]
eraging. We used the time of the insult (status epilepticus)
to estimate the posterior probability of animal r to be in
to study the transition from baseline to silent, and we also
condition c,
used the onset of the first spontaneous seizure to study the
π(k)(c)=Prob(x(k) =c|y(k)), c=0,1,2,3. (8) transitionfromlatenttochronic.AsshowninFigs.7-leftand
r r r
Fig. 8, this led to two different median probability profiles
Because the hidden Markov model was trained on the for each condition. In each figure, the origin corresponds to
epileptic animals, we regularized the probability transition the physiological landmark used to realign the time series.
matrixP toallowatransition,whichwasneverobservedin Tohelpvisuallycomparethetwofigures,weaddedthelatent
thetrainingdata,fromchronictobaselinewithaverysmall periodtoFig.7-left,andthesilentperiodtoFig.8.Wenote
probability. that the median probability for the silent and latent periods
Figure 6 displays the four probabilities π(k) of being in are different in both figures: indeed, the large variation in
baseline, silent, latent, and chronic state for a non-epileptic the onset of the first spontaneous seizure leads to different
animal (top), and an epileptic animal (bottom). We plot the realignment results.
seizure count (clipped at 1) in magenta for the epileptic Epileptic animals. The baseline and chronic periods were
animal. The biomarker clearly identifies the baseline pe- welldefinedforthe12epilepticanimals(seeFig.7-leftand
riod (blue), and is able to detect changes in the evoked 8), and corresponded to the period before status epilepticus,
potential h(k) that are indicative of neuronal alterations and the period after the first spontaneous seizure, respec-
leading to epilepsy, before the first spontaneous recurring tively.Thesensitivitywasfoundtobe1forboththebaseline
1 baseline 1
silent baseline
0.5 latent 0.5 silent
latent
chronic
0 0
-10 0 10 20 30 -10 0 10 20 30
days days
Figure 7. Median probability π(k)(c) of being in condition c for the epileptic animals (left), and the non epileptic animals (right). Day 0 (black bar)
correspondstostatusepilepticusacrossallanimals.
1
silent
latent
0.5 chronic
seizures
0
-10 0 10 20 30 40 50 60 70 80
days
Figure 8. Median probability π(k)(c) of being in condition c for the epileptic animals. Day 0 (black bar) corresponds to the first spontaneous seizure
acrossallanimals.Seizurecount(clippedat1)isdisplayedinlightgrey.
and chronic periods. The specificity was 1 for the baseline [2] K. Lukasiuk and A. J. Becker, “Molecular biomarkers of
period, and 0.59 for the chronic period, indicating that the epileptogenesis,” Neurotherapeutics, vol. 11, no. 2, pp. 319–
probability π(k)(3) rises to 1 before the first spontaneous 323, 2014.
seizure (see Fig. 8). Indeed, the biomarker is able to detect
[3] S. R. Shultz, T. J. O’Brien, M. Stefanidou, and R. I.
changes in the evoked potential h(k) that are indicative Kuzniecky, “Neuroimaging the epileptogenic process,” Neu-
of neuronal alterations leading to epilepsy, before the first rotherapeutics, vol. 11, no. 2, pp. 347–357, 2014.
spontaneous seizure.
[4] R. Staba, M. Stead, and G. Worrell, “Electrophysiological
Non-epileptic animals. Five of the 17 animals, which
biomarkers of epilepsy,” Neurotherapeutics, vol. 11, no. 2,
receivedtheconvulsant(pilocarpine),neverdevelopedspon- pp. 334–346, 2014.
taneousseizures.Aftertheinjectionsofparaldehyde,andthe
subsequent recovery from status epilepticus, these animals [5] C.Huneau,P.Benquet,G.Dieuset,A.Biraben,B.Martin,and
F.Wendling,“Shapefeaturesofepilepticspikesareamarker
returned to a baseline condition. These animals were inter-
of epileptogenesis in mice,” Epilepsia, vol. 54, no. 12, pp.
esting,astheyallowedustoaddressacommonlimitationof
2219–2227, 2013.
machinelearningbasedbiomarkerswithregardtospecificity
of the biomarker to the general population. The biomarker [6] J.EngelandF.daSilva,“High-frequencyoscillations–where
we are and where we need to go,” Prog. Neurobiol., vol. 98,
was able to predict the return to baseline condition (see
no. 3, pp. 316–318, 2012.
Fig. 7-right) in spite of the fact that the training dataset
never included non-epileptic rats. This result is significant [7] L.M.delaPrida,R.J.Staba,andJ.A.Dian,“Conundrums
since it indicates that the low-dimensional clusters formed of high-frequency oscillations (80–800 hz) in the epileptic
by embedding the time-delay wavelet coordinates have a brain,” Journal of Clinical Neurophysiology, vol. 32, no. 3,
pp. 207–219, 2015.
universalaspectthatcanbereliablyusedtodecodethestatus
of the disease. [8] L. Chauvie`re, N. Rafrafi, C. Thinus-Blanc, F. Bartolomei,
Control animals. The 5 rats that received no drugs where M.Esclapez,andC.Bernard,“Earlydeficitsinspatialmem-
classifiedasbeinginthebaselineconditionwithprobability oryandthetarhythminexperimentaltemporallobeepilepsy,”
J. Neurosci., vol. 29, no. 17, pp. 5402–5410, 2009.
1 at all time (results not shown). The 2 rats that received all
drug injections except for pilocarpine, were also classified [9] M. Hein and J.-Y. Audibert, “Intrinsic dimensionality esti-
as being in the baseline conditions at all time, except for a mation of submanifolds in Rd,” in Proceedings of the 22nd
brief suppression of amplitude, classified as silent (results internationalconferenceonMachinelearning. ACM,2005,
pp. 289–296.
not shown).
[10] M. Ghil, M. Allen, M. Dettinger, K. Ide, D. Kondrashov,
ACKNOWLEDGEMENT M. Mann, A. W. Robertson, A. Saunders, Y. Tian, F. Varadi
et al., Spectral dimensionality reduction. Springer, 2006.
FGM was supported in part by NSF DMS 1407340.
[11] M.TippingandC.Bishop,“Mixturesofprobabilisticprinci-
REFERENCES palcomponentanalyzers,”Neuralcomputation,vol.11,no.2,
pp. 443–482, 1999.
[1] A. Pitka¨nen and K. Lukasiuk, “Mechanisms of epileptogen-
[12] O. Cappe´, E. Moulines, and T. Ryde´n, Inference in Hidden
esis and potential treatment targets,” The Lancet Neurology,
Markov Models. Springer, 2009.
vol. 10, no. 2, pp. 173–186, 2011.