State and parameter estimation of a neural mass model from electrophysiological signals during the status epilepticus.

YNIMG-12030; No. of pages: 13; 4C: 3, 4, 5, 6, 7, 8, 9, 10, 11 NeuroImage xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

State and parameter estimation of a neural mass model from electrophysiological signals during the status epilepticus Armando López-Cuevas a,⁎, Bernardino Castillo-Toledo a, Laura Medina-Ceja b, Consuelo Ventura-Mejía b a b

Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, CINVESTAV, Unidad Guadalajara, Av. del Bosque 1145, Col. El Bajío, Zapopan, 45015 Jalisco, Mexico Laboratorio de Neurofisiología y Neuroquímica, Departamento de Biología Celular y Molecular, CUCBA, Universidad de Guadalajara, Mexico

a r t i c l e

i n f o

Article history: Accepted 26 February 2015 Available online xxxx Keywords: Electrical stimulation device Closed loop stimulation Cubature Kalman filter Epilepsy Neural mass model Pilocarpine Population model Real time parameter estimation Status epilepticus

a b s t r a c t Status epilepticus is an emergency condition in patients with prolonged seizure or recurrent seizures without full recovery between them. The pathophysiological mechanisms of status epilepticus are not well established. With this argument, we use a computational modeling approach combined with in vivo electrophysiological data obtained from an experimental model of status epilepticus to infer about changes that may lead to a seizure. Special emphasis is done to analyze parameter changes during or after pilocarpine administration. A cubature Kalman filter is utilized to estimate parameters and states of the model in real time from the observed electrophysiological signals. It was observed that during basal activity (before pilocarpine administration) the parameters presented a standard deviation below 30% of the mean value, while during SE activity, the parameters presented variations larger than 200% of the mean value with respect to basal state. The ratio of excitation– inhibition, increased during SE activity by 80% with respect to the transition state, and reaches the lowest value during cessation. In addition, a progression between low and fast inhibitions before or during this condition was found. This method can be implemented in real time, which is particularly important for the design of stimulation devices that attempt to stop seizures. These changes in the parameters analyzed during seizure activity can lead to better understanding of the mechanisms of epilepsy and to improve its treatments. © 2015 Elsevier Inc. All rights reserved.

Introduction Epilepsy is a brain disorder characterized by recurrent seizures (Fisher et al., 2005) affecting 2–5% of the world's population (Angus and Parsons, 2008). Epilepsy is generated by abnormal, hypersynchronic neuronal activity in the brain; its onset can involve several regions (generalized seizures) or just a circumscribed brain region (focal seizures) (Engel, 2006). The causes of epilepsy are multifactorial, including injuries, infections, abnormal brain development, unbalance in neurotransmitters, and brain tumors among others. Epilepsy treatments include pharmacology or surgical methods but even with advances in medicine approximately 30% of patients remain with seizures (Perucca, 2007). Some seizures can persist for long time, causing severe brain damage or even death. This state of long-lasting seizures is known as status epilepticus (SE). According to the International League Against Epilepsy (Commission on Classification and Terminology of the International League Against Epilepsy, 1981) SE is defined as a prolonged seizure or recurrent seizures without full recovery between them and in the literature it is assumed that SE is a seizure lasting longer than 30 min (Brodie, 1990). The long-lasting SE effects include severe ⁎ Corresponding author at: Center for Research and Advanced Studies, CINVESTAV, Mexico. Fax: +52 3337773609. E-mail address: [email protected] (A. López-Cuevas).

clinical conditions resulting in disorders of various organs and systems as well as neurological damage that could lead to death (Jenssen et al., 2006). Termination of an SE episode is related to neuroprotection and better prognosis in patients. Therefore it is important to search in detail the neuronal mechanisms that lead to the SE episodes in order to carry out better diagnosis and treatment of SE given the poor knowledge about the pathophysiological mechanisms of SE that lead to the loss of the ability to spontaneous seizure self-termination. The possibility of estimating the internal changes that lead to seizures during SE would make it possible to design devices able to deliver drugs or electrical stimulation automatically in real time to reduce or eventually eliminate such pathological activity. Moreover, if a model that describes the electrical activity in a certain region of the brain can be obtained, then a feedback control strategy could be designed to maintain the brain activity within normal parameters. Computational modeling in neuroscience is a well-established and increasing research area that has contributed to better understanding of the mechanisms underlying certain brain activities or phenomena. When combined with experimental studies it can be a powerful tool in brain research, especially in brain disorders (Goodfellow et al., 2012; Lytton, 2008; Traub et al., 2010). In this work, we use neural mass modeling (Jansen and Rit, 1995) combined with in vivo experiments to investigate whether a change in the states and parameters of the model can be observed during the transition to SE, where SE is a

http://dx.doi.org/10.1016/j.neuroimage.2015.02.059 1053-8119/© 2015 Elsevier Inc. All rights reserved.

Please cite this article as: López-Cuevas, A., et al., State and parameter estimation of a neural mass model from electrophysiological signals during the status epilepticus, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.02.059

2

A. López-Cuevas et al. / NeuroImage xxx (2015) xxx–xxx

self-sustained seizure activity. By using a cubature Kalman filter (Arasaratnam and Haykin, 2009) we are able to estimate on-line both the state and parameter changes. The paper is organized as follows: In section 2, the experimental methods carried out for data acquisition and the animal model of SE used in this study are described, the obtained signals are analyzed in frequency and time domains and differences between basal and ictal signals are put in evidence; a detailed description of the neural mass model used here is given, and a comparison between signals produced by the model and actual signals is done. In section 3 a model-based state and parameter identification from electrophysiological signals is presented, discussion and conclusions are drawn in sections 4 and 5 respectively.

Histology

Materials and methods

The intracranial EEG recordings showed high-voltage spikes during the first few minutes after pilocarpine injection (right and left DG and CA1 areas). The onset of SE in both animals was observed 5 ± 3 min after pilocarpine injection in the right hemisphere and propagated to the left hemisphere 1 h later and it was characterized by high-voltage spikes. The complete recording of a representative experiment is shown in Fig. 1. There are four signals, from top to bottom: right DG, right CA1, left DG and left CA1. The figure shows four vertical lines, from left to right, the first (red) indicates the time when pilocarpine was administrated, the second line (black dashed) indicates the time of status epilepticus onset in the right hemisphere, the third line (black dashed) shows the time of status epilepticus onset in the left hemisphere, the fourth vertical line shows the moment of SE cessation. There are two important things to observe in Fig. 1; the first one is that SE begins in the right hemisphere in both DG and CA1 area, a few minutes after pilocarpine administration, nevertheless, in the left hemisphere SE starts almost an hour later. This is because injection was administrated on the right lateral ventricle (see Experimental setting). The second issue is that SE ceases in the four channels approximately in the time 10,000 s. A time–frequency analysis of the signal was carried out, empirical mode decomposition (Huang et al., 1998) was applied to the signals in order to separate frequency bands and then, a Hilbert–Huang transformation was carried out to obtain the time–frequency spectrum of the signals. Additionally, Fourier transform was applied to selected segments of the signal to obtain their main frequency components in order to get information for the design of the model. The Hilbert spectrum of the complete recording in the right and left DG and CA1 areas is shown in Fig. 2. As can be noted the power increases during seizure activity, in the Theta band (4–12 Hz) but also in higher frequencies. The main frequency components in three different segments of the recording, corresponding to basal activity, and SE activity are shown in Fig. 3, the segments shown in this figure correspond to the selected segments (marked with asterisks) in Fig. 1.

Pilocarpine model of epilepsy Seizures are generated by abnormal hypersynchronic neuronal activity in the brain; its onset can involve several regions or just a circumscribed brain region. Pilocarpine is a muscarinic agonist used to induce seizures to simulate several characteristics of SE. This experimental model was chosen because the electrophysiological activity from pilocarpine treated animals resembles to that in SE condition in humans (Medina-Ceja et al., 2014; Cavalheiro et al., 1996). To test and support our model with real data, the following procedures were carried out. Male Wistar rats (190–200 g) were maintained individually in a temperature-controlled room (22 ± 2°C) on a 12 h light/dark cycle (lights on at 7:00 a.m.), with ad libitum access to food and water. All experimental procedures were designed to minimize animal suffering, and the experimental protocol was in accordance with the Rules for Research in Health Matters (Mexican Official Norms NOM-062-ZOO-1999, NOM-033-ZOO-1995) and it was approved by the local Animal Care Committee. The electrophysiological data obtained in two animals with SE were used in the present study. These animals were anesthetized with isofluorane (Sofloran, PISA, Laboratories, Mexico) in 100% oxygen and then secured in a Stoelting stereotaxic frame such that lambda and bregma were in the same horizontal plane. Fixed recording microelectrodes, four 60 μm diameter pairs of tungsten wires with a vertical tip separation of 1.5 mm, were implanted bilaterally at symmetrical locations in the posterior hippocampus (CA1, anteroposterior, AP, − 5.0 mm relative to bregma; mediolateral, ML, ± 5.0 mm; deepvertical, DV, −5.5 mm) and anterior hippocampus (dentate gyrus DG, AP, − 3.5 mm; ML ± 2.00 mm; DV − 4.0 mm). The microelectrodes were attached to a pin connector and fastened to the skull with dental cement. After a week of recovery, the rats were allowed to move freely and their behavior was recorded. Five 4-channel MOSFET small amplifiers were placed on the cable connector to eliminate movement artifacts, and the electrical activity in the hippocampus was recorded using a 7D polygraph with eight amplifiers (Grass Technologies, RI, USA) and a wide-band (0.1–3 kHz). The sensitivity was 75 μV/cm per channel and a 5 kHz/channel sampling rate was used with 12 bit precision. Experiments were performed using an iMac A1048 (Apple, USA) and MP150 software (BIOPAC Systems, CA, USA). Matlab™ routines were used to analyze the signals (López-Cuevas et al., 2013). Each electrode signal accounts for the electrical activity of a population of neurons in a vicinity of the implanted microelectrode. A hole was drilled in the rats' skull, above the right lateral brain ventricle at the following stereotaxic coordinates relative to bregma: AP −4.1 mm, L −5.2 mm, V 7 mm. A single dose of pilocarpine hydrochloride (2.4 mg in a total volume of 2 μl; Sigma-Aldrich, USA) was injected using an injection pump that was attached to the stereotaxic frame (flow 1 μl/min; Stoelting Co. IL. USA), after which the animals were scored according to the Racine scale (Racine, 1972). Animals with a score of 4/5 were considered to have developed SE.

After each experiment, the animals were anesthetized with sodium pentobarbital, perfused transcardially with 100 ml of normal saline (0.9)% in 0.12 M buffer/CaCl2, followed by 300 ml of 4% paraformaldehyde in 0.12 M buffer/CaCl2 at pH 7.3. The brain of the animals was then removed, and 50 μm thick coronal sections were cut and stained with cresyl violet to confirm the correct positioning of the microelectrodes. Also, animals were injected with direct blue to directly corroborate the position of the needle in the right lateral ventricle. Signal analysis

Neural mass model The present work proposes a model based on real data obtained from intracranial EEG recordings of experimental rats by using the neural mass model methodology (Molaee-Ardekani et al., 2010; Zavaglia et al., 2006; David and Friston, 2003). This methodology (also known as population model) allows for the representation of local field potentials (LFP) as the result of the interaction among populations of neurons, where some populations are formed by purely excitatory neurons and other populations are formed by purely inhibitory neurons. Then, the neural mass models and the cubature Kalman filter (CKF) (Arasaratnam and Haykin, 2009) are combined to estimate in real time the states and parameters of the model. It is assumed that certain parameters of the model are time-varying and that the variation of these parameters eventually leads to epileptiform discharges. Some recent studies, which have been carried out from a theoretical



Pilocarpine

SE onset

5

−5 0

*

* 1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

Right

0 −5 0

R−SE onset 10

*

0 −10 0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

5

Left hemisphere

Amplitude (mV)

5

hemisphere

0

3

0 −5 0

Time (s) Fig. 1. Complete records of the experiment and the: from top to bottom, right dentate gyrus, right CA1, left dentate gyrus and left CA1. The vertical lines in the figure represent the time of pilocarpine injection (red line), status epilepticus onset (black dashed lines) and cessation of status epilepticus (dashed red line). The asterisks show three representative segments that are shown in Fig. 3.

approach in the frame of bifurcation of dynamical systems, have analyzed how certain parameter variations affect the model behavior (Spiegler et al., 2010). In other studies, more detailed models have been used and compared to real data with good results, nevertheless, estimation of parameters using genetic algorithms is done off-line

(Wendling et al., 2005). In (Friston et al., 2003) a Bayesian approach and dynamic causal modeling are used to estimate the effective connectivity among populations using functional magnetic resonance imaging (fMRI) time series, in (Havlicek et al., 2011) the CKF is utilized to estimate states of an fMRI model.

Right Dentate Gyrus

50 40 30 20 10 0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

6000

7000

8000

9000

10000

6000

7000

8000

9000

10000

6000

7000

8000

9000

10000

Frequency (Hz)

Right CA1

50 40 30 20 10 0

1000

2000

3000

4000

5000 Left Dentate Gyrus

50 40 30 20 10 0

1000

2000

3000

4000

5000 Left CA1

50 40 30 20 10 0

1000

2000

3000

4000

5000

Time (s) Fig. 2. Normalized Hilbert spectrum of the signals. Red color indicates the maximum energy and blue color the lowest energy.


4


0

1

2

0 0

3

1

0 −2 0

1

2

20

30

40

50

60

10

20

30

40

50

60

10

20

10

15

20

60

d)

0.5

0 0

3

1

e)

2 0 −2 0

10

2

Power ( V )

c)

2

4

0.5

f)

2

Amplitude (mV)

4

b)

2

Power ( V )

2

−2 0

Amplitude (mV)

1

a)

Power ( V )

Amplitude (mV)

4

1

2

0.5

0

3

Time (s)

0

Frequency (Hz)

Fig. 3. Three representative types of activity (left) and their respective frequency spectrum (right). The segments in this figure correspond to the segments marked with asterisks in Fig. 1: a, b activity from basal state in right DG and its frequency spectrum, respectively (red asterisk in Fig. 1). c, d activity from SE in left dentate gyrus and its frequency spectrum respectively (black asterisk in Fig. 1). e, f activity from SE in right dentate gyrus and its respective frequency spectrum (green asterisk).

Neural mass models represent field potentials in a local area of the brain, in our case in an area adjacent to the implanted microelectrodes. It is assumed that the resulting signal is the result of interacting neural populations, specifically excitatory and inhibitory populations (Jansen and Rit, 1995). The model can be visualized as formed by three principal components, the first of them the external input which represents activity coming from other areas (populations) of the brain. This input is usually represented by random noise. The second component is the populations of neurons presented in that particular area whose average membrane potentials are modeled by differential equations as a function of the population firing rates. The third component is the output of the model which is a combination of the states of the model. In Fig. 4 a generic scheme of a neural mass

model is shown, where p(t) is the external input, the main excitatory population is delimited by a dashed oval, the inhibitory population is depicted at the bottom delimited by a dashed rectangle, and a secondary excitatory population appears at the top of the figure. Each rectangular block labeled h(t) accounts for an impulse response block whose inputs are firing rates and its output is a postsynaptic potential (PSP), similar to dendrites that receive excitatory or inhibitory synapses and transform them in postsynaptic membrane potentials. These blocks carry out a convolution between the incoming action potential sequence, i.e., firing density (Lopes da Silva et al., 2003) and synaptic impulse response function, the result are PSP's that are added linearly. The blocks representing a sigmoid, convert the postsynaptic potential values to a pulse density i.e., the firing rate of the population, through a nonlinear sigmoid function (Freeman, 1987; Wilson and Cowan, 1972). The synaptic impulse response can be represented by a second order differential equation as

y¼ z 2 z¼ Aaxðt Þ þ 2azðt Þ þ a yðt Þ

ð1Þ

where x(t) is the input, y(t) is the output and parameters A and a shape the synaptic response. These parameters have different values for each neurotransmitter, i.e., GABA, Glutamate. The sigmoidal function that transforms postsynaptic potential in firing rate is described by the equation

sigmðvÞ ¼ Fig. 4. Schematic representation of a general population model. Rectangles with the letter hx represent the conversion from pulse density to postsynaptic potentials and rectangles with a sigmoid figure inside represent the nonlinear conversion from postsynaptic potentials to pulse density. Abbreviations: FR, firing rate; PSP, postsynaptic potential. C1, C2, C3 and C4 are parameters of the model.

2e0 1 þ erðv−v0 Þ

ð2Þ

where, e0 represents the maximum firing rate, v0 the PSP for which a 50% firing rate is achieved, and r the steepness of the sigmoidal transformation.



Dentate gyrus model

Table 1 Parameters of the DG model.

The DG is a cortical region that is part of a larger functional brain system called the hippocampal formation. It lies between the entorhinal cortex and the cornu ammonis area and it is believed to play an important role in preprocessing information coming from cortical areas that ultimately leads to the production of episodic memories (Amaral et al., 2007). The dentate area possesses two principal excitatory cells, the dentate granule cells and the mossy cells, and at least eight types of inhibitory interneurons (Morgan and Soltesz, 2008). In the present study, the DG was modeled as formed by four subpopulations of neurons, the granule cells (GC) and mossy cells (MC) (both excitatory) and the pyramidal basket cells (BC) and hilar cells (HC) (both inhibitory). A schematic representation is shown in Fig. 5, as can be observed, the inhibitory populations have feedback loops which is consistent with morphological studies and this recurrent inhibitory connection allows the model to have oscillations in the gamma band (30–70 Hz). MC connects to GC and vice-versa, BC connects principally to GC and HC connects mainly to GC and BC (continuous arrows represent excitatory connections and dashed arrows represent inhibitory connections). The equations that describe the dynamics of the model are a modified version of Wendling and coworkers' model (MolaeeArdekani et al., 2010) and are given by:

vgc ¼ igc

5

2

igc ¼ Aðt ÞaðSigmðvmc −ðC 4 vbc þ C 6 vhc ÞÞÞ−2aigc a vgc vmc ¼ imc 2 imc ¼ Aðt Þa pðt Þ þ C 2 Sigm C 1 vgc −2aimc −a vmc

C3

C4

C5

C6

C7

C8

a

b

g

A, B, G

173

28

34

86

54

3

14

100

50

300

Time varying

may not hold under certain pathological conditions; for example hippocampal sclerosis may change the number of synapses between excitatory and inhibitory populations (Marucci et al., 2010). Fig. 6 shows different types of activity that the model is able to reproduce and their respective frequency spectrum (compare with Fig. 3). CA1 model For the CA1 area, the model proposed by Wendling (MolaeeArdekani et al., 2010) was taken. This model was designed according to the organization and connectivity of CA1 area. It has one principal excitatory population (pyramidal cells), one excitatory interneuron population and two inhibitory interneuron populations, one for slow inhibition and other for fast inhibition. The model reliably reproduces different types of activity in the CA1 area, specifically, it is able to reproduce ictal, interictal and fast onset activity, for this reason this model was chosen for our experiments. The equations used for CA1 area are given below

2 ipy ¼ Aðt Þa Sigm ve −vsi −v f i −2aipy −a vpy

C2

28

vpy ¼ ipy

vbc ¼ ibc 2 ibc ¼ Bðt Þb Sigm C 3 vgc −C 8 vhc −2bibc −b vbc

C1

ð3Þ

vhc ¼ ihc 2 ihc ¼ Gðt Þg Sigm C 5 vgc −C 7 vbc −2bihc −b vhc yDG ¼ ðvmc −ðC 4 vbc þ C 6 vhc ÞÞ:

ve ¼ ie 2 ie ¼ Aðt Þa pðt Þ þ C 2 Sigm C 1 vpy −2aie −a ve

vsi ¼ isi 2 isi ¼ Bðt ÞbC 4 Sigm C 3 vpy −2bisi −b vsi

ð4Þ

vfi ¼ ifi

2 i f i ¼ Gðt ÞgC 7 Sigm C 5 vpy −C 6 vsi2 −2gi f i −b v f i

vsi2 ¼ isi2

2 isi2 ¼ Bðt ÞbSigm C 3 vpy −2bisi2 −b vsi2 yCA1 ¼ ðvmc −ðC 4 vbc þ C 6 vhc ÞÞ:

The states vx are the postsynaptic potentials of the population indicated by the suffix, i.e., vgc is the postsynaptic membrane potential of the granule cell population, yDG is the output of the system, it represents the signal that would be seen in an electrode. The parameters A, B, G are the excitatory, slow inhibitory and fast inhibitory gains respectively and are considered to be time varying; a, b, g are the lumped representation of the sum of the reciprocal of the time constant of passive membrane and all other spatially distributed delays in the dendritic network, p(t) is the input to the system and represents the incoming activity from external population activity. The parameters C1, C2,..., C8 represent the average number of synapses from one population to another. To be consistent with the structure of DG, relations among synapses from one population to another were averaged from results reported in the literature (Morgan and Soltesz, 2008). The values of the parameters are summarized in Table 1. It is important to notice that these relations

State and parameter estimation

Fig. 5. Schematic drawing of the dentate gyrus connection. Continuous arrows represent excitatory connections and dashed arrows represent inhibitory connections, rectangles represent population of specific types of neurons, the external input represents afferent activity from the entorhinal cortex. Abbreviations: MC, mossy cells; GC, granule cells; BC, pyramidal basket cells and HC, hilar cells.

Biological systems have the property of internal self-regulated environment and keep it stable. This property is called homeostasis, and it has been shown to play a central role in the regulation of normal activity within the brain (Chakravarthy et al., 2009a). In this work it was assumed that relations among internal parameters keep the brain from experiencing and spreading seizures, but when there is a pathological condition, homeostasis fails and these relations are no longer sustained; accordingly, it is hypothesized that when estimating key parameters in the model from a real pathological signal, parameter relations should change during an ictal event. We particularly focus on three parameters, the relation of the excitation and inhibition gains (Wendling et al., 2002) and the parameter K, the strength connection from a population in one area to a population in another area, i.e., from right DG to right CA1, or, from right DG to left DG. Here the estimation of states and

The states vpy, vsi and vfi represent postsynaptic potentials of the pyramidal cell (excitatory), slow inhibitory and fast inhibitory populations respectively, ve represents the postsynaptic potentials from excitatory interneuron population, the state vsi2 represents the fraction of the slow inhibitory population that projects to the fast inhibitory population. The variable yCA1 is the output of the model and is equivalent to the local field potential of an electrode. Parameters A, B and G represent excitation, slow inhibition and fast inhibition gains. The parameter values are given in Table 2.


6


1 Power (µV)

Amplitude (mV)

4 2 0 −2 0

1

2

0

3

10

20

30

40

50

60

1 Power (µV)

Amplitude (mV)

4 2 0

0.5

0

−2 0

1

2

3

10

20

30

40

50

60

1 Power (µV)

4 Amplitude (mV)

0.5

2 0 −2 0

1

2

3

0.5

0 0

10

20

30

40

50

60

Frequency (Hz)

Time (s)

Fig. 6. Three segments of simulated signals (left) and their respective frequency spectrum (right). This figure shows the different types of activity that the model is able to reproduce (compare to Fig. 3).

parameters of the model was carried out with the CKF which performs efficient joint state and parameter estimation (Togneri and Deng, 2003) and it is specially well suited for nonlinear systems (López-Cuevas et al., 2011). Briefly, the CKF is a Bayesian, recursive predict-update process. This recursion allows for on-line estimation. In the Bayesian nonlinear filtering paradigm it is necessary to solve integrals of the form: (nonlinear function) × (Gaussian density), which are usually intractable. The CKF overcomes this issue by evaluating integrals numerically using cubature rules (Arasaratnam and Haykin, 2009). For a detailed explanation of the method see (Havlicek et al., 2011). In general a nonlinear system can be put in terms of process and measurement equations as: xk ¼ f ðxk−1 Þ þ vk−1 yk ¼ hðxk Þ þ wk

ð5Þ

where xk ∈ ℝn is the state of the system at time k, yk is the measurement of the system or observed variable, vk, wk are zero-mean Gaussian noise and f(⋅) is a nonlinear function. The task of estimating the states of a nonlinear system consists thus in estimating xk using only the measurement, namely, the signals from the implanted microelectrodes. Additionally to the states, it is also required to estimate the parameters of the system, so that the concatenated state is described in the following way

xk f ðxk−1 Þ þ vk−1 ¼ pk pk−1 þ ξk−1 yk ¼ hðxk Þ þ ωk

result of imbalance between excitation and inhibition (Engel, 1996), in this work we choose to track the changes in the parameters A, B and G that can be considered as the amount of excitation and inhibition, respectively. In several studies it has been shown that an increase in excitatory postsynaptic currents (EPSC) and a decrease in inhibitory postsynaptic currents (IPSC) occur in pilocarpine treated rats (Zhan et al., 2010); probably because of synaptic reorganization, axon sprouting and loss of specific types of interneurons, besides, in (Morales-Villagrán et al., 2008a, 2008b) it was demonstrated that extracellular changes in glutamate and GABA occur during seizures; therefore, it is reasonable to expect a change in the parameters during SE. In the case of estimation of the states and parameter of the DG, vector xk, and pk take the following form 0

1 vpy B C B ipy C B C B v C B e C B C B ie C B C B vsi C C ðxk Þ ¼ B B i C B si C Bv C B fi C B C B ifi C B C Bv C @ si2 A i 0 si21 A ðpk Þ ¼ @ B A G

ð6Þ

ð7Þ

where p is the parameter vector of the system. This method is known as joint state and parameter estimation. Since epilepsy is believed to be the

Table 2 Parameters of the CA1 model. C1

C2

C3

C4

C5

C6

C7

a

b

g

A, B, G

135

108

33.75

33.75

13.5

13.5

108

100

28.5

200

Time varying

and the measurement yk corresponds to the signal in the implanted microelectrode for experiments 2 and 3, and to an artificially generated signal for experiment 1. Whether artificial or real signal, the measurement is the only information provided to the CKF, i.e., there is no information available on the states or parameters since it would require



having additional sensors or electrodes. The estimation for CA1 was performed in a similar way. Results

7

undergoes a bifurcation which resembles a transition from interictal to ictal state.

Experiment 2

In this work three different experiments were carried out: in the first one, signals were generated artificially with the present model and parameters were varied during simulation; the task in experiment 1 consisted of identifying correctly the parameters and the states of the model by using the CKF and with only the measured variable of the model and no extra information. The second experiment consisted of estimating states and parameters from real signals of intracranial microelectrodes recordings (see experimental setting section) and analyzing whether a change occurred in the parameters during SE and drug injection. The third experiment consisted of identifying strength coupling variation among populations in different microelectrode locations; special emphasis was put on populations from opposite hemispheres. The three experiments carried out and the respective results are described in the following section. Experiment 1 In the first case, parameter A was varied during the simulation of the model. The variation was not intended to have a biological meaning but to observe the performance of the method and the behavior of the model to variations of parameter A, for this experiment the model of DG was utilized. Fig. 7 shows the estimation of the output y and state vgc. It is worth noticing that the parameter was varied sinusoidally, from 0 to 20 s. Immediately after, abrupt changes in the parameter were introduced. As can be observed from Fig. 8, the estimator is capable of tracking the parameter variation in a suitable way. It is equally important to remark that the system, during this parameter variation,

In this experiment, the states and parameters of the neural mass model are estimated from real electrophysiological signals measured through microelectrodes implanted in experimental animals with SE induced by pilocarpine. The goal was to estimate the states and parameters of the DG and CA1 models and analyze changes occurring in the parameter during transition from basal state to SE. Figs. 9 and 10 show the identification of the right DG and left CA1 signals respectively. The figures show three sub-plots that correspond to the entire signal (top), selected segment A (middle) and selected segment B (bottom). The real signal is shown in black and the estimated signal in red. The mean squared errors (MSE) of the entire experiment are 0.0021 for the CA1 and 0.045 for the DG. The simulations were performed individually on each model, i.e., connections among CA1 and DG were not taken into account in this experiment. The evolution of the estimated parameters during all the experiments is shown in Fig. 11 for left DG. In the upper sub-plot, parameter A is shown in magenta color and the signal from the microelectrode is shown in gray, the actual signal was scaled for comparison purposes. The three horizontal black lines in the figures represent the mean of parameter A in the basal state (before pilocarpine injection) and the mean ± standard deviation (black dashed lines). As can be observed from Fig. 11, during basal activity, parameter A presented a standard deviation of 30% of the mean value, while during SE activity, the parameter presented variations larger than 200% of the mean value with respect to the basal state. It also can be observed that, in the transition period between the pilocarpine injection and SE onset, parameter A presented variations up to 80% of the mean value. Also, it is interesting to observe

a)

b) 20 y y

15

est

Amplitude (mV)

Amplitude (mV)

20

10 5 0 −5 0

15

y yest

10 5 0

5

10

15

20

25

30

35

40

−5 0

45

1

c)

2

3

d)

20

20

Vgc

Vgc V

gc est

Amplitude (mV)

Amplitude (mV)

V

10

0 0

5

10

15

20

4

25

Time (s)

30

35

40

45

gc est

10

0 0

1

2

3

4

Time (s)

Fig. 7. State estimation from experiment 1. This figure shows the estimation from simulated data, the black signal represents the simulated data and the red signals represent estimated data. Plots a and c show the signal corresponding to output y and state vgc. Plots b and d show an expanded view of the selected segments in plots a and c respectively.


8


Parameter A 27 real estimate

25 23

Amplitude (mV)

21 19 17 15 13 11 9

0

5

10

15

20

25

30

35

40

45

Time (s) Fig. 8. Parameter estimation from experiment 1. The parameter A was varied during the experiment, the simulated parameter is shown in blue and the estimated in red.

that at time 9000 s parameter A reaches its lowest value, just before the end of SE. In the bottom sub-plot of Fig. 11, parameters B (blue) and G (red) are shown. It can be observed that during the basal state parameter B, which represents inhibition from basket cells in the DG, has variations of the 10% of its mean value while in the transition period (after the pilocarpine injection) the value of parameter B presented variations up to 35% of the mean basal value while during the SE parameter B presents its larger variation, up to 80% of the mean basal value. In a similar manner, the variations of parameter G can be observed, which represents inhibition from hilar cells.

In Fig. 12, the mean values of parameters A, B and G are shown in the upper plot. The mean values were calculated over four distinct stages of the experiment: basal state, transition state (after pilocarpine), SE and cessation. This graphic illustrates the average changes that the different parameters present. In the lower subplot of Fig. 12, the ratio between excitation and inhibition is shown, specifically, the ratios A=B and A=G where Ā, B, G represent the mean values of A, B and G respectively. From this figure, it can be observed that inhibition from hilar cells (parameter G) presented its highest average value in the transition period (after pilocarpine and before SE onset), perhaps this increase in inhibition (with respect to the basal state) represents

Right Dentate Gyrus Signal Identification Amplitude (mV)

4

Segment B

Segment A

2 0 −2 −4 0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

646

647

648

649

650

9266

9267

9268

9269

9270

Segment A Amplitude (mV)

1 0.5 0 −0.5 −1

641

642

643

644

645 Segment B

Amplitude (mV)

2 1 0 −1 −2

9261

9262

9263

9264

9265

Time (s) Fig. 9. Experiment 2 results. Right DG signal identification. Real signal is shown in black and estimated signal in red. The upper plot shows the entire signal identification, the MSE is 0.045. The middle and bottom plots show a magnified view of the (randomly) selected segments A and B respectively.



9

Left CA1 signal identification

Amplitude (mV)

2.5 Real Estimated

Segment B

Segment A

0

−2.5 0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

Segment A

Amplitude (mV)

0.5

0

−0.5

916

917

918

919

920

921

922

923

924

925

926

Segment B

Amplitude (mV)

0.5

0

−0.5 5352

5353

5354

5355 Time (s)

5356

5357

5358

Fig. 10. Experiment 2 results. Left CA1 signal identification. Real signal is shown in black and estimated signal in red. The upper plot shows the entire signal identification, the MSE is 0.0021. The middle and bottom plots show a magnified view of the (randomly) selected segments A and B respectively.

a homeostasis mechanism to prevent SE propagation and it may be the reason that the SE onset in the left DG is approximately 1 h after pilocarpine injection. Also, it can be observed that during SE the average Pilocarpine Injection

value of G is lower than during the basal state, in addition the ratio A=G presents its highest value during this period which means that there is an increase in the excitation–inhibition ratio. Interestingly the largest

Parameter A estimation Left DG Parameter A

Amplitude (mV)

20 15 10 5 0 −5 −10

Basal activity

Status epilepticus

After pilocarpine

Cessation

−15 1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

Parameters B and G estimation 70 Parameter G Parameter B Left DG

Amplitude (mV)

60 50 40 30 20 10 0 −10

1000

Status epilepticus

After pilocarpine

Basal activity 2000

3000

4000

5000 Time (s)

6000

7000

8000

Cessation 9000

10000

Fig. 11. Experiment 2 results. Parameter identification. In this experiment three parameters were estimated from real data in the left DG, parameters A, B and G. Upper subplot, Parameter A (magenta). Lower subplot: Parameter B (red) and Parameter G (blue). The mean and standard deviations of the parameters during the basal state are shown for comparison, mean is presented with black lines and the mean ± the standard deviation in black dashed lines. The signal of the left DG is shown in gray color.


10

A. López-Cuevas et al. / NeuroImage xxx (2015) xxx–xxx Pilocarpine Injection

Averaged Parameters

70 Parameter G Parameter B Parameter A

60

Amplitude

50 40 30 20 10 0 −10

Basal activity 1000

Status epilepticus

After pilocarpine 2000

3000

4000

5000

6000

7000

8000

Cessation 9000

10000

Averaged Excitation−Inhibition ratio 0.4 A/B A/G

Amplitude

0.3 0.2 0.1 0 −0.1 1000

Status epilepticus

After pilocarpine

Basal activity −0.2

2000

3000

4000

5000 Time (s)

6000

7000

8000

Cessation 9000

10000

Fig. 12. Experiment 2 results. Segment-wise averaged parameters. In this figure the parameters were averaged over four different segments: basal state, transition state, SE and cessation. Upper subplot: the averaged parameters are shown, parameter Ā (magenta), parameter B (red), parameter G (blue). Lower subplot: averaged excitation inhibition ratio. A=B is shown in red and A=G in blue color.

average value for parameter B (inhibition from basket cells) is presented during SE, but perhaps is not strong enough to stop seizures. It appears that cessation is achieved because a reduction in the excitation, in fact the ratio A=B and A=G reaches the lowest average value in this period. Experiment 3 Another characteristic of epilepsy is that populations of neurons can synchronize even over long distances or separate hemispheres. It has been reported that a possible mechanism for such synchronization is the strengthened coupling between neurons in distant populations. In this experiment, changes in a coupling parameter between populations from different areas were analyzed when signals were taken from two different microelectrodes. For this experiment 2 populations were modeled to be coupled by a parameter K (Wendling et al., 2000; Chakravarthy et al., 2009b). In this case the coupling parameter K is modeled as follows right right ¼ igc right right left right right right 2 right −2aigc −a vgc igc ¼ Aðt Þa Sigm vmc þ KyDG − C 4 vbc þ C 6 vhc right right vmc ¼ imc right right right 2 right −2aimc −a vmc imc ¼ Aðt Þa pðt Þ þ C 2 Sigm C 1 vgc right right vbc ¼ ibc right right right right 2 right ibc ¼ Bðt Þb Sigm C 3 vgc −C 8 vhc −2bibc −b vbc right right vhc ¼ ihc right right right right 2 right −2bihc −b vhc ihc ¼ Gðt Þg Sigm C 5 vgc −C 7 vbc right right right right yDG ¼ vmc − C 4 vbc þ C 6 vhc

vgc

ð8Þ

left

left

vgc ¼ igc

left left right left left left 2 left −2aigc −a vgc igc ¼ Aðt Þa Sigm vmc þ KyDG − C 4 vbc þ C 6 vhc

left

left

vmc ¼ imc

left left left 2 left imc ¼ Aðt Þa pðt Þ þ C 2 Sigm C 1 vgc −2aimc −a vmc left left vbc ¼ ibc left left left left 2 left ibc ¼ Bðt Þb Sigm C 3 vgc −C 8 vhc −2bibc −b vbc left left vhc ¼ ihc left left left left 2 left ihc ¼ Gðt Þg Sigm C 5 vgc −C 7 vbc −2bihc −b vhc left left left left yDG ¼ vmc − C 4 vbc þ C 6 vhc :

ð9Þ

The superscripts right and left indicate the hemisphere under consideration. Fig. 13 shows parameter K from the right DG to the left DG. In the first and last subplots, the signals of the right and left DG are shown respectively. Coupling between the two populations is shown in the middle subplot. It is interesting to note that the coupling between the two populations increased during the propagation of the seizure activity to the left hemisphere, even when in the right hemisphere the seizure activity started earlier. The coupling from Right CA1 to left CA1 was estimated in a similar manner, Fig. 14 shows the coupling of the left and right CA1. Discussion Epilepsy is a disorder that affects 1% of the world population. Despite the advent of new antiepileptic drugs in recent years, approximately 30% of epileptic patients are resistant to medication (Stacey and Litt, 2008). Regarding this matter, neuroscientists have been developing alternative therapies. One of these therapies is deep brain stimulation,



11

Right Dentate Gyrus Amplitude (mV)

4 2 0 −2 1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

Amplitude (dimensionless)

Coupling parameter K Coupling Mean (basal state) Mean ± std (basal state)

200 100 0 1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

7000

8000

9000

10000

Left Dentate Gyrus Amplitude (mV)

4 2 0 −2 1000

2000

3000

4000

5000

6000

Time (s) Fig. 13. Experiment 3 results. Coupling between two different areas. In this figure, the coupling parameter K is estimated for the right and left DG. The middle sub-plot shows the estimated parameter K (red). The upper and lower sub-plots show the signal from right DG and left DG respectively (blue).

which consist of stimulating specific areas of the brain by electrical pulses in order to reduce seizure appearance. A major drawback in the existing devices is that they function in open loop mode, meaning that a fixed rate of stimulating pulses is pre-programmed, independently whether the patient is undergoing a seizure or not (Cox et al., 2014; Han et al., 2014; McIntyre et al., 2014). An alternative approach consists of developing a device that permits stimulation only when needed and with specific intensity. This

is known as a closed loop approach. In order to achieve closed loop operation on the device, it is necessary that a model describes, in mathematical terms, the electrical activity of a specific area of the brain able to generate seizures in order to apply an adequate stimulation. The results of this work are suitable for this task, in the sense that a mathematical model is provided along with a method to fit the model online. The results here can also be used in designing implantable devices that deliver medication at the right time, given that

Amplitude (mV)

Right CA1

5 0 −5 0

1000

2000

3000

4000

5000

6000

7000

8000

9000

10000

7000

8000

9000

10000

7000

8000

9000

10000

Amplitude (dimensionless)

Coupling parameter K

400 200 0 0

1000

2000

3000

4000

5000

6000

Amplitude (mV)

Left Dentate Gyrus

5 0 −5 0

1000

2000

3000

4000

5000

6000

Time (s) Fig. 14. Experiment 3 results. Coupling between two different areas. In this figure, the coupling parameter K is estimated for the right and left CA1. The middle sub-plot shows the estimated parameter K (red). The upper and lower sub-plots show the signal from right CA1 and left CA1 respectively (blue).


12


parameters related to overall excitation or inhibition can be tracked online (Lohman et al., 2005). In the present study, we applied the CKF estimation approach on data from EEG recording episodes. The proposed method is able to track changes in the states and parameters of the model from electrophysiological signals obtained from an experimental model of SE. Although the time scale of the SE is, by definition, longer than 30 min, which is large in comparison to the time scale of a regular seizure (in the range of seconds or minutes), the method is still able to track changes in shorter time scales, as shown in experiment 1 of this work, where the dynamics of the simulated parameters were fast (see Fig. 8). Therefore, the model proposed here, may serve, in principle, to estimate parameters not only during SE but in regular seizures as well, although more experimental data and further experiments are needed. In addition, according to EEG and histopathological studies a seizure lasting 5 or more minutes is enough to be considered as serious because the SE and the chances for self-termination are too small, therefore it requires the same treatment as seizures lasting 30 min (SE) (Meierkord et al., 2010; Millikan et al., 2009). With regard to excitation–inhibition balance, it is considered that the pathophysiological mechanism that underlies SE induced by pilocarpine is the activation of N-methyl D-aspartate (NMDA) receptors initiated via a muscarinic type 1 receptor that increases intracellular calcium and release of glutamate from glutamatergic neurons generating an excitotoxicity process that leads to neuronal death (Smolders et al., 1997). Also, a decrease has been observed in the synaptic inhibition during SE, such as a decrease in the number of postsynaptic gammaaminobutyric acid (GABA) type A receptors, whereas increases in postsynaptic glutamatergic receptors or release of glutamate further upsets the balance between excitation and inhibition (Medina-Ceja et al., 2000; Naylor, 2010; Reddy and Kuruba, 2013). These previous data agree with the present results in which an increase in the ratio of excitation/inhibition was observed during SE, parameter A=G. In addition, it is interesting to observe that previous to the onset of SE parameter G, inhibition from hilar interneurons increases in order to stop seizures developed by the SE (fast inhibition), when this compensatory mechanism does not work a second attempt is triggered to stop SE, the slow inhibition from basket cells (interneurons from DG) parameter B. It has been found that during the SE induced by pilocarpine, an increase in the tonic current of GABA and change in the reversal potential for GABA current to depolarizing maintain the homeostasis of fast-spiking basket cell activity in order to limit increases in dentate excitability (Yu et al., 2013), this physiological behavior could be happening during the second attempt to stop SE. Therefore, the failure to stop seizure development by the SE could be due to limited availability of GABA A receptors, as previous data show that the GABA A receptors move from the synaptic membrane to the cytoplasm where they are functionally inactive (Chen et al., 2007) and the subunit composition of glutamate and GABA receptors in patient and animal models of status epilepticus promote self-sustaining seizures (Loddenkemper et al., 2014). Also this fact explains the tendency of seizures to resist inhibition (slow and fast) and become self-sustaining during the SE.

Conclusions The estimation of the states and parameters of the neural mass models, both simulated and during real intracranial EEG recordings obtained during SE induced by pilocarpine is studied in this work. This estimation is based on the cubature Kalman filter which is especially well suited for nonlinear systems. Both, synthetic and real data were utilized in the experiments to validate the model. The estimation performance is good, given that the MSE of the estimated signals was below 0.05. It is worth noticing that this method can be implemented in real time, which is particularly important for the design of stimulation devices that attempt to stop seizures. The use of a neural mass model

could allow for the design of feedback control strategies able to deliver electrical stimulation or anti-epileptic drugs in an optimal way. In addition, the parameters of the model related to excitation and inhibition were also estimated from synthetic and real data. It was observed that during basal activity (before pilocarpine administration) the parameters presented a standard deviation below 30% of the mean value, while during SE activity, the parameters presented variations larger than 200% of the mean value with respect to the basal state. The ratio of excitation–inhibition, A=G, increased during SE activity by 80% with respect to the transition state, and reaches the lowest value during cessation. The study of internal changes occurring during seizure activity can lead to better understanding of the mechanisms of epilepsy and to improve its treatment. Acknowledgments This work was partially supported by the CONACYT grants 127858, 215648 and 106179. References Amaral, D., Scharfman, H., Lavenex, P., 2007. The dentate gyrus: fundamental neuroanatomical organization (dentate gyrus for dummies). Prog. Brain Res. 163, 3–22. Angus, L.H., Parsons, L.M., 2008. Epilepsy: epidemiology, classification and natural history. Epilepsia 36, 571–578. Arasaratnam, I., Haykin, S., 2009. Cubature Kalman filters. IEEE Trans. Autom. Control 54, 1254–1269. Brodie, M.J., 1990. Status epilepticus in adults. Lancet 336, 551–552. Cavalheiro, E., Santos, N., Priel, M., 1996. The pilocarpine model of epilepsy in mice. Epilepsia 37, 1015–1019. Chakravarthy, N., Tsakalis, K., Sabesan, S., Iasemidis, L., 2009a. Homeostasis of brain dynamics in epilepsy: a feedback control systems perspective of seizures. Ann. Biomed. Eng. 37 (3), 565–585. Chakravarthy, N., Sabesan, S., Tsakalis, K., Iasemidis, L., 2009b. Controlling epileptic seizures in a neural mass model. J. Comb. Optim. 17, 98–116. Chen, J.W., Naylor, D.E., Wasterlain, C.G., 2007. Advances in the pathophysiology of status epilepticus. Acta Neurol. Scand. 115 (4 Suppl.), 7–15. Commission on Classification and Terminology of the International League Against Epilepsy, 1981. Proposal for revised clinical and electroencephalographic classification of epileptic seizures. Epilepsia 22, 489–501. Cox, J.H., Seri, S., Cavanna, A.E., 2014. Clinical utility of implantable neurostimulation devices as adjunctive treatment of uncontrolled seizures. Neuropsychiatr. Dis. Treat. 10, 2191–2200. David, O., Friston, K., 2003. A neural mass model for MEG/EEG: coupling and neuronal dynamics. NeuroImage 20, 1743–1755. Engel Jr., J., 1996. Excitation and inhibition in epilepsy. Can. J. Neurol. Sci. 23 (3), 167–174. Engel Jr., J., 2006. Report of the ILAE classification core group. Epilepsia 47, 1558. Fisher, R.S., Van Emde, B., Blume, W., Elger, C.E., Genton, P., Lee, P., Engel Jr., J., 2005. Epileptic seizures and epilepsy: definitions proposed by the International League Against Epilepsy (ILAE) and the International Bureau for Epilepsy. Epilepsia 46, 470. Freeman, W.J., 1987. Simulation of chaotic EEG patterns with a dynamic model of the olfactory system. Biol. Cybern. 56, 139–150. Friston, K.J., Harrison, L., Penny, W., 2003. Dynamic causal modelling. NeuroImage 19 (4), 1273–1302. Goodfellow, M., Schindlerb, K., Baier, G., 2012. Self-organised transients in a neural mass model of epileptogenic tissue dynamics. NeuroImage 59 (3), 2644–2660. Han, C., Hu, W., Stead, M., Zhang, T., Zhang, J., Worrell, G.A., Meng, F., 2014. Electrical stimulation of hippocampus for the treatment of refractory temporal lobe epilepsy. Brain Res. Bull. 109C, 13–21. Havlicek, M., Friston, K.J., Jan, J., Brazdil, M., Calhoun, V., 2011. Dynamic modeling of neuronal responses in fMRI using cubature Kalman filtering. NeuroImage 56 (4), 2109–2128. Huang, N., Shen, Z., Long, S., Wu, M., Shih, H., Zheng, Q., Yen, N., Tung, C., Liu, H., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. A 454, 903–995. Jansen, B.H., Rit, V., 1995. Electroencephalogram and visual evoked potentials generation in a mathematical model of coupled cortical column. Biol. Cybern. 68, 275–283. Jenssen, S., Gracely, E.J., Sperling, M.R., 2006. How long do most seizures last? A systematic comparison of seizures recorded in the epilepsy monitoring unit. Epilepsia 47, 1499–1503. Loddenkemper, T., Talos, D.M., Cleary, R.T., Joseph, A., Sanchez Fernandez, I., Alexopoulos, A., Kotagal, P., Najm, I., Jensen, F.E., 2014. Subunit composition of glutamate and gamma-aminobutyric acid receptors in status epilepticus. Epilepsy Res. 108 (4), 605–615. Lohman, R.J., et al., 2005. Validation of a method for localised microinjection of drugs into thalamic subregions in rats for epilepsy pharmacological studies. J. Neurosci. Methods 146, 191–197. Lopes da Silva, F., Blanes, W., Kalitzin, S., Parra, J., Suffczynski, P., Velis, D., 2003. Dynamical diseases of brain systems: different routes to epileptic seizures. IEEE Trans. Biomed. Eng. 50 (5), 540–548.


A. López-Cuevas et al. / NeuroImage xxx (2015) xxx–xxx López-Cuevas, A., Castillo-Toledo, B., Medina-Ceja, L., Ventura-Mejía, C., 2011. Nonlinear identification of epileptiform signals using an artificial neural network trained with the cubature Kalman filter. IEEE World Congress on Engineering and Technology (CET), Shanghai, p. 21012. López-Cuevas, A., Castillo-Toledo, B., Medina-Ceja, L., Ventura-Mejía, C., Pardo-Peña, K., 2013. An algorithm for on-line detection of high frequency oscillations related to epilepsy. Comput. Methods Prog. Biomed. 110 (3), 354–360. Lytton, W., 2008. Computer modelling of epilepsy. Nat. Rev. Neurosci. 9, 626–637. Marucci, G., Rubboli, G., Giulioni, M., 2010. Role of dentate gyrus alterations in mesial temporal sclerosis. Clin. Neuropathol. 29 (1), 32–35. McIntyre, C.C., Chaturvedi, A., Shamir, R.R., Lempka, S.F., 2014. Engineering the next generation of clinical deep brain stimulation technology. Brain Stimul. http://dx.doi.org/ 10.1016/j.brs.2014.07.039. Medina-Ceja, L., Morales-Villagran, A., Tapia, R., 2000. Action of 4-aminopyridine on extracellular amino acids in hippocampus and entorhinal cortex: a dual microdialysis and electroencephalographic study in awake rats. Brain Res. Bull. 53, 255–262. Medina-Ceja, L., Pardo-Peña, K., Ventura-Mejía, C., 2014. Evaluation of behavioral parameters and mortality in a model of temporal lobe epilepsy induced by intracerebroventricular pilocarpine administration. Neuroreport. 25 (11), 875–879. Meierkord, H., Boon, P., Gocke, K., Shorvon, S., Tinuper, S., Holtkamp, M., 2010. EFNS guidelines on the management of status epilepticus in adults. Eur. J. Neurol. 17, 348–355. Millikan, D., Rice, B., Silbergleit, R., 2009. Emergency treatment of status epilepticus: current thinking. Emerg. Med. Clin. North Am. 27, 101–113. Molaee-Ardekani, Behnam, Benquet, Pascal, Bartolomei, Fabrice, Wendling, Fabrice, 2010. Computational modeling of high-frequency oscillations at the onset of neocortical partial seizures: from altered structure to dysfunction. NeuroImage 52, 1109–1122. Morales-Villagrán, A., Sandoval-Salazar, C., Medina-Ceja, L., 2008a. An analytical flow injection system to measure glutamate in microdialysis samples based on an enzymatic reaction and electrochemical detection. Neurochem. Res. 33, 1592–1598. Morales-Villagrán, A., Medina-Ceja, L., López-Pérez, S.J., 2008b. Simultaneous glutamate and EEG activity measurements during seizures in rat hippocampal region with the use of an electrochemical biosensor. J. Neurosci. Methods 168 (1), 48–53. Morgan, R., Soltesz, I., 2008. Nonrandom connectivity of the epileptic dentate gyrus predicts a major role for neuronal hubs in seizures. PNAS 105, 6179–6184. Naylor, D.E., 2010. Glutamate and GABA in the balance: convergent pathways sustain seizures during status epilepticus. Epilepsia 51 (Suppl. 3), 106–109.

13

Perucca, E., 2007. Development of new antiepileptic drugs: challenges, incentives, and recent advances. Lancet Neurol. 6, 793–804. Racine, R.J., 1972. Modification of seizure activity by electrical stimulation. II. Motor seizure. Electroencephalogr. Clin. Neurophysiol. 32 (3), 281–294. Reddy, D.S., Kuruba, R., 2013. Experimental models of status epilepticus and neuronal injury for evaluation of therapeutic interventions. Int. J. Mol. Sci. 14, 18284–18318. Smolders, I1., Khan, G.M., Manil, J., Ebinger, G., Michotte, Y., 1997. NMDA receptormediated pilocarpine-induced seizures: characterization in freely moving rats by microdialysis. Br. J. Pharmacol. 121 (6), 1171–1179. Spiegler, A., Kiebel, S., Atay, F., Knösche, T., 2010. Bifurcation analysis of neural mass models: impact of extrinsic inputs and dendritic time constants. NeuroImage 52, 1041–1058. Stacey, W., Litt, B., 2008. Technology insight: neuroengineering and epilepsy, designing devices for seizure control. Nat. Clin. Pract. Neurol. 4 (4), 190–201. Togneri, R., Deng, L., 2003. Joint state and parameter estimation for a target-directed nonlinear dynamic system model. IEEE Trans. Signal Process. 51 (12), 3061–3069. Traub, R., Whittington, M., Cunningham, M., 2010. Epileptic fast oscillations and synchrony in vitro. Epilepsia 51, 28. Wendling, F., Bellanger, J.J., Bartolomei, F., Chauvel, P., 2000. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biol. Cybern. 83 (4), 367–378. Wendling, F., Bartolomei, F., Bellanger, J., Chauvel, P., 2002. Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Eur. J. Neurosci. 15, 1499–1508. Wendling, F., Hernandez, A., Bellanger, J., Chauvel, P., Bartolomei, F., 2005. Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG. J. Clin. Neurophysiol. 22 (5), 343–356. Wilson, H., Cowan, J., 1972. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12. Yu, J., Proddutur, A., Elgammal, F.S., Ito, T., Santhakumar, V., 2013. Status epilepticus enhances tonic GABA currents and fast-spiking basket cells depolarizes GABA reversal potential in dentate fast-spiking basket cells. J. Neurophysiol. 109, 1746–1763. Zavaglia, M., Astolfi, L., Babiloni, F., Ursino, M., 2006. A neural mass model for the simulation of cortical activity estimated from high resolution EEG during cognitive or motor tasks. J. Neurosci. Methods 157, 317–329. Zhan, R., Timofeeva, O., Nadler, J., 2010. High ratio of synaptic excitation to synaptic inhibition in hilar ectopic granule cells of pilocarpine-treated rats. J. Neurophysiol. 104 (6), 3293–3304.


Neural population coding: combining insights from microscopic and mass signals.

Status epilepticus: a new rodent model.

Basin structure of optimization based state and parameter estimation.

Status epilepticus and refractory status epilepticus management.

Linking a neural mass model with a 3D model of the human brain to reproduce EEG signals.

A rodent model of focally evoked self-sustaining status epilepticus.

[Electrographic "status epilepticus" during afternoon sleep].

Status epilepticus developing during lacosamide monotherapy.

Estimating a dynamic state to relate neural spiking activity to behavioral signals during cognitive tasks.

Super-refractory nonconvulsive status epilepticus secondary to fat embolism: A clinical, electrophysiological, and pathological study.

Variational Bayesian Parameter Estimation Techniques for the General Linear Model.

Artificial neural networks in bioprocess state estimation.

Enzyme histochemistry of the rat hippocampus during experimental status epilepticus.

A simulation model of AIDS in San Francisco: I. Model formulation and parameter estimation.

Identifiability and parameter estimation.

Transport parameter estimation from lymph measurements and the Patlak equation.

Model-based MR parameter mapping with sparsity constraints: parameter estimation and performance bounds.

A Pharmacokinetics-Neural Mass Model (PK-NMM) for the Simulation of EEG Activity during Propofol Anesthesia.

Emergency Department Management of Pediatric Unprovoked Seizures and Status Epilepticus in the State of Illinois.

A hybrid parameter estimation algorithm for beta mixtures and applications to methylation state classification.

Improvement in Neural Respiratory Drive Estimation From Diaphragm Electromyographic Signals Using Fixed Sample Entropy.

Neural network programming in bioprocess variable estimation and state prediction.

Neural correlates of success and failure signals during neurofeedback learning.

Underestimated Rate of Status Epilepticus according to the Traditional Definition of Status Epilepticus.