Patient-Specific Neural Mass Modeling - Stochastic ... - Levin Kuhlmann

2 downloads 0 Views 857KB Size Report
3010, Australia, and The Bionics Institute, 384-388 Albert St. East Melbourne, VIC,. 3002, Australia. E-mail: {deanrf,levink,grayden,dnesic}@unimelb.edu.au, ...
November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

1

Patient-Specific Neural Mass Modeling - Stochastic and Deterministic Methods D.R. Freestone, L. Kuhlmann, M. Chong, D. Neˇsi´ c and D.B. Grayden NeuroEngineering Laboratory, Department of Electrical and Electronic Engineering, and The Center for Neural Engineering, The University of Melbourne, Parkville, VIC, 3010, Australia, and The Bionics Institute, 384-388 Albert St. East Melbourne, VIC, 3002, Australia E-mail: {deanrf,levink,grayden,dnesic}@unimelb.edu.au, [email protected] P. Aram Department of Automatic Control and Systems Engineering, The University of Sheffield, Mappin Street, Sheffield, S1 3JD, United Kingdom E-mail: [email protected] R. Postoyan CRAN-ENSEM - 2 avenue de la fort de Haye, 54516 Vandoeuvre-ls-Nancy, France E-mail: [email protected] M.J. Cook Department of Neurology, St. Vincent’s Hospital Melbourne, Fitzroy VIC 3065, Australia E-mail: [email protected] Deterministic and stochastic methods for online state and parameter estimation for neural mass models are presented and applied to synthetic and real seizure electrocorticographic signals in order to determine underlying brain changes that cannot easily be measured. The first ever online estimation of neural mass model parameters from real seizure data is presented. It is shown that parameter changes occur that are consistent with expected brain changes underlying seizures, such as increases in postsynaptic potential amplitudes, increases in the inhibitory postsynaptic time-constant and decreases in the firing threshold at seizure onset, as well as increases in the firing threshold as the seizure progresses towards termination. In addition, the deterministic and stochastic estimation methods are compared and contrasted. This work represents an important foundation for the development of biologically-inspired methods to image underlying brain changes and to develop improved methods

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

2 for neurological monitoring, control and treatment. Keywords: Estimation; Neural Mass Model; EEG; Epilepsy

1. Introduction This chapter describes two methods for forming patient-specific mesoscopic neural mass models: a stochastic method and a deterministic method. The aim of these methods is to form a bridge between clinical and computational neuroscience that will facilitate the application of control engineering tools1 to develop new therapies for neurological conditions, such as epilepsy. The human brain is arguably the most complicated system known to man. The development of a complete theory of its function is one of the greatest challenges faced by researchers today. To address this challenge, researchers have used both theoretical and experimental frameworks to develop and test hypotheses. Experimental frameworks are typically designed to uncover causal relationships between the system’s properties or parameters and its function (reverse engineering the brain). Theoretical studies are typically used in one of two ways. The first is to explain data acquired in an experiment and the second is to predict system behaviour. This work is aimed at developing two more applications for theoretical studies and computational models. The first is to use the models as filters, providing insight into the physiology of the brain by estimating the states (i.e. neural activity) and the parameters (e.g. synaptic strength) of the modelled physiology. The second is to use the model as a tool for developing new therapies, where they can be used to provide feedback to a system that can systematically deliver an intervention (e.g. electrical stimulation) in a robust, controlled manner. Presently, it is almost half a century since Hodgkin and Huxley shared the Nobel prize (in 1963 with Eccles for Physiology and Medicine) for their influential model that helped to establish the field of theoretical neuroscience. Over this time period, both experimental and theoretical methods have developed considerably. Experimental approaches have been able to isolate neural function using various forms of manipulation in greater detail. In parallel, theoretical and computational models have been able to explain and provide new hypotheses for an ever increasing assortment of neural phenomena. Through the development of these models, a fundamental set of parametric equations has been established that explain neuronal responses to sensory input at varying spatiotemporal scales, from small patches of membrane to networks of neural ensembles.

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

3

Although the predictive explanatory power of the theoretical models is rather vast when they are tuned to mimic experimental conditions, a major limitation has been in establishing a rigorous method for choosing model parameters in more general situations. For example, over the last decade, it has come to light that there is significant variability in neural systems across subjects despite seemingly similar network activity.2 We expect an analogous situation in disease states and pathological activity such as epileptic seizures. Therefore, for models to be clinically useful, they must be subject-specific. By assimilating experimental or clinical measurements with mathematical models, one can infer unmeasured or latent system properties or parameters from standard clinical (electrophysiological) recordings. Inference of patient-specific parameters has the potential to revolutionise the treatment of disease. Typically, one can measure the blood-oxygen level dependent signal through functional magnetic resonance imaging (fMRI) or the electromagnetic fields of the brain through electroencephalography (EEG) or magnetoencephalography (MEG). In clinical neurology, these measurement modalities give a ‘where and when’ indication of the presence of a pathology, but provide little information on the causative mechanisms for disease. A classic example of this is in epilepsy management, where the disease is diagnosed using EEG by identifying electrographic seizures. If electrical seizures are recorded with EEG, a medication plan is prescribed. The choice of medication is not based on the data acquired with the EEG, but based upon the clinician’s experience and the patient’s circumstance, then the treatment plan is often modified based on trial and error. This process of searching for the best drug combination is often required and similar electrographic seizures can have fundamentally different mechanisms of initiation.3 This means that epilepsy is difficult to treat. Ideally, in epilepsy monitoring we would like to image the concentration of ions, the synaptic dynamics, the connectivity strength and structure, or other parameters to better inform a treatment plan by understanding physiological changes that lead to seizures. This is also highly desirable for epileptic seizure prediction which has shown limited progress to date.4,5 The signal processing and control theory literature generally provides two approaches to state and parameter estimation of models using measured data (otherwise known in engineering as ‘system identification’): stochastic and deterministic. In deterministic approaches,6–8 there is no universal technique that can be applied to a given general nonlinear system. Rather, the synthesis of deterministic approaches is often more specific to the form

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

4

of the system being estimated. Often deterministic assumptions are initially quite strict when first devising an estimator; however, these assumptions can usually be relaxed in order to deal with noise issues. Stochastic approaches1,9,10 on the other hand can be applied more easily to arbitrary systems; however, these approaches depend on initialization. Moreover, the convergence of the estimates to the true values is not guaranteed for every trajectory. In this chapter, we first describe the basic structure of neural mass models and give specific reference to the formulation of Jansen and Rit with which we will perform estimation.11 Then we show how two different estimation methods, one stochastic and one deterministic, can perform state and parameter estimation to illustrate their usefulness in inferring underlying and unmeasured physiological variables using only limited physiological measurements. Specifically, we apply a deterministic approach, referred to as an adaptive observer,12,13 to estimation of the states and parameters using simulated data. Then we apply a stochastic approach based on the unscented Kalman filter (UKF)14,15 to estimation of the states and parameters using both simulated data and real electrocorticography (ECoG) data of a seizure. Finally, we contrast and compare the stochastic and deterministic approaches in the discussion. 2. Neural Mass Model To define a standard neural mass model, we begin by defining the postsynaptic potential of population n as a result of an input firing rate from population m as their convolution Z αmn t hmn (t − t0 )gm (vm (t0 )) dt0 (1) vn (t) =vr,n + τmn −∞ Z αmn t vn (t) − vr,n = hmn (t − t0 )gm (vm (t0 )) dt0 (2) τmn −∞ Z αmn t v˜n (t) = hmn (t − t0 )gm (vm (t0 )) dt0 , (3) τmn −∞ where αmn is the gain for the post-synaptic response kernel, denoted by hmn (t), from neural population m to n, and τmn is the membrane time constant. Typically, αmn and τmn are constants (particularly for current based synapses), but later in this chapter we will consider them as time varying quantities. Also, gm (vm (t0 )) describes the input firing rate as a function of the pre-synaptic membrane potential. The resting membrane potential of the post-synaptic population is denoted by vr,n , vn (t) is the

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

5

post-synaptic membrane potential and v˜n (t) is the deviation of the membrane from the resting potential. For the network of neural masses that we are considering in the chapter, the index n (post-synaptic) may represent either the pyramidal (p), excitatory interneuron (spiny stellate) (e), external (x) or inhibitory interneuron (i) populations. The post-synaptic response kernel, hmn (t), typically takes one of three different forms: one first-order and two second-order. The first-order form has an instantaneous rise time and a decay defined by a single time constant.16 The second-order kernels have finite rise and decay times, with the difference being with one form having separate time constants17 (biexponential) for the rise (synaptic time constant) and decay (membrane time constant), whereas, the other form is defined using a single time constant variable by18   t , (4) hmn (t) = η(t)t exp − τmn where η(t) is the Heaviside step function. This is the form that we shall use; however the framework holds for other forms. This convolution can conveniently be written as αmn gm (vm (t)), (5) D˜ vn (t) = τmn where the linear differential operator is D=

d2 2 d 1 + + 2 . 2 dt τmn dt τmn

(6)

This allows the dynamics of the neural mass to be described by the differential equation d2 v˜n (t) 2 d˜ vn (t) 1 αmn + + 2 v˜n (t) = gm (vm (t)). dt2 τmn dt τmn τmn

(7)

This second-order ODE can be written as two coupled first-order ODEs by defining d˜ vn (t) . (8) dt Recasting the system in this way allows formation of a state-space model in a canonical format. This gives the system: zn (t) =

d˜ vn (t) =zn (t) dt 2 1 dzn (t) αmn = gm (vm (t)) − zn (t) − 2 v˜n (t). dt τmn τmn τmn

(9) (10)

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

6

There is a sigmoidal relationship between the mean membrane potential and firing rate of each of the populations. This sigmoid nonlinearity may take different forms, for example the cumulative density function (error function) or the logistic / hyperbolic tangent. Typically, the logistic function form is used, which is defined by 1 1 + exp (ςn (v0n − v˜n (t))) 1 g (˜ vn (t)) = 1 + exp (ςn (v0n + vrn − vn (t))) 1 g (vn (t)) = , 1 + exp (ςn (˜ v0n − vn (t))) g (˜ vn (t)) =

(11) (12) (13)

1 α τ h(t) τ ∝ α

0.5 0

0

mean input firing rate

0.05 Time (s)

0.1

Firing Rate (APs/s)

Amplitude (mV)

where v˜0n = v0n + vrn . Note that in this formulation, we are absorbing the maximal firing rate, which is typically a linear coefficient on the sigmoid, into the PSP gain (αmn ). This removes a redundant parameter that can not be recovered by estimation methods. The quantities ςn and v0n describe the slope of the sigmoid (variance of firing thresholds within the populations) and the mean firing threshold, respectively. These quantities are usually considered as constants, but again later in this chapter they will be treated as time varying. The parameter v˜0n describes the deviation from the resting membrane potential, which becomes our lumped threshold parameter. For ease of notation we can drop the tilde remembering the resting membrane potential resides within this term. Figure 1 depicts a standard neural mass.

1 0.5 g (v ) v0 0

post-synaptic potential

0

10 PSP (mV)

20

mean output firing rate

Fig. 1. Graphical Representation of the Neural Mass Model. The neural mass model converts an input firing rate to a mean post-synaptic potential by a convolution with the post-synaptic response kernel. The membrane potential is converted to output firing rate by the sigmoidal activation function.

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

7

This neural mass maps from a mean pre-synaptic firing rate to a postsynaptic mean membrane potential, which in turn determines the output firing rate of the post-synaptic population. The terms that are usually considered parameters of the model include τ , α, v0 , and ς. These can be set to model different neural populations, such as pyramidal neurons, spiny stellate cells, and fast and slow inhibitory interneurons (GABAa and GABAb ). The neural populations can then be configured to represent the circuitry of a cortical column and networks of cortical columns. Various kinds of neural mass models have been developed.11,19–21 The parameters of the neural masses not only define the population type, but also the behaviour the model exhibits. For example, for certain parameter combinations we get a model of a cortical column that will generate alpha type activity and for another set of parameters we get another model that will exhibit epileptic behaviour.22 Therefore, we consider this neural mass as a family of models, which we define as ˙ x(t) =fθ (x(t), u(t))

(14)

y(t) =Cx(t) + e(t),

(15)

where x(t) ∈ Rnx is a state vector representing the postsynaptic membrane potentials, generated by each population synapse, and their derivatives, nx is the number of states, u(t) represents the system input, which may be from afferent connections, other brain regions or electrical stimulation or other model inaccuracies. The function fθ (·) describes the dynamics, where θ ∈ Rnθ determines the model type and the behaviour it exhibits. The EEG is denoted by y(t), C is the observation matrix, and e(t) is the observation noise. The model focused on in the following estimation sections is the formulation by Jansen and Rit.11 Given space limitations we refer the reader to the article by Jansen and Rit11 where the original state-space equations are presented. The parameters that were described in this section are related to the Jansen and Rit paper by αpi αep αxp αpe = = = (16) A= 2e0 c1 2e0 c3 2e0 c2 2e0 αip B= (17) 2e0 c4 1 1 1 1 a= = = = (18) τpe τpi τep τxp 1 (19) b= , τip

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

8

where e0 is a parameter that scales the maximum firing rate, A and B are synaptic gains for excitation and inhibition respectively, and a and b are the reciprocals of the synaptic time constants for excitation and inhibition, respectively. By making the assumption that all excitatory synapses share the same time constants and by defining the connectivity constants, c1 , c2 , c3 and c4 , the network of neural masses on the left side of Figure 2 can be reduced to the system depicted on the right hand side.

input u(t)

input u(t)

spiny stellate population

ve (t)

vp (t)

ECoG pyramidal population

g(vi (t))

vi (t)

interneuron population

g(ve (t))

vp3 (t)

vp1 (t)

inhibitory feedback

excitatory feedback

ve (t) vp2 (t)

c2

+

∗Aa.he (t)

ECoG g(vp (t))

+

vp (t)

c1

g(ve (t))

∗Aa.he (t)

∗Bb.hi (t) g(vp (t)) c4

c3

g(vi (t)) vi (t)

Fig. 2. Model of a Cortical Column. The model shows three interconnected neural masses, which are pyramidal neurons, excitatory spiny stellate cells, and inhibitory interneurons. The specific subtype of neural population is defined by the parameters that describe the post-synaptic response kernels.

2.1. Synthetic Data The synthetic data used to test the estimation schemes employed parameters that generate alpha rhythm like activity.11 The model was discretized using Euler’s method with a time step of 1 ms. The input for the model was Gaussian noise with a standard deviation of 60 pulses per second and a mean of 90 or 220 pulses per second for the deterministic and stochastic estimators, respectively. The observation noise was also Gaussian with a mean of zero and a standard deviation of 20% of the size of the standard deviation of the simulated pyramidal membrane potential.

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

9

3. Deterministic Estimation: Adaptive Observer This section describes the adaptive observer23 and provides adequate details and equations for implementation. Proofs that demonstrate that the observer’s parameter and state estimates converge to the true values for the 23 Jansen and  Rit model are provided in Postoyan et al. In the following a a vector is denoted (a, b), for all a ∈ Rna , b ∈ Rnb . b The adaptive observer considered here has been designed for systems that have the state-space form, x˙ 0 = M0 x0 + φ0 (y)θ x˙ 1 = M1 x1 + φ1 (x0 , u)θ (20)

y = C1 x1 ,

where x0 ∈ Rn0 and x1 ∈ Rn1 are the states, θ ∈ Rp is a vector of constant and unknown parameters, y ∈ Rny is the measurement and u ∈ Rnu is the input. Note that we are assuming that the measurement noise is zero for the derivation. The Jansen-Rit model11 can be expressed in the state-space form of equation (20) by taking the states to be x0 = (x01 , x02 ) ∈ R2 and x1 = (x11 , . . . , x14 ) ∈ R4 , where x01 , x11 , x13 are membrane potential contributions of the pyramidal neurons and the excitatory and inhibitory interneurons respectively, and x02 , x12 , x14 are their respective derivatives. According to the notation of Jansen and Rit,11 x0 = (y0 , y3 ) and x1 = (y1 , y4 , y2 , y5 ). The measured output (ECoG) is y ∈ R, u ∈ R is the excitatory input from neighbouring columns, which is assumed to be known and θ = (A, B) ∈ Θ ⊆ R2 is a vector of unknown parameters, where A and B represent the synaptic gains of the excitatory and inhibitory neuronal populations, respectively. The matrices in (20) are defined as  C1 = 1 0 −1 0 , M0 = Ma , M1 = diag(Ma , Mb ),

 where Ma =

0 1 2 −a −2a



 , Mb =

0 1 −b2 −2b

 .

The inverse post-synaptic response time constant parameters, a, b ∈ R≥0

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

10

are assumed to be known. The nonlinear terms in (20) are given by   0 0 φ0 (y) = , a2e0 g(y) 0   0 0  ac2 2e0 g(c1 x01 ) + au  0 . φ1 (x0 , u) =    0 0 0 bc4 2e0 g(c3 x01 ) The connectivity parameters, c1 , c2 , c3 , c4 ∈ R≥0 are assumed to be known, g denotes the sigmoid function defined in equation (13) and 2e0 ∈ R≥0 defines the maximal firing rate. It is worthwhile to note that vp , ve and vi in the right side of Figure 2 correspond to y = x11 − x13 , c1 x01 , and c3 x01 , respectively. For ease of notation, we write (20) in the following form x˙ = Mx + φ(y, u, x)θ y = Cx,

(21)

where x = (x0 , x1 ), M = diag(M0 , M1 ), C = (0, C1 ) and φ = (φ0 , φ1 ). The nonlinear terms, φ0 : R → Rn0 × Rp and φ1 : Rn0 × Rnu → Rn1 × Rp , are globally Lipschitz and bounded. We consider the following adaptive observer for system (21): ˆ˙ x yˆ ˆθ˙ ˙ Υ ˙ Q

ˆ + Γ(y − yˆ) ˆ )θ = Mˆ x + φ(y, u, x = Cˆ x ¯ = Γ(y − yˆ) ˆ ), = MΥ + ∆φ(y, u, x = dQ − dQΥ> C> CΥQ,

(22) with Υ(0) = 0 with Q(0) = Q> (0) > 0,

¯ Γ ¯ = QΥ> C> and ∆ = diag(In , 1 In ). The quantity where Γ = ∆−1 ΥΓ, 0 d 1 d > 0 is a design parameter that trades off convergence speed for accuracy. ˆ the estimate of θ. The ˆ , denotes the estimate of x and θ The vector, x (n0 +n1 )×p variable Υ ∈ R is initialized to Υ(0) = 0. An essential assumption required for our adaptive observer to work is as follows. For any signals u, y, x ˆ that belong to L∞ , there exist a1 , a2 ∈ R≥0 , T ∈ R≥0 such that the solution to ˙ = MΥ + ∆φ(y, u, x ˆ ) with Υ(0) = 0, Υ

(23)

satisfies, for all t ≥ 0, a1 I2 ≤

R t+T t

Υ> (τ )C> CΥ(τ )dτ ≤ a2 I2 .

(24)

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

11

This condition (24) is known as the persistency of excitation of the signal CΥ(t) and is a well-known condition in the control theory literature.24 Inequality (24) is hard to verify in general, however, simulations demonstrate that this condition is satisfied for the Jansen and Rit model. 3.1. Demonstration of Adaptive Observer using ECoG Data Simulation Here we demonstrate online parameter estimation for two parameters of the Jansen and Rit model, A and B, that represent the synaptic gains of the excitatory and inhibitory neuronal populations, respectively. The design parameter of the adaptive observer was set to d = 1.5. In Figure 3 it can be seen that the parameter estimates converge to the true estimates within 10 seconds, even though the initial estimates begin at either 60% above or 60% below the actual value in the two different simulations. 4. Stochastic Estimation In this section, we shall switch to discrete time notation, indicated by using the superscript t to index the samples. Furthermore, we shall use the same notation f to describe the system, bearing in mind that the discrete version has different properties to the continuous version. An alternative way of describing the model in equation 14 is to allow the parameters to be time varying, but on a much slower time scale than the states that have previously been defined. By making the assumption that the parameters are varying slowly, we can form the augmented state-space model  xt+1 =f xt , θ t + εt (25) θ t+1 =0 + ϕt t

t

(26) t

y =Cx + e ,

(27)

where the new disturbance term, ϕt , which is an i.i.d. Gaussian, essentially adds a small amount of uncertainty to the trivial dynamics of the parameters. By augmenting the state vector with parameters that are to be estimated, we can perform state and parameter estimation simultaneously. This allows the Kalman filter to provide small corrections to the model predictions with each arrival of new data.25 Our goal is to estimate the states, xt , and the slowly varying parameters, θ t , given knowledge of the biophysics of the mass action of the brain and the noisy ECoG mea-

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

12 vp(t)

y(t)

Est. From Above

Est. from Below

Actual

y(t) (mV)

a) 2 1 0 10

10.5 Time (s)

11

A (mV)

b) 4

2

B (mV)

c) 30 20 10 0

5

10

15

Time (s)

Fig. 3. Adaptive Observer Estimation Results from Synthetic ECoG Data. a) An example of a 1 second segment of the synthetic ECoG data that was used, illustrating the level of noise added to the model output. Black line: the membrane potential of the pyramidal population, y. Gray dots: the ECoG signal that was used for estimation. The estimates of the (b) excitatory gain parameter and the (c) inhibitory gain parameter from when the initial estimate was 60% above and below the actual value, respectively. Dotted line: the actual parameter value. Dot-dash line: the estimation results when the estimated value was initilalized at 60% above. Dashed line: the results when initialized at 60% below.

surements. Therefore, we shall define the augmented state vector as  > xtθ = xt θ t .

(28)

Our goal is to find  t 1 2  t ˆ t+ x θ = E xθ |y , y , . . . , y ,

(29)

which is known as the a posteriori augmented state estimate. This can be found by optimally combining the a priori augmented state estimate, which is defined as  t 1 2  t−1 ˆ t− x , (30) θ = E xθ |y , y , . . . , y with the noisy measurements. The a priori estimate is an estimate of the same quantity as the a posteriori estimate, but with only using the prediction from the model and not using the information from the most recent

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

13

observation at time t.26 The way the a priori estimate of the augmented state is corrected using the noisy observation is based on the level of confidence we have in our measurement (the noise level) and the inaccuracy of our model (the disturbance level). This is the basic philosophy of the Kalman filter. We can define the a posteriori and a priori augmented state estimate error covariances as   t > ˆ t− =E (xt − x ˆ t− ˆ t− P (31) θ θ )(xθ − x θ )  t  t+ t+ t+ t > ˆ =E (x − x ˆ )(x − x ˆ ) , P (32) θ

θ

θ

θ

respectively. For linear systems, the famous Kalman filter provides the solution to the estimation problem. However, for the nonlinear neural mass model, the integration that is required for solving the expectations has no closedform solution. Therefore, an approximate solution is required for efficient filtering (of high dimensional systems). An appropriate method for this approximation is the unscented transform, which leads to the unscented Kalman filter.27 The unscented transform (UT) is a method for approximating the statistics of a random variable that undergoes a nonlinear transformation. Consider transforming an n-dimensional random variable, xθ (augmented state vector), through a nonlinear function, xtθ = f xt−1 . Assume θ t−1 t−1 t−1 xθ ∼ N (¯ xθ , Pxθ ). To map the statistics through the nonlinear transform, we first define the so-called (2nxθ + 1) sigma vectors, which form the sigma matrix X , as X0t−1 =¯ xt−1 θ Xit−1 Xit−1

(33)

 q t−1 =¯ xt−1 + (n + λ) P , x xθ θ θ i  q t−1 (n + λ) P =¯ xt−1 − xθ xθ θ

i = 1, . . . , nxθ ,

i = nxθ + 1, . . . , 2nxθ ,

(34) (35)

i−nxθ

where the quantity λ = γ 2 (nxθ + κ) − nxθ (usually a small positive value e.g. 1 ≥ γ ≥ 10−4 ) is a scaling parameter γ. The constant γ determines ¯ t . The constant, κ, the spread of the sigma points around the mean, x is a secondary scaling  parameter, which is usually set to 0 or 3 − nxθ .  q t−1 (nxθ + λ) Pxθ is the ith column of the matrix square root (e.g. lower i

triangle Cholesky factorisation).

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

14

The sigma vectors are propagated through the nonlinear system,  Xit = f Xit−1 i = 1, . . . , 2nxθ , (36) and the mean and covariance for the transformed variable, which link back to equations 30 and 31, are approximated by 2nxθ

¯ tθ x

=ˆ xt− θ

X



(m)

Wi

Xit

(37)

i=0 2nxθ

ˆ t− ≈ =P xθ

Ptxθ

X

(c)

Wi

ˆ t− Xit − x θ



ˆ t− Xit − x θ

>

(38)

i=0

where the weights, Wi , are (m)

W0

(c)

W0 (m)

Wi

(c)

= Wi

λ nxθ + λ  λ + 1 − γ2 + β = nxθ + λ 1 = i = 1, . . . , 2nxθ , 2 (nxθ + λ)

=

(39) (40) (41)

and β is a variable that includes prior knowledge of the distribution of xθ . By using the UT to propagate the state and parameter estimates and errors through time, the standard Kalman filter update equations (since the observation equation is modeled as being linear) can be used. The UT captures the properties of the nonlinear transform of a Gaussian random variable to 3rd order Taylor series expansion for all nonlinearities (and higher with the appropriate choice of γ and β 28 ). 4.1. Demonstration of UKF using ECoG Data Simulation In this section, we demonstrate how the stochastic estimation framework is able to estimate the neural mass model parameters when their values are initialized 60% above and below their actual values. The results are shown in Figure 4 along with a 1 second of the data that was used. The results show that the parameter values converge to the actual values for both initial condition cases. The discrepancy between the time it takes to converge when initialised from above compared to below is similar to the deterministic case. This effect may be due to the saturating effects of the sigmoidal nonlinearity. Alternatively, it may be due to the smaller time constants (larger reciprocals) allowing the system to respond more readily to changes. It indicates that is is advantageous to initialize parameters to values in the upper region of the physiological range.

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

15 Est. From Above b) a (1/s)

y(t)

a) 20 y(t) (mV)

vp(t)

10

100 50

50

50.5 d)

5 4

b (1/s)

A (mV)

Actual

150

0 49.5

c)

Est. from Below

3

80 60 40

2 20 f)

e)

v0 (mV)

B (mV)

30 20

8 6 4

10 0

50 Time (s)

100

0

50 Time (s)

100

Fig. 4. UKF Estimation Results from Synthetic Data. a) An example of the synthetic data that was used, illustrating the level of noise added to the model output. Black line: the membrane potential of the pyramidal population. Gray dots: the ECoG signal that was used for estimation. Subplots b) to f) show the estimated parameters from when the initial guess was 60% above and below the actual value. Dotted line: the actual parameter value. Dot-dash line: the estimation results when the estimated value was initilalized at 60% above. Dashed line: the results when initialized at 60% below. b) Estimate of the reciprocal of the excitatory time constant. c) Estimate of the excitatory gain parameter. d) Estimate of the reciprocal of the inhibitory synaptic time constant. e) Estimate of the inhibitory gain parameter. f) Estimate of the sigmoid threshold parameter.

4.2. Demonstration of using Real ECoG Data 4.2.1. ECoG Data The data that was used in this section was collected from a patient undergoing assessment for resective surgery at St. Vincent’s Hospital Melbourne, and is used with patient consent and permission from the Human Research Ethics Committee (protocol HRECA-006/08). Data was recorded from a sub-dural grid electrode with electrode spacing of 0.5 cm and surface diameter of 2 mm (Ad-Tech Medical). The most focal differential channel was used, relative to the seizure onset zone. The data was sampled at a rate of 1kHz and low-pass filtered with a cut-off of 95 Hz. The data was re-scaled by a factor of 1 × 104 to have an amplitude in the same order of magnitude

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

16

to the output of the neural mass model. A 5 minute epoch of the data was used, with approximately 200 seconds before electrographic seizure onset and 100 seconds post-onset. 4.2.2. Parameter Estimates A major challenge in applying the method lies in initializing the estimator. The initial parameter estimates were chosen based on a trial and error approach. The Jansen and Rit parameters to simulate the alpha rhythm were initially used. If any one parameter diverged rapidly then the initial condition was adjusted in the direction of divergence in small steps until the direction changed. After some tuning, the parameters gave a steady value until approximately 40 seconds prior to seizure onset, where the inhibitory population synaptic time constant began to diverge from the base-line level. Figure 5 shows the results from the ECoG data. As expected, the gains for both excitatory and inhibitory, A and B, populations increased during the seizure. Although it is not shown explicitly, the ratio of excitation to inhibition was fairly constant prior to the seizure, increased during the seizure, and returned to pre-seizure levels after the seizure. The excitatory synaptic time constant became shorter during the seizure, which explains higher frequency activity during the event. Interestingly, the inhibitory time constant decreased approximately 40 seconds prior to the onset, peaked shortly after the onset, and then decreased during the seizure. This provides a new intuition into higher frequency oscillations prior to seizures and the chirping effect observed during seizures. Finally, the threshold parameter decreased at seizure onset and then increased at seizure termination. We speculate that this may be a result in a shift in the resting membrane potential due seizure related changes in ionic concentrations; however, this remains to be tested. 5. Discussion We have presented a novel method developed by our group, the adaptive observer, and compared it to a novel application of an existing method, the augmented state UKF, for state and parameter estimation of ECoG data using a neural mass model. These methods can be used for inferring underlying brain changes based on a single ECoG recording. The stochastic estimation method relies on an extension of the UKF, while the deterministic estimation method relies on Lyapunov and persistency of excitation-based approaches for adaptive observer design. The adaptive observer estimates

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

17

a)

b) 120 a (1/s)

ECoG (mV)

10 0 −10

100 80 60

c)

d)

3

b (1/s)

A (mV)

4

2 1

150 100

e)

f) 25 20

3 v0 (mV)

B (mV)

200

15 10 5

2 1 0 −1

100

200 Time (s)

300

100

200

300

Time (s)

Fig. 5. UKF Estimation Results from ECoG Data. a) The ECoG data from a focal channel recorded from a grid array. b) Estimate of the reciprocal of the excitatory synaptic time constant. c) Estimate of the excitatory gain parameter. d) Estimate of the reciprocal of the inhibitory time constant. e) Estimate of the inhibitory gain parameter. f) Estimate of the sigmoid threshold parameter.

are guaranteed to converge to the true state and parameter values of the neural mass model provided the assumption of persistency of excitation holds. Although convergence is not guaranteed for the stochastic method, the UKF deals well with uncertainty and model error because it takes these into account implicitly. Due to the Lyapunov-based design, the adaptive observer has a degree of robustness to uncertainty, and nominal robustness for small model-system mismatch. However, the adaptive observer needs to be analysed further in terms of robustness to uncertainties. Here we have shown in simulations that the adaptive observer is robust to a large degree of measurement noise although this is not considered in the observer design. The adaptive observer was applied to the real ECoG data, although it was not shown in the chapter, and the estimation results were unstable. This is may be due to the need for the adaptive observer to know the input to the cortical column. Further work is needed to specifically design an adaptive observer that is robust to input uncertainty. The flexibility of the UKF to uncertainties means that it performs well

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

18

at estimating a larger number of parameters than the adaptive observer, and that it works effectively on real data. This result of online estimation of neural mass model parameters applied to real seizure ECoG data is the first of its kind. Moreover, the parameter changes occurring at seizure onset, such as an increase in postsynaptic potential amplitudes (A,B) and an initial decrease in threshold v0 , appear consistent with what might be expected during seizures. The online estimation presented here contrasts with and compliments the excellent offline estimation techniques that have also been applied to epilepsy data.29,30 The true meaning and interpretability of parameter changes requires further study. In particular, different parameter combinations may produce similar ECoG seizure dynamics and the question needs to be answered: which parameter combination is the true combination? One way to get around this problem is to limit the complexity of the computational models to allow for reliable estimation and at the same time reduce the dimensionality of the parameter space. For more complex models, higher fidelity data acquisition and greater signal sampling is required for coping with uncertainties. To put it another way, the models can be considered tools for interpolating missing or hidden data. To capture more complex features of the system, richer data is required or the interpolation will be overly simplified. Signal scaling, re-referencing, offsets, and DC-filtering also need to be considered in future estimator designs. Such biologically motivated online estimator designs will likely revolutionise how we can image the underlying activity of the brain and devise improved methods for neurological monitoring, control and treatment. 6. Acknowledgements This research was partly supported by the Australian Research Council (Linkage Project LP100200571). The Bionics Institute acknowledges the support it receives from the Victorian State Government through the Operational Infrastructure Support Program. References 1. S. Schiff, Neural Control Engineering: The Emerging Intersection Between Control Theory and Neuroscience (The MIT Press, Cambridge, MA, 2011). 2. E. Marder and J. Goaillard, Nature Reviews Neuroscience 7, 563 (2006). 3. F. Marten, S. Rodrigues, P. Suffczynski, M. Richardson and J. Terry, Physical Review E 79, p. 21911 (2009). 4. F. Mormann, R. Andrzejak, C. Elger and K. Lehnertz, Brain 130, 314 (2007).

November 24, 2012

2:7

WSPC - Proceedings Trim Size: 9in x 6in

ws-procs9x6

19

5. L. Kuhlmann, D. Freestone, A. Lai, A. Burkitt, K. Fuller, D. Grayden, L. Seiderer, S. Vogrin, I. Mareels and M. Cook, Epilepsy Research 91, 214 (2010). 6. H. Khalil, Nonlinear Systems, 3rd edn. (Prentice-Hall, Englewood Cliffs, New Jersey, U.S.A., 2002). 7. G. Bastin and M. Gevers, IEEE Transactions on Automatic Control 33, 650 (1988). 8. A. Zemouche and M. Boutayeb, Systems and Control Letters 58, 282 (2009). 9. O. David, S. Kiebel, L. Harrison, J. Mattout, J. Kilner and K. Friston, NeuroImage 30, 1255 (2006). 10. D. Freestone, P. Aram, M. Dewar, K. Scerri, D. Grayden and V. Kadirkamanathan, NeuroImage 56, 1043 (2011). 11. B. Jansen and V. Rit, Biological Cybernetics 73, 357 (1995). 12. Q. Zhang, Automatic Control, IEEE Transactions on 47, 525 (2002). 13. G. Besan¸con, Systems & control letters 41, 271 (2000). 14. S. Julier and J. Uhlmann, A new extension of the Kalman filter to nonlinear systems, in Int. Symp. Aerospace/Defense Sensing, Simul. and Controls, Orlando, FL, 1997. unscented kalman filter. 15. R. van der Merwe (2003). 16. O. Faugeras, J. Touboul and B. Cessac, Frontiers in computational neuroscience 3 (2009). 17. F. Freyer, J. Roberts, R. Becker, P. Robinson, P. Ritter and M. Breakspear, The Journal of Neuroscience 31, 6353 (2011). 18. G. Deco, V. Jirsa, P. Robinson, M. Breakspear and K. Friston, PLoS Comput Biol 4, p. e1000092 (2008). 19. F. L. da Silva, A. Hoek, H. Smith and L. Zetterberg, Cybernetic 15, 27 (1974). 20. F. Wendling, A. Hernandez, J. Bellanger, P. Chauvel and F. Bartolomei, J. Clin. Neurophysiol. 22, 343 (2005). 21. O. David and K. Friston, NeuroImage 20, 1743 (2003). 22. F. Wendling, F. Bartolomei, J. Bellanger and P. Chauvel, European Journal of Neuroscience 15, 1499 (2002). 23. R. Postoyan, M. Chong, D. Neˇsi´c and L. Kuhlmann, Parameter and state estimation for a class of neural mass models, in IEEE Conference on Decision and Control , (Maui, Hawaii, 2012). 24. P. Ioannou and J. Sun, Robust Adaptive Control (Prentice-Hall, New Jersey, 1996). 25. E. W. Ericwan, R. V. D. Merwe and A. T. Nelson, 666 (2000). 26. R. Kalman et al., Journal of basic Engineering 82, 35 (1960). 27. S. Haykin, Kalman filtering and neural networks (Wiley Online Library, 2001). 28. S. Julier, The scaled unscented transformation, in American Control Conference, 2002. Proceedings of the 2002 , 2002. 29. A. Blenkinsop, A. Valentin, M. Richardson and J. Terry, Eur J Neurosci. 36, 2188 (2012). 30. F. Wendling, F. Bartolomei, F. Mina, C. Huneau and P. Benquet, Eur J Neurosci. 36, 2164 (2012).