Inferring effective connectivity in epilepsy using dynamic causal modeling

Inferring effective connectivity in epilepsy using dynamic causal modeling

JID:IRBM AID:378 /FLA [m5+; v1.211; Prn:9/10/2015; 8:57] P.1 (1-10) Disponible en ligne sur ScienceDirect www.sciencedirect.com IRBM ••• (••••) •••...

1008KB Sizes 0 Downloads 89 Views

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.1 (1-10)

Disponible en ligne sur

ScienceDirect www.sciencedirect.com IRBM ••• (••••) •••–•••

Original article

Inferring effective connectivity in epilepsy using dynamic causal modeling W. Xiang a,b,e , C. Yang a,b,e , J.-J. Bellanger c,d,e , H. Shu a,b,e,∗ , R. Le Bouquin Jeannès c,d,e a LIST, School of Computer Science and Engineering, Southeast University, Nanjing, 210096, China b Key Laboratory of Computer Network and Information Integration, Southeast University, Nanjing, 210096, China c INSERM, U1099, Rennes, F-35000, France d LTSI, Université de Rennes 1, Rennes, 35042, France e Centre de Recherche en Information Biomédicale Sino-Français (CRIBs), France

Received 23 January 2015; received in revised form 11 September 2015; accepted 11 September 2015

Abstract This study deals with effective connectivity analysis among distant neural ensembles recorded with intracerebral near field electrodes during seizures in the brain of epileptic patients. Our goal is to analyze the ability of Dynamic Causal Modeling (DCM) approach to detect causal links when the underlying model is a well-known neural population model dedicated to the simulation of epileptic activities in hippocampus. From the state-space description of the system obtained by coupling a pair of such models, a linearization around the equilibrium state leads to a transition matrix and a parametrized description of the power spectral density matrix for the corresponding pair of output ElectroEncephaloGraphic (EEG) signals in the two-population model. Then, the parameters of this global model are estimated in a Bayesian framework from simulated EEG signals by the Expectation Maximization (EM) algorithm, and Log Bayes Factors are employed to discriminate among the possible effective connectivity hypotheses. Simulation results show that DCM can identify and distinguish the independence, unidirectional or bidirectional interactions between two epileptic populations. © 2015 AGBM. Published by Elsevier Masson SAS. All rights reserved. Keywords: Epilepsy; Effective connectivity; Dynamic causal modeling; Physiological model; Power spectral density; EM algorithm

1. Introduction In recent years, due to the high incidence of epilepsy, its treatment has become a hot research topic. About 30% patients are drug-resistant and this pathology affects all their life [1]. Drug-resistant epilepsies are often partial with an epileptogenic zone (EZ) located in a particular brain region. In this case, in order to remove the EZ to avoid seizures, a surgical treatment should be considered. In order to investigate the organization of the epileptogenic zone and its surrounding perturbed cerebral sites, a good characterization of the causal relationships between different intracerebral electroencephalographic (iEEG) brain recordings using depth electrodes can help in improving chirurgical planning prior to surgery [2]. In this perspective, * Corresponding author.

E-mail address: [email protected] (H. Shu). http://dx.doi.org/10.1016/j.irbm.2015.09.001 1959-0318/© 2015 AGBM. Published by Elsevier Masson SAS. All rights reserved.

a number of methods aiming at evaluating statistical causal or unidirectional relationships between EEG signals have been already proposed. They are classified into two categories, modelfree based approaches which were firstly developed and modelbased ones. These techniques fall under the scope of effective connectivity [3], which denotes some causal effect of one neural system over another one in neuroscience. As it seems a priori useful, physiological hypotheses have been taken into account in Dynamic Causal Modeling (DCM) methodology first introduced by K. Friston to analyze effective connectivity in fMRI responses in cognitive tasks [4] and for which a number of studies can be found in the literature [5–8]. Subsequently, this technique has been extended to characterize the distributed neuronal architectures in EEG signals underlying both event related potentials in time domain [9,10] and cross spectral density in frequency domain [11,12]. DCM includes a biophysical realistic model to generate brain data and statistical inference from data, to study how activities in distinct brain areas affect

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.2 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

2

each other. From a practical point of view, to select a particular method it is necessary to consider a compromise between sensitivity, robustness and process execution time. As far as epilepsy is concerned, in the first kind of approaches to infer connectivity, Wiener–Granger causality index [13], coherence based approaches [14], nonlinear correlation index [15] and transfer entropy [16] are the most representative tools. In this paper we present a first investigation of a DCM approach involving neural populations models taking physiological evidences into account. In this work, we are only concerned with the characterization of steady state effective connectivity between two different potentially epileptic populations corresponding to ictal onset where a narrowband tonic activity is observed. In this contribution, the two following aspects are considered: (1) the inputs of the physiology-based model [17] used to generate the epileptic signals are neither controlled nor measured, and we assume that these inputs are caused by endogenous fluctuations in cortical activity considered as white noise; (2) at the seizure onset, the epileptic populations present a fast narrow band activity around 25 Hz. Therefore, this paper is based on DCM for Steady State Response (SSR), which is designed to fit spectral data features [11,12]. This paper comprises five sections: in Section 2, we consider the physiology-based model to specify the cross spectral density of epileptic signals. The Variational Bayesian Expectation Maximization (VBEM) algorithm [18] is considered to estimate the model parameters and then the free energy is used to get the model evidence in Section 3. In Section 4, in order to address the validity of DCM-SSR, we compute theoretical autoand cross-spectra of two populations from a given set of chosen population model parameters and coupling parameters (imposing the type of connectivity) which provide a “ground truth” for our simulations. Finally, concluding remarks and some perspectives are given in Section 5. 2. Methods and materials 2.1. Physiology-based model Dynamic Causal Modeling [4] considers distributed brain networks as a connected set of neuronal populations, each population being described by a set of differential equations. For the hippocampal model [17] considered here, given two populations j and k, the corresponding equations for population j can be written by the following set of 16 equations:

(10)

x˙i,j (t) = xi+5,j (t) i = 1, 2, . . . , 5   x˙6,j (t) = Aj aS x13,j (t) − 2ax6,j (t) − a 2 x1,j (t)   x˙7,j (t) = Aj aC2 S x14,j (t) − 2ax7,j (t) − a 2 x2,j (t)   + Aj a uj (t) + Kj k x11,k (t)   x˙8,j (t) = Bj bC4 S x15,j (t) − 2bx8,j (t) − b2 x3,j (t)   x˙9,j (t) = Gj gC7 S x16,j (t) − 2gx9,j (t) − g 2 x4,j (t)   x˙10,j (t) = Bj bS x15,j (t) − 2bx10,j (t) − b2 x5,j (t)

(11)

x˙11,j (t) = x12,j (t)

(1..5) (6) (7) (8) (9)

Fig. 1. Interactions between different neuronal subpopulations of physiologybased model [17,19].

(12)

  x˙12,j (t) = Ad,j ad S x13,j (t) − 2ad x12,j (t) − ad2 x11,j (t)

(13)

x˙13,j (t) = x7,j (t) − x8,j (t) − x9,j (t)

(14)

x˙14,j (t) = C1 x6,j (t)

(15)

x˙15,j (t) = C3 x6,j (t)

(16)

x˙16,j (t) = C5 x6,j (t) − C6 x10,j (t)

(1)

where S(·) is a nonlinear sigmoidal function, which transforms a presomatic time varying membrane potential x(t) to an output firing rate:   S x(t) =

2e0 1 + exp(−r(V0 − x(t)))

(2)

and e0 , r and V0 are constants. As shown in Fig. 1, each population is composed of 3 subpopulations, two of them being inhibitory interneurons (with fast and slow inhibitions) and the third one corresponding to excitatory interneurons with intrinsic connections, whose strengths are parameterized by C = {C1 , C2 , C3 , C4 , C5 , C6 , C7 }. The extrinsic coupling parameter Kj k represents the strength of the effective directional connection from population k to population j (Fig. 2). The input uj (t) summarizes the influence of distant afferent neurons and is assumed to be a stationary white Gaussian noise with constant mean mj and Power Spectral Density (PSD) equal to a positive constant αj . From empirical adjustments [17,20] a common choice for the values of mj and mk is 90, and 900 for αj and αk . This means that the exogenous influences uj (t) and uk (t) are statistically identical and assumed to be independent. A less constrained choice which relaxes the constraint αj = αk = 900 will be considered in the experimental section. The ratios 1/a, 1/b, 1/g and 1/ad correspond to the time constants in synaptic kinetics (respectively equal to 1/100, 1/30, 1/250 and 1/33). To focus on the parameters of interest, only A, Ad (equal to A in our model), B, G, K and α are expected to vary during a transition from a normal process to an epileptic seizure and have to be estimated. However, these parameters are assumed to be constant on the observation interval obtained beforehand from a segmentation procedure (automatic or interactive). The remaining parameters

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.3 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

3

Fig. 2. Extrinsic connections between two populations j and k. Table 1 Priors on free biophysical parameters [20] in one population. Centered Gaussian Scale parameters θ = [θA , θB , θG , θK , θα ]T are introduced (θi(i=A,B,G,K,α) ∼ N (μθi , −1 θi )). These parameters are logtransformed to enforce positivity. Parameter

Description

Prior distribution

A B G K α

Synaptic gain for excitation Dendritic slow inhibition Somatic fast inhibition Extrinsic connections Input white noise PSD (optional)

A = 5 exp(θA ) θA ∼ N (0, 0.5) B = 3 exp(θB ) θB ∼ N (0, 0.5) G = 20 exp(θG ) θG ∼ N (0, 0.5) K = 500 exp(θK ) θK ∼ N (0, 1) α = 5 exp(θα ) θα ∼ N (0, 0.25)

are supposed to be fixed and known as in [17,19]. The a priori density probabilities (priors) introduced for the Bayesian inference correspond to log-normal distribution as reported in Table 1. The observed EEG signal yj (t) is supposed to be proportional to the component x13,j (t) after cancellation of its temporal average x13,j :   (3) yj (t) = l x13,j (t) − x13,j where l is the instrumentation gain (we set l = 1). The mean value is subtracted as the instrumentation preamplifier frequency response cancels at the null frequency. The output can be interpreted as intracranial electroencephalographic signals similar to those recorded with proximal electrodes in the brain (local field potentials). 2.2. Model based cross power spectral density Aiming at modeling two populations j and k, the differential equation system and the observation equation respectively described by (1) and (3) must be duplicated and can be regrouped in a global differential system including 2 × 16 one-order equations and a two-component observation vector:   ˙ = f x(t), θ + Du(t) x(t)   (4) y(t) = L x(t) − x

where f (x(t), θ ) = [fj (x(t), θj )T , fk (x(t), θk )T ]T is the nonlinear global state derivative function and θ = [θjT , θkT ]T is the global vector of parameters we want to estimate. The state vector x(t) = [x1:16,j (t), x1:16,k (t)]T ∈ 32×1 includes all state components implied in the pair of populations, x is the time average of x(t). The input two-dimensional vector u(t) = [uj (t), uk (t)]T ∈ 2×1 includes two independent white Gaussian noises whose means and standard deviations are typically equal to 90 and 30 respectively. These moment values were firstly adopted in accordance with [17] (where they were chosen empirically from simulations and from real observations). The input matrix D = [Dj , Dk ] ∈ 32×2 contains only two non-zero entries, one in line (7) (see (1)) for the population j and another in line (7) + (16) for the population k (see (1)) (in the global differential system which represents the populations pair). The observed two-dimensional output y(t) = [yj (t), yk (t)]T ∈ 2×1 is obtained from the product of the gain matrix L = [Lj , Lk ] ∈ 2×32 with the state vector (after cancellation of the time average). This matrix also contains only two non-zero entries, one in column (13) (see (3)) for population j and another in column (13) + (16) (see (3)) for population k. In order to simplify the expression of the PSD matrix of the bi-dimensional random process y (which results from a nonlinear operator acting on the random process u) we linearize the differential system around the equilibrium state x 0 (θ, mj , mk ) defined as the solution in 32 , when it exists,

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.4 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

4

of 0 = f (x 0 , θ ) + D[mj , mk ]T , and compute the PSD of the perturbation δx(t)  x(t) − x 0 resulting from the action of the centered version of u(t), δu(t)  u(t) − [mj , mk ]T on the linearized system. Note that x 0  x. Clearly, if mj and mk are free parameters, they influence not only the value of the equilibrium point x 0 (θ, mj , mk ) when it exists, but also its existence. Now, in the present study, we comply with the two following constraints 1) mj = mk = m and 2) m = 90 to reduce the estimation task of the algorithm introduced in Section 3, and because we considered that the mean levels of the surrounding background stimulations could be approximately equal for both populations and m = 90 is the standard choice as indicated above. We do not link the parameters mj , mk , αj , and αk to the parameters Kj k and Kkj since we consider that the modifications of the interaction levels between the two populations principally come from hyper synchronization phenomena inside one (or the two) population(s) without modifying notably the level of the global background stimulation. As in [21], we consider the following linearization of (4):   d δx (t) = (x 0 )δx(t) + Dδu(t) dt y(t) = Lδx(t) (5) (x(t),θ) where (x 0 ) = ∂f ∂x(t) |x(t)=x 0 ∈ 32×32 is the Jacobian matrix of the f (·, θ ) computed at x = x 0 . The Laplace transform of (5) allows to write:

sδX(s) = (x 0 )δX(s) + DδU (s) Y (s) = LδX(s)  δX(s) = (sI − )−1 DδU (s) ⇒ Y (s, θ ) = L(sI − )−1 DδU (s)

(6)

which leads, after identification to the 2 × 2 linear Laplace transfer function matrix:   Hjj (s, θ ), Hj k (s, θ ) H (s, θ ) = Hkj (s, θ ), Hkk (s, θ ) −1

D∈ . (7) √ When evaluated at s = −1ω, H gives the frequency transfer function. So, the power spectral density matrix for the pair of signals (yj , yk ) is given by:   gjj (ω, θ )gj k (ω, θ ) g(ω, θ ) = gkj (ω, θ )gkk (ω, θ )   √ √ αj 0 (8) g(ω, θ ) = H ( −1ω, θ ) H ∗ ( −1ω, θ ) 0 αk = L(sI − )

2×2

where superscript ∗ denotes the conjugate transpose operator. To detail one of the entries of g(ω, θ ), the cross spectral density (CSD) is: √ √ gj k (ω, θ ) = Hjj ( −1ω, θ )αj Hj k ( −1ω, θ )∗ √ √ + Hj k ( −1ω, θ )αk Hkk ( −1ω, θ )∗ . (9) Note that, if we assume the statistical stationarity of the solution of (5) as a prerequisite for introducing the PSD of the

output y, we must also assume that none of the (x 0 ) eigenvalues has a non-negative real part. Here, from a physiological point of view, this hypothesis holds since we consider that the observed neural populations do not behave together as unstable systems able to produce spontaneous oscillations without any external excitation from afferent neurons. Now, as indicated below, some steps in the parameter estimation algorithm described later in Section 3 can lead to values of (θ, mj , mk ) getting out of the stability domain. In this case, the steady state trajectory of the nonlinear differential system ˙ = f (x(t), θ ) + D[mj , mk ]T generally stays on a limit cyx(t) ˙ = f (x(t), θ ) + Du(t) is cle and the stationary solution of x(t) a randomly perturbed version of this cycle. This corresponds to an aperiodic output y which admits a diffuse narrow band PSD without pure harmonic components and which may look like the narrow band PSD of real epileptic activities such as those analyzed here. Now, to compute a PSD function from (8) the transfer function H (s) = L(sI − )−1 D must be stable, i.e. the eigenvalues of the matrix  (see (6)) must be strictly included in the left part of the s-plane. So, if there is no equilibrium state x 0 (θ, mj , mk ), it is impossible to determine the PSD matrix components following the linearization approach. This problem can be bypassed considering the two following steps: Step 1: Choose x 0 = x where x is an estimation of the asymptotic temporal average of the steady state trajectory for the system given in (4) and compute (x 0 ) = (x) (according to the fact that this trajectory is not far from x 0 ). Step 2: Compute the eigenvalues of (x) and eventually stabilize (5) if some eigenvalues have non-negative real parts (in this case, they are set to small negative values). Now, considering that we infer parameters from a measured PSD matrix g(ω) ˜ (an estimation from the observed temporal data), we introduce the likelihood function [11]: g(ω ˜ l ) = g(ωl , θ) + ε(ωl ),   εjj (ω)εj k (ω) ε(ω) = εkj (ω)εkk (ω)

l = 1, .., N (10)

To get g(ω), ˜ we used a 12-pole AR modeling as in the SPM Spectral Toolbox [11]. Here, for any (n, m) ∈ (j, k) × (j, k), εnm is modeled as a realization of a complex i.i.d. (independent and identically distributed) Gaussian sequence nm (ωl ), l = 1, .., N such that:     Re nm (ωl ) ∼ N 0, nm (λ)     Im nm (ωl ) ∼ N 0, nm (λ)

nm (λ) = exp(λnm )V (ωl )

(11)

where λmn is an unknown hyperparameter, V is a constant and the observation errors jj (ωl ), kk (ωl ), j k (ωl

) are assumed to be independent for any (l, l , l

). (5)–(11) specify the PSD between two populations given by parameters θ , so the power spectral density matrix is an analytic function of the parameters θ . The likelihood function p(g(ω)|θ ˜ ) can be obtained as above and the prior density p(θ) corresponds to a normal distribution N (μθ , −1 θ ) (Table 1), the joint density for the measured PSD matrix and the parameters

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.5 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

vector is p(g(ω), ˜ θ ) = p(g(ω)|θ ˜ )p(θ ). Hence, a full probabilistic generative model for PSD data is specified in terms of biophysical parameters. 3. Model inversion In DCM, Bayesian inversion provides parameter estimation and allows model comparison for competing models’ hypotheses [9]. Concerning parameter estimation, for a given measured power spectral density g(ω) ˜ the posterior parameter distribution p(θ |g(ω), ˜ m) can be theoretically computed for any model m defined by its network architecture (different directions of information flow between the two populations can be considered by constraining one of the parameters Kj k or Kkj , or both of them, to be equal to zero, so leading to 4 distinct cases). Now, obtaining the theoretical expression of p(θ |g(ω), ˜ m) is too cumbersome, and the well-known Monte Carlo Markov chain methodology (which samples θ values) requests a large transient period before convergence, i.e. before the distance be˜ m) tween the Markov chain state distribution and p(θ |g(ω), becomes sufficiently small. So, we used the standard variational Bayesian scheme [18] which searches the best approximation of p(θ |g(ω), ˜ m) in a given class of posterior θ distributions constrained to be normal distributions with independent θ components. More precisely this method consists in maximizing iteratively the free energy Fm :      Fm = ln p g(ω)|m ˜ − KL q(θ), p θ |g(ω), ˜ m

(12)

under the constraint that the computed approximate posterior density q(θ) corresponds to a normal distribution N (v, v ), KL denoting the Kullback–Leibler divergence between the true and approximate posteriors. The VBEM (Variational Bayes Expectation Maximization Algorithm) furnishes an optimal (or near optimal) estimate of (v, v ) to make q(θ) ≈ p(θ |g(ω), ˜ m). This approximation which can be obtained iteratively turns out to be often efficient. Therefore, Fm is a lower bound on the log model evidence ln p(g(ω)|m), ˜ which is used to evaluate the relative goodness of models. Besides, the free energy Fm can be interpreted as a combination of a measure of model fit and a measure of model complexity [22]. ˜ l ), l = 1, .., N , the free enGiven ε(ωl , v) = g(ωl , v) − g(ω ergy is [12,18]: Fm = F (v, λ, v ) 1 1 = − Re(ε) (λ)−1 Re(ε)T − Im(ε) (λ)−1 Im(ε)T 2 2



1 4N 1 − ln

(λ) − ln 2π − (v − μ) −1 (v − μ)T 2 2 2 1 1 − ln | | + ln | v | (13) 2 2 where (λ) denotes the (diagonal) covariance matrix built according to (11). The VBEM [18] is carried out to estimate the unknown q(θ):

5

initial v = μθ (populations’ priors are shown in Table 1) until convergence E-step: until convergence ∂g(ω, v) Gv = ∂v ( v )−1 = Re(Gv )T (λ)−1 Re(Gv ) T −1 + Im(G Im(Gv ) + θ  v ) (λ)  T −1 )

(λ) Re(ε) Re(G v −1

v = −( v ) + Im(Gv )T (λ)−1 Im(ε) + θ (v − μθ end M-step: until convergence ∂Fm Fλ = ∂λ ∂ 2 Fm Fλλ = ∂λ2 −1

λ = −Fλλ Fλ end Fm = F (v, λ, v ) (use of (13)) end (14) According to the VBEM algorithm, we can get the estimation of v, λ and v , and using (13), we obtain the free energy values of the different model structures. Then, we use, as an index, the Log Bayes Factor ln Brs :     ln Brs  ln p g(ω)|m ˜ = r − ln p g(ω)|m ˜ =s ≈ F r − Fs

(15)

For two models (model r and model s) to detect the optimized model (the most likely model is the one with largest free energy). This process is called model selection or Bayesian model comparison [22]. For ln Brs > 3, there is strong evidence [23] in favor of model r over model s. 4. Experimental results 4.1. Results To simulate EEG observations, we considered three scenarios: the first one (scenario 1) corresponds to independent populations, the second one (scenario 2) to unidirectional coupled connectivity and the third one (scenario 3) to a bidirectional coupled connectivity scheme. We introduced four models, corresponding to four kinds of connectivity direction (Fig. 2): – Model 1: no connection (Kkj = Kj k = 0) – Model 2: connection from population j to population k (Kkj > 0, and Kj k = 0) – Model 3: connection from population k to population j (Kj k > 0, and Kkj = 0) – Model 4: bidirectional connection (Kkj > 0 and Kj k > 0). Scenarios 1, 2 and 3 correspond to an instance of Models 1, 2 and 4, respectively. For each scenario the four models were tested. Note that when a population is autonomous (i.e. not influenced by the other one), as it is the case for one population

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.6 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

6

in Models 2 and 3 and for both populations in Model 1, the set of θ parameters of the global two-population model includes at most one component θK (at least one coefficient K (Kj k or Kkj ) is known and equal to zero). For the three scenarios, the parameters vectors values were respectively: [Aj Ak Bj Bk Gj Gk Kkj Kj k αi αk ]1 = [5 3.67 3 2.3 20 22.45 0 0 900 900] [Aj Ak Bj Bk Gj Gk Kkj Kj k αi αk ]2 = [5 3.5 3 3.5 20 36 300 0 900 900] [Aj Ak Bj Bk Gj Gk Kkj Kj k αi αk ]3 = [5.3 5.3 2.6 2.6 25 25 300 300 900 900] These parameters values were selected to obtain simulated EEG signals whose spectra were consistent with those observed at the beginning of temporal epilepsy seizures and correspond to the models introduced in [24]. All models were stable, i.e. for each of them, when setting αj and αk to zero, after a transient period, the state trajectory tends towards a fixed point. The case of autonomous nonlinear neural oscillators (i.e. spontaneously oscillating without any external stimulation) is not considered in the models we were interested in. For each scenario, in Fig. 3 are displayed one simulated pair of EEG time signals and the spectra and cross spectra either measured by the SPM Spectral Toolbox or obtained theoretically using (8) from the true parameters values. Fig. 4 illustrates results obtained for one trial with scenario 1. On the top left, we displayed the means (grey boxes)

and ±1.6 std intervals (pink bars) for Model 1. Posterior parameters distributions are depicted as well as the true values of the parameters (blue crosses). The convergence of VBEM algorithm is depicted on the top right of Fig. 4 where are plotted free energies values over iterations. Generally the free energy evolution does not decrease. Now, in the case of Model 3, we can see on the figure an example of an irregular evolution arising when the estimated parameters vector leads to an unstable model which must be stabilized as explained in Section 2.2. On the bottom left of Fig. 4, the spectra and cross spectra estimated by the SPM Spectral Toolbox on the simulated observations are compared with the corresponding theoretical spectra and cross spectra computed by (8) using either the true parameters or the estimated ones. Finally, on the bottom right of Fig. 4, for the same trial, we represented the difference between the four free energy values obtained with the 4 models and the lowest one, corresponding here to Model 4. In this example, Model 1 is the selected model, which is in accordance with the tested scenario. For the 3 scenarios, 11 independent trials were simulated, and the four free energies corresponding respectively to the 4 models were maximized. For a given simulated observation (a pair of simulated EEG signals) the free energies F1 , F2 , F3 , F4 corresponding respectively to the 4 models were computed and ranked. The selected index model was given by the value of m corresponding to the maximal value in the set {Fm , m = 1, 2, 3, 4}. In Table 2, we reported not only the model selection results but also the sum of the free energies over trials (in brackets). The rule given by Fr ≥ Fs + 3 ⇒ (Fr is preferred to Fs ) was applied to obtain the results displayed in this table. Tables 3, 4 and 5 give mean values and standard deviations of

Fig. 3.1. Model 1 (Scenario 1)

Fig. 3. Simulated signals and corresponding spectra. (a) Simulated EEG signals, (b) PSD measured from the simulated signals (Mu is the measured PSD of signal yu , Mj k is the measured cross PSD between signal yj and yk ), (c) theoretical PSD obtained from the true parameters using (8) (Tu is the theoretical PSD of signal yu , Tj k is the theoretical cross PSD between signal yj and yk ).

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.7 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

7

Fig. 3.2. Model 2 (Scenario 2)

Fig. 3.3. Model 4 (Scenario 3)

Fig. 3. (continued)

the parameters obtained with Model 1 for scenario 1, Model 2 for scenario 2 and Model 4 for scenario 3 respectively, each computed on 11 trials. 4.2. Discussion First of all, it is clear that the number of trials (namely 11) to test the VBEM algorithm is very low to assess accurately the

means and dispersions of the estimated parameters and to provide entire credence to the models’ classification performance. Nevertheless, these results give interesting trends which can be discussed compared with results obtained in [15] where the same models were introduced to test the possibility of detecting different types of relations (no relation, uni/bi-directional relation) when using a nonlinear cross-correlation function. This function corresponds to the computation of a nonlinear corre-

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.8 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

8

Fig. 4. An illustration on one trial in the case of independent signals. (a) True and estimated parameters [θAj θAk θBj θBk θGj θGk θαj θαk ] of the physiological model (Model 1). The blue crosses are the true values, the conditional or posterior means are depicted as grey bars, and the pink bars correspond to 90% confidence intervals. (b) VBEM convergence for the four models (scenario 1). (c) (i) Real part of the power spectral densities. (ii) Imaginary part of the power spectral densities. (iii) Modulus of the power spectral densities. Mj : PSD measured from the simulated signal yj . Tj : theoretical PSD obtained from the true parameters. Ej : theoretical PSD obtained from the estimated parameters. Mj k : cross PSD measured from the simulated signals yj and yk . Tj k : theoretical cross PSD obtained from the true parameters. Ej k : theoretical cross PSD obtained from the estimated parameters. (d) Result of Bayesian model comparison. Log Bayes Factors are plotted relative to the worst model (Model 4).

Table 2 Model selection results on 11 trials. Scenario

Accuracy (sum of free energies over trials)

1 2 3

Model 1

Model 2

Model 3

Model 4

7/11 (−3092.5) 1/11 (−4167.2) 1/11 (−4803.6)

1/11 (−4053.6) 7/11 (−4077.5) 2/11 (−5215.8)

1/11 (−3355.6) 1/11 (−4921.9) 0/11 (−5458.6)

2/11 (−3853.7) 2/11 (−4639.7) 8/11 (−3895.5)

Table 3 Estimation of the parameters P for Model 1. Mean values (Mean_P) and standard deviations (Std_P) of the parameters obtained with Model 1 (scenario 1). For comparison, the true values for scenario 1 are also given (True_P). Parameter P

Model 1 Aj

Ak

Bj

Bk

Gj

Gk

αj

αk

True_P Mean_P Std_P

5 4.897 0.1222

3.67 5.05 0.1350

3 3.01 0.0273

2.3 3.02 0.0852

20 19.88 0.1072

22.45 20.156 0.3353

900 1188.79 399.2

900 499.87 95.70

JID:IRBM AID:378 /FLA

[m5+; v1.211; Prn:9/10/2015; 8:57] P.9 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

9

Table 4 Estimation of the parameters P for Model 2. Mean values (Mean_P) and standard deviations (Std_P) of the parameters obtained with Model 2 (scenario 2). For comparison, the true values for scenario 2 are also given (True_P). Parameter P

Model 2 Aj

Ak

Bj

Bk

Gj

Gk

Kkj

αj

αk

True_P Mean_P Std_P

5 4.8780 0.2808

3.5 5.800 1.9185

3 3.017 0.0401

3.5 2.4237 0.5322

20 19.988 0.5262

36 26.598 3.9071

300 897.566 864.904

900 1198.08 470.504

900 593.15 249.760

Table 5 Estimation of the parameters P for Model 4. Mean values (Mean_P) and standard deviations (Std_P) of the parameters obtained with Model 4 (scenario 3). For comparison, the true values for scenario 3 are also given (True_P). Parameter P

Model 4 Aj

Ak

Bj

Bk

Gj

Gk

Kkj

Kj k

αj

αk

True_P Mean_P Std_P

5.3 4.8434 0.7008

5.3 4.6521 0.3330

2.6 2.1905 0.2619

2.6 2.3992 0.3785

25 24.9951 2.3234

25 25.8266 1.7104

300 283.592 85.4786

300 356.775 94.0359

900 1333.65 643.912

900 1442.06 476.725

(t+τ )−f (X(t))) lation index h2XY defined by h2XY (τ ) = 1 − var(Yvar(Y (t+τ )) where τ is the delay variable and f (.) is a piecewise affine function approximating the nonparametric nonlinear regression function from X(t) to Y (t + τ ). The values of this function are in the interval [0; 1]. In [15], when h2XY (τ ) presents a significant maximum positive value for a significantly positive τ value, it was inferred that X drives Y . If the maximal value is too small, the hypothesis of unidirectional influence from the population generating X to the population generating Y is rejected. A conclusion of this simulation experimentation was that the independent case (no link) and the bidirectional case were very difficult to discriminate for these physiological models well suited to the representation of epileptic narrow band fast activities (currently observed in temporal epilepsy). Our preliminary results show that this discrimination is easier using the VBEM approach. Now the estimation of the means of the population’s excitation noises, which has an impact on model stability, was not addressed in this contribution. We tried to set heuristically these means to fixed values, so that the set of qualitative dynamics obtained when tuning the other estimated parameters allows us to cover the case of proximal field real EEG observations.

5. Conclusion In this paper, we first focused on the VBEM algorithm to compute from frequency domain measures the posterior distribution of the parameters of an epileptic physiological model representing a pair of connected (or not connected) neural populations, then we used Log Bayes Factors to detect the different kinds of connectivity in these models. According to our preliminary experiments, independence, unidirectional and bidirectional connectivities were correctly identified and the corresponding estimated parameters are close to their true values. From the point of view of model discrimination efficiency this is an encouraging result in comparison with some approach such as the nonlinear correlation measure. In an applicative perspective, some progress should be done to decrease com-

putational time. This question will be addressed in the future, as we plan to extend this method to multiple populations, test it on real epileptic EEG signals, and compare it with the nonlinear correlation approach and the Granger causality index based on a non-physiological linear modeling. Acknowledgements Research supported by the NSFC under Grant 31400842, the SRFDP under Grants 20130092120035 and 20120092120036, the SRF for ROCS, SEM, the TF for SOCS, MPC, BY2014127 -11 and by Qing Lan Project. References [1] Engel J Jr. Outcome with respect to epileptic seizures. In: Surgical treatment of the epilepsies. 1993. p. 609–21. [2] Yang C, Le Bouquin-Jeannès R, Faucon G, Shu H. Detecting information flow direction in multivariate linear and nonlinear models. Signal Process 2013;93(1):304–12. [3] Vicente R, Wibral M, Lindner M, Pipa G. Transfer entropy—a model-free measure of effective connectivity for the neurosciences. J Comput Neurosci 2011;30:45–67. [4] Friston KJ, Harrison L, Penny W. Dynamic causal modelling. NeuroImage 2003;19:1273–302. [5] Penny WD, Stephan KE, Mechelli A, Friston KJ. Comparing dynamic causal models. NeuroImage 2004;22:1157–72. [6] David O, Guillemain I, Saillet S, Reyt S, Deransart C, Segebarth C, et al. Identifying neural drivers with functional MRI: an electrophysiological validation. PLoS Biol 2008;6:2683–97. [7] Stephan KE, Kasper L, Harrison LM, Daunizeau J, den Ouden HE, Breakspear M, et al. Nonlinear dynamic causal models for fMRI. NeuroImage 2008;42:649–62. [8] Rowe JB, Hughes LE, Barker RA, Owen AM. Dynamic causal modelling of effective connectivity from fMRI: are results reproducible and sensitive to Parkinson’s disease and its treatment?. NeuroImage 2010;52:1015–26. [9] David O, Kiebel SJ, Harrison LM, Mattout J, Kilner JM, Friston KJ. Dynamic causal modeling of evoked responses in EEG and MEG. NeuroImage 2006;30:1255–72. [10] Friston K, Henson R, Phillips C, Mattout J. Bayesian estimation of evoked and induced responses. Hum Brain Mapp 2006;27:722–35.

JID:IRBM AID:378 /FLA

10

[m5+; v1.211; Prn:9/10/2015; 8:57] P.10 (1-10)

W. Xiang et al. / IRBM ••• (••••) •••–•••

[11] Moran RJ, Stephan KE, Seidenbecher T, Pape HC, Dolan RJ, Friston KJ. Dynamic causal models of steady-state responses. NeuroImage 2009;44:796–811. [12] Friston KJ, Bastos A, Litvak V, Stephan KE, Fries P, Moran RJ. DCM for complex-valued data: cross-spectra, coherence and phase-delays. NeuroImage 2012;59:439–55. [13] Bressler SL, Seth AK. Wiener–Granger causality: a well established methodology. NeuroImage 2011;58:323–9. [14] Yang C, Le Bouquin-Jeannès R, Faucon G. Determining the flow direction of causal interdependence in multivariate time series. In: 18th European signal processing conference (EUSIPCO). 2010. p. 636–40. [15] Wendling F, Chauvel P, Biraben A, Bartolomei F. From intracerebral EEG signals to brain connectivity: identification of epileptogenic networks in partial epilepsy. Front Syst Neurosci 2010;4:154. http://dx.doi.org/10.3389/fnsys.2010.00154. [16] Yang C, Le Bouquin-Jeannès R, Bellanger J-J, Shu H. A new strategy for model order identification and its application to transfer entropy for EEG signals analysis. IEEE Trans Biomed Eng 2013;60:1318–27. [17] Wendling F, Hernandez A, Bellanger J-J, Chauvel P, Bartolomei F. Interictal to ictal transition in human temporal lobe epilepsy: insights

[18]

[19] [20]

[21]

[22] [23] [24]

from a computational model of intracerebral EEG. J Clin Neurophysiol 2005;22:343–56. Friston K, Mattout J, Trujillo-Barreto N, Ashburner J, Penny W. Variational free energy and the Laplace approximation. NeuroImage 2007;34:220–34. Yang C. Contribution à l’analyse de connectivité effective en épilepsie. PhD Thesis. Université Rennes 1; 2012. Wendling F, Bellanger J-J, Bartolomei F, Chauvel P. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biol Cybern 2000;83:367–78. Moran RJ, Kiebel SJ, Stephan KE, Reilly RB, Daunizeau J, Friston KJ. A neural mass model of spectral responses in electrophysiology. NeuroImage 2007;37:706–20. Penny WD. Comparing dynamic causal models using AIC, BIC and free energy. NeuroImage 2012;59:319–30. Kass RE, Raftery AE. Bayes factors. J Am Stat Assoc 1995;90:773–95. Wendling F, Chauvel P. Transition to ictal activity in temporal lobe epilepsy: insights from macroscopic models. In: Stolesz I, Staley K, editors. Computational neuroscience in epilepsy. 2008. p. 356–86.