Wentao Xiang

Ahmad Karfoul

Chunfeng Yang

Huazhong Shu

Régine Le Bouquin Jeannès
email: regine.le-bouquin-jeannes@univ-rennes1.fr

Investigation of two neural mass models for DCM-based effective connectivity inference in temporal epilepsy

Keywords: Epilepsy, Neural mass model, Effective connectivity, Power spectral density, Spectral dynamic causal modelling --

Background and objective: Recently, spectral Dynamic Causal Modelling (DCM) has been used increasingly to infer effective connectivity from epileptic intracranial electroencephalographic (iEEG) signals. In this context, the Physiology-Based Model (PBM), a neural mass model, is used as a generative model. However, previous studies have highlighted out the inability of PBM to properly describe iEEG signals with specific power spectral densities (PSDs). More precisely, PSDs that have multiple peaks around β and γ rhythms (i.e. spectral characteristics at seizure onset) are concerned. Methods: To cope with this limitation, an alternative neural mass model, called the complete PBM (cPBM), is investigated. The spectral DCM and two recent variants are used to evaluate the relevance of cPBM over PBM.

Results: The study is conducted on both simulated signals and real epileptic iEEG recordings. Our results confirm that, compared to PBM, cPBM shows (i) more ability to model the desired PSDs and (ii) lower numerical complexity whatever the method.

Conclusions: Thanks to its intrinsic and extrinsic connectivity parameters as well as the input coming into the fast inhibitory subpopulation, the cPBM provides a more expressive model of PSDs, leading to a better understanding of epileptic patterns and DCM-based effective connectivity inference.

Introduction

Epilepsy is a neurological disease that affects approximately 1% of the world population and is characterized by the repetition of seizures whose frequency and duration are variable [START_REF] Scheffer | Ilae classification of the epilepsies: position paper of the ilae commission for classification and terminology[END_REF]. The treatment for epileptic patients generally begins with medication. Recent studies have shown that up to 70% of children and adults with epilepsy can be successfully treated with Anti-Epileptic Drugs (AEDs) [START_REF] Brodie | Staged approach to epilepsy management[END_REF][START_REF] Wu | Automatic lateralization of temporal lobe epilepsy based on MEG network features using support vector machines[END_REF]. However, there are still about 30% of epileptic patients who cannot be adequately controlled with AEDs, which means that (i) seizures are not (completely) suppressed and/or (ii) patients suffer from strong side-effects due to AEDs [START_REF] Van Mierlo | Functional brain connectivity from EEG in epilepsy: Seizure prediction and epileptogenic focus localization[END_REF]. These patients are often diagnosed as drug-resistant or refractory epileptic patients [START_REF] Zigouris | Aed strategy after refractory epilepsy surgery[END_REF][START_REF] Engel | Outcome with respect to epileptic seizures[END_REF]. For these drug-resistant epilepsies, a surgical intervention can be considered as an alternative solution to remove the Epileptogenic Zone (EZ). This zone has been defined as the area of cortex primarily responsible for the generation of clinical seizures [START_REF] Nair | A critical review of the different conceptual hypotheses framing human focal epilepsy[END_REF]. The standard pre-surgery evaluation includes clinical review, brain imaging with Magnetic Resonance Imaging (MRI), ElectroEncephaloGraphy (EEG) or intracranial EEG (iEEG) signals recordings, neuropsychological and psychiatric assessments [START_REF] Duncan | Brain imaging in the assessment for epilepsy surgery[END_REF]. Thanks to its high temporal resolution, the iEEG technique permits to capture very fast neuronal dynamics related to the epileptic seizure. This technique is considered as gold standard for EZ identification and allows a more precise brain exploration in spite of its invasiveness. In fact, seizure activity is not limited to the EZ but also may be modulated -or be mediated by -distributed dynamics in other involved regions.

To better understand the organization of epileptic seizures in terms of onset and propagation, one may resort to effective connectivity [START_REF] Rubinov | Complex network measures of brain connectivity: uses and interpretations[END_REF] (i.e. causal effects of one neural system over another).

A variety of approaches have been proposed to infer effective connectivity. These methods can be divided into two categories, model-free and model-based ones. In [START_REF] Penny | Modelling functional integration: a comparison of structural equation and dynamic causal models[END_REF], authors stressed that effective connectivity analysis tends to be more hypothesis-driven (model-based) than data-driven (model-free). In this paper, we mainly focus on model-based approaches, which rely on two aspects: the generative model (to describe the observed signals) and the model comparison (to seek the best model structure among a set of predefined model structures based on the model evidence). Recently, among model-based techniques, Dynamic Causal Modelling (DCM) [START_REF] Friston | Dynamic causal modelling[END_REF] is increasingly used in quantifying effective connectivity in the scope of epilepsy. Epileptic seizures always show some particular spectral characteristics, i.e. fast activities also referred as rapid discharges corresponding to power spectral densities (PSDs) with several peaks around β and γ frequency bands (15 ∼ 50 Hz) at the seizure onset [START_REF] Wendling | From intracerebral EEG signals to brain connectivity: identification of epileptogenic networks in partial epilepsy[END_REF]. Therefore, researchers have been tempted to use spectral DCM [START_REF] Moran | Dynamic causal models of steady-state responses[END_REF][START_REF] Friston | DCM for complex-valued data: cross-spectra, coherence and phase-delays[END_REF] to infer causality among brain regions in the context of epilepsy [START_REF] Cooray | Characterising seizures in anti-NMDA-receptor encephalitis with dynamic causal modelling[END_REF][START_REF] Papadopoulou | Dynamic causal modelling of seizure activity in a rat model[END_REF][START_REF] Xiang | Inferring effective connectivity in epilepsy using dynamic causal modeling[END_REF][START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF] to cite a few, under the assumption of the generative model (i.e. neural mass model)-Canonical MicroCircuit (CMC) model proposed in [START_REF] Bastos | Canonical microcircuits for predictive coding[END_REF][START_REF] Wei | Bayesian fusion and multimodal dcm for eeg and fmri[END_REF] or Physiology-Based Model (PBM) proposed in [START_REF] Wendling | From intracerebral EEG signals to brain connectivity: identification of epileptogenic networks in partial epilepsy[END_REF][START_REF] Wendling | Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition[END_REF][START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF][START_REF] Shayegh | A model-based method for computation of correlation dimension, lyapunov exponents and synchronization from depth-eeg signals[END_REF].

Each population in the CMC model comprises four subpopulations: (i) superficial pyramidal, (ii) spiny stellate, (iii) inhibitory interneuron and (iv) deep pyramidal. More precisely, in [START_REF] Cooray | Characterising seizures in anti-NMDA-receptor encephalitis with dynamic causal modelling[END_REF][START_REF] Papadopoulou | Tracking slow modulations in synaptic gain using dynamic causal modelling: validation in epilepsy[END_REF] authors used spectral DCM in the CMC model to infer causal coupling between brain areas in order to evaluate the relative contribution of changes in intrinsic (excitatory and inhibitory) connectivity and endogenous input. This study was conducted using EEG recordings acquired from two female patients with recurrent partial seizures.

Authors concluded that changes in excitatory/inhibitory connections were able to explain the observed seizure activity and that these changes were reproduced over seizures. Later, authors in [START_REF] Cooray | Dynamic causal modelling of electrographic seizure activity using Bayesian belief updating[END_REF] proposed a Bayesian belief updating scheme to accelerate the convergence of the spectral DCM. The efficiency of the proposed approach in terms of model inversion compared to the standard DCM was reported based on EEG seizure recordings. In addition, authors in [START_REF] Papadopoulou | Dynamic causal modelling of seizure activity in a rat model[END_REF] used a Bayesian model reduction-based scheme in [START_REF] Friston | Bayesian model reduction and empirical bayes for group DCM studies[END_REF] to identify the key synaptic (e.g. synaptic time constant in deep pyramidal cells) and connectivity parameters (e.g. self-inhibitory connection in inhibitory subpopulation) that underlie seizure onset in a rat model with induced epilepsy. These different studies [START_REF] Cooray | Characterising seizures in anti-NMDA-receptor encephalitis with dynamic causal modelling[END_REF][START_REF] Papadopoulou | Tracking slow modulations in synaptic gain using dynamic causal modelling: validation in epilepsy[END_REF][START_REF] Papadopoulou | Dynamic causal modelling of seizure activity in a rat model[END_REF], stated that the inhibitory subpopulation plays an important role in explaining the seizure activity under CMC model. As emphasized in [START_REF] Wendling | Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition[END_REF], the onset of seizures is often characterized by the appearance of a fast activity in mesial structures, but this activity cannot be modeled in Jansen's model [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns[END_REF]. Authors in [START_REF] Wendling | Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition[END_REF] proposed another extension of Jansen's model, known as PBM. This model takes into account: (i) the respective role of (somatic) fast and (dendritic) slow inhibitory interneurons and (ii) the possible inhibition-related functional reorganizations that take place in the epileptic tissue. The PBM is particularly suitable for studying neuronal dynamics (especially fast activity) during epileptic seizures. With such classification, authors in [START_REF] Wendling | From intracerebral EEG signals to brain connectivity: identification of epileptogenic networks in partial epilepsy[END_REF][START_REF] Shayegh | A model-based method for computation of correlation dimension, lyapunov exponents and synchronization from depth-eeg signals[END_REF] argued that the PBM can be used to simulate iEEG signals in the hippocampus during the transition from interictal to fast ictal activity. As also stated in [START_REF] Grinenko | A fingerprint of the epileptogenic zone in human epilepsies[END_REF], the fast inhibitory interneurons acting on pyramidal cells can be viewed as a prominent mechanism in seizure triggering.

PBM was considered to infer effective connectivity in the context of epilepsy using spectral DCM [START_REF] Xiang | Inferring effective connectivity in epilepsy using dynamic causal modeling[END_REF] and, to deal with the local maximization issue of the free energy, two strategies have been proposed in [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF], the L-DCM where a local adjustment strategy scheme was utilized and the D-DCM where a deterministic annealing scheme was employed. However, PBM is difficult to fit PSDs with multiple peaks around β and γ rhythms and presents a high computational complexity due to the calculation of equilibrium point in Jacobian matrix [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF] related to the predicted PSD.

Later, Zavaglia et al. [START_REF] Zavaglia | A neural mass model to simulate different rhythms in a cortical region[END_REF] as well as Ursino et al. [START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF] proposed a complete Physiology-Based Model, abbreviated hereafter by cPBM, to simulate the multiple/multimodal PSD in cortical regions for simple motor tasks by adding a new negative self-loop as well as a new afferent input to the fast inhibitory interneurons. Also, the cPBM comprises a new extrinsic connectivity parameter modelling direct connections from pyramidal neurons in one population to fast inhibitory interneurons in another population. As indicated in [START_REF] Zavaglia | A neural mass model to simulate different rhythms in a cortical region[END_REF][START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF][START_REF] Cona | A neural mass model of interconnected regions simulates rhythm propagation observed via TMS-EEG[END_REF], the cPBM is much more suitable to produce two distinct resonance peaks and may be more consistent with the physiological reality than the PBM. To the best of our knowledge, cPBM has never been introduced in the context of epileptic seizures and DCM-based framework. In addition, contrary to PBM, cPBM sounds to be more attractive from a computational point of view as the system identification is performed through a linearization of the generative model and output equations around the null equilibrium point [START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF]. Therefore, thanks to the abovementioned properties that cPBM enjoys, this model is considered in the current study.

Following [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF] and thanks to the ability of cPBM to simulate intracranial EEG signals having PSDs with multiple peaks, the motivation of this study is to investigate a potential relationship between the estimation of effective connectivity and the underlying generative model. Therefore, the rest of this paper is organized as follows. A detailed description of the cPBM model is given in section 2. Section 3 presents a full explanation of the cPBM model in the framework of spectral DCM, along with two recently proposed DCM algorithms, L-DCM and D-DCM algorithms [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF]. Section 4 is dedicated to numerical experiments where a comparative study between DCM, L-DCM and D-DCM algorithms in the cPBM framework was firstly conducted using simulated iEEG data. Secondly, the three algorithms were compared in both PBM and cPBM frameworks using real iEEG signals. Finally, section 5 presents discussion.

The neural mass model

Complete Physiology-Based Model

Similarly to the PBM [START_REF] Wendling | Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition[END_REF][START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF], each population in the cPBM consists of four neural subpopulations: pyramidal neurons (P p ), excitatory interneurons (P e ), slow inhibitory interneurons (P s ) and fast inhibitory interneurons (P f ), respectively. Each subpopulation represents a group of neurons of the same type, which approximately share the same membrane potential and can be lumped together. These subpopulations are interconnected with excitatory {C ep ll , C pe ll , C sp ll , C f p ll } and inhibitory {C ps ll , C f s ll , C pf ll , C f f ll } intrinsic connections, as shown in Fig. 1. The differential equations governing cPBM [START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF] describe the dynamic of a set of hidden states (e.g. x l (t) = [x 1,l (t), • • • , x 14,l (t)] T ) for a given population pop l , l = l (l is related to another population, pop l ):

ẋ1,l (t) = x 2,l (t) ẋ2,l (t) = G e,l ω e,l S(x 3,l (t)) -2ω e,l x 2,l (t) -ω 2 e,l x 1,l (t) ẋ3,l (t) = C pe ll x 5,l (t) -C ps ll x 8,l (t) -C pf ll x 11,l (t) ẋ4,l (t) = x 5,l (t) ẋ5,l (t) = G e,l ω e,l (S(x 6,l (t)) + K pp ll S(x 3,l (t -τ ))/C pe ll + u p l (t)/C pe ll ) -2ω e,l x 5,l (t) -ω 2 e,l x 4,l (t) ẋ6,l (t) = C ep ll x 2,l (t) ẋ7,l (t) = x 8,l (t) ẋ8,l (t) = G s,l ω s,l S(x 9,l (t)) -2ω s,l x 8,l (t) -ω 2 s,l x 7,l (t) ẋ9,l (t) = C sp ll x 2,l (t) ẋ10,l (t) = x 11,l (t) ẋ11,l (t) = G f,l ω f,l S(x 14,l (t)) -2ω f,l x 11,l (t) -ω 2 f,l x 10,l (t) ẋ12,l (t) = x 13,l (t) ẋ13,l (t) = G e,l ω e,l (u f l (t) + K f p ll S(x 3,l (t -τ ))) -2ω e,l x 13,l (t) -ω 2 e,l x 12,l (t) ẋ14,l (t) = C f p ll x 2,l (t) -C f s ll x 8,l (t) -C f f ll x 11,l (t) + x 13,l (t) (1)
Comparing PBM with cPBM as shown in Fig. 1, we can note some important differences as summarized in Table 1 [START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF][START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. The values of the parameters specifying Eq. ( 1) are given in Table 2. According to Fig. 1, the two extrinsic random inputs, u p l (t) = σ p l w p l (t) and u f l (t) = σ f l w f l (t) in population pop l , are independent and stand for white noises where w p l (t) and w f l (t) follow Gaussian probability distributions N (0, α p l ) and N (0, α f l ) respectively. To accommodate the amplitude of u p l (t) and u f l (t), we set the diffusion parameters σ p l and σ f l to [START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF] with (a) interactions between four neuronal subpopulations (e.g. pyramidal neurons (P p ), excitatory interneurons (P e ), slow inhibitory interneurons (P s ) and fast inhibitory interneurons (P f ), respectively). (b) block diagram for each hidden state as described in Eq. ( 1), where excitatory and inhibitory intrinsic connections {C ep ll , C pe ll , C sp ll , C f p ll } and {C ps ll , C f s ll , C pf ll , C f f ll }, respectively. S (•) is the sigmoid function, u p l (t) and u f l (t) are the inputs, y l (t) is the output.

20

√ 2 and 10 as shown in Table 2.

The output of cPBM for the l-th population is given by:

y l (t) = x 3,l (t) (2)
Table 1. Comparison of PBM and cPBM. Note that u f l (t) = σ f l w f l (t) and

w f l (t) ∼ N 0, α f l . Case PBM cPBM Sigmoid (S (x (t))) 2e0 1+exp(r0(V0-x(t))) 2e0 1+exp(r0(V0-x(t))) -e 0 Self-loop in P f (C f f ll ) × √ Input in P f (u f l (t)) × √ Extrinsic connection from P p to P f (K f p ll ) × √ Delay (τ ) × √
where y l (t) reflects the mean membrane potential in pyramidal cells which induces the measured iEEG signal, recorded in an invasive way with high temporal resolution.

Table 2 lists the parameters priors used in Eq. ( 1) and Eq. ( 2), such prior distributions are introduced for model inversion (e.g. parameters estimation described in section 3.1.2) based on a Bayesian approach. For simplicity, these parameters can be divided into two categories depending whether the variance σ 2 θi is equal to zero or not: (i) fixed parameters (for which the prior variances are equal to 0), i.e. Φ 1 = {C ep ll , ..., C f f ll , σ p l , σ f l , e 0 , r, v 0 , τ }, and (ii) free parameters i.e. Φ 2 = {ω e,l , ω s,l , ω f,l , G e,l , G s,l , G f,l , α p l , α f l , K pp ll , K f p ll }, which, instead of being estimated directly, are estimated through the vector θ l = θω e,l , θω s,l , θω f,l , θG e,l , θG s,l , θG f,l , θ

α p l , θ α f l , θ K pp ll , θ K f p ll T
. The i-th component, θ i , of vector θ l is defined such that the parameter P i = V i exp(θ i ) with P i ∈ Φ 2 (see Table 2). In fact, the free parameters Φ 2 are expected to vary during a transition from background brain activity to an epileptic seizure activity. Table 2. Model parameters for pop l , l = l , corresponding to cPBM, where the values in red are extracted from [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF], the values in blue from [START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF], and the values in black correspond to our proposal. Prior distributions of P i = V i exp(θ i ), where V i is a fixed default parameter value and Therefore, for a given pair of observed signals, identifying the best model structure or family in the model space is the issue we have to deal with. Such issue can be solved by the DCM-based methods (see section 3).

θ i ∼ p(θ i ) = N µ i , σ 2 θi . If σ 2 θi =
C f p ll Pp → P f 54 N (0, 0) C f s ll Ps → P f 27 N (0, 0) C pf ll P f → Pp 54 N (0, 0) C f f ll P f → P f 27 N (0,

Methods

Spectral DCM

The predicted PSD

The cPBM for two potentially coupled populations represented by Eq. ( 1) and (2) can be described in a set of (2×n s ) first-order delay differential equations (n s = 14 is the number of state variables) and a two-component output vector

y(t) = [y j (t), y k (t)] T , as follows:    ẋ (t) = f (x(t, Ψ), θ) + Du(t) y(t) = Qx(t) (3)
where

x(t) = vec([xj(t), x k (t)] T ) ∈ R 2ns×1 , x l (t) = [x 1,l (t), • • • , x ns,l (t)] T ∈ R ns×1 is the state vector associated to pop l , l ∈ {j, k}. The nonlinear function f (x(t, Ψ), θ) is given by f (x(t, Ψ), θ) = vec([fj (x(t, Ψ), θj) , f k (x(t, Ψ), θ k )] T ) ∈ R 2ns×1 . The
k pop j pop f j u t       p j u t       p k u t       f k u t      
k pop j pop f j u t       p j u t       p k u t       f k u t      
(K jk = K pp jk , K f p jk T , K kj = K pp kj , K f p kj T ).
components of f l (x(t, Ψ), θ l ) ∈ R ns×1 , l ∈ {j, k}, are the restrictions of the right terms of Eq. ( 1) on the restricted domains induced respectively by u p l (t) = 0 and u f l (t) = 0. Note that the first equation in Eq. ( 3) is a vector representation of the following component-wise differential equation:

ẋi (t) = f i (x(t, Ψ(i, :)), θ) + (Du) i (t) (4)
where

∀i ∈ {1, • • • , 2n s }, x i (t), f i (x(t, Ψ(i, :)), θ) and (Du) i (t) denote the i-th component of the vectors x(t), f (x(t, Ψ), θ)
and Du(t), respectively. For the sake of convenience, Ψ(i, :) denotes the i-th row of the delay matrix

Ψ ∈ R 2ns×2ns and x n (t -Ψ i,n
) stands for the n-th component of the 2n sth dimensional vector x(t, Ψ(i, :)). According to Eq. ( 1) and Eq. ( 4), the entries of Ψ are zeros excepted the four entries Ψ 9,6 , Ψ 25,6 , Ψ 10,5 and Ψ 26,5

which are equal to τ (e.g. Ψ 9,6 = τ indicates that x 6 (t -Ψ 9,6 ) is induced by ẋ9 (t) in Eq. ( 4)). As far as the input matrix D ∈ R 2ns×4 is concerned, all its entries are zeros excepted for D 9,1 , D 10,2 , D 25,3 and D 26,4 which are equal to G e,j w e,j /C pe jj , G e,k w e,k /C pe kk , G e,j w e,j and G e,k w e,k , respectively. The input vector u(t) is equal to u p j (t),

u p k (t), u f j (t), u f k (t)
T . The output matrix Q ∈ R 2×2ns is an all-zero matrix excepted for two entries:

Q 1,5 = Q 2,6 = 1.
As we need an analytical expression (at least approximate) of the output PSD, a linearization of Eq. ( 3) around the equilibrium state, following [START_REF] Moran | A neural mass model of spectral responses in electrophysiology[END_REF][START_REF] Moran | Dynamic causal models of steady-state responses[END_REF],

using first-order Taylor expansion is given as follows:

   ẋ(t) = ˜ (θ) x(t) + D 1 w(t) y(t) = Qx (t) (5)
where 3) is equal to the null vector since the sigmoid function and w(t) are centered and the mean of is equal to zero as indicated in [START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF]). D 1 = (I + Ψ • (x 0 , θ)) -1 Ddiag(σ) with diag(σ) denoting a square diagonal matrix whose diagonal corresponds to the elements of the vector σ

˜ (θ) (I + Ψ • (θ)) -1 (θ) [33] (• is Hadamard product) and w(t) = w p j (t), w p k (t), w f j (t), w f k (t) T . Note that, (θ) = ∂f (x(t),θ) ∂x(t) | x(t)=0 ∈ R 2ns×2ns is the Jacobian matrix of f (•, θ), computed at x(t) = 0 (i.e. the equilibrium point of Eq. (
= σ p j , σ p k , σ f j , σ f k T .
Then, the Fourier transform to both sides of Eq. ( 5) is applied:

   i2πν x (ν) = ˜ (θ) x (ν) + D 1 w (ν) y (ν) = Q x (ν) = H(ν, θ) w (ν) (6)
where

H(ν, θ) = Q i2πνI -˜ (θ) -1
D 1 , ν denotes the frequency variable and

H(ν, θ) ∈ R 2×4
is the transfer function (I is the identity matrix). The analytical expression of the PSD obtained from the parameters (named predicted PSD)

is denoted by G (ν, θ) =   G jj (ν, θ) G jk (ν, θ) G kj (ν, θ) G kk (ν, θ)   ∈ R 2×2 and computed from
Eq. ( 6), as follows [START_REF] Moran | Dynamic causal models of steady-state responses[END_REF]:

G (ν, θ) = H (ν, θ) G w (ν, θ) H H (ν, θ) (7)
where

Gw (ν, θ) =         G pp w,jj (ν, θ) G pp w,jk (ν, θ) G pf w,jj (ν, θ) G pf w,jk (ν, θ) G pp w,kj (ν, θ) G pp w,kk (ν, θ) G pf w,kj (ν, θ) G pf w,kk (ν, θ) G f p w,jj (ν, θ) G f p w,jk (ν, θ) G f f w,jj (ν, θ) G f f w,jk (ν, θ) G f p w,kj (ν, θ) G f p w,kk (ν, θ) G f f w,kj (ν, θ) G f f w,kk (ν, θ)         =         α p j 0 0 0 0 α p k 0 0 0 0 α f j 0 0 0 0 α f k         (8)
denotes the PSD matrix associated with the white Gaussian noise vector w(t).

Since the inputs w p j (t), w p k (t), w f j (t) as well as w f k (t) in pop j and pop k are assumed to be independent in/between populations, G w (ν, θ) is diagonal with diagonal entries, G ξξ w,ll (ν, θ) = α ξ l , ξ ∈ {p, f } and l ∈ {j, k}, denoting the auto PSD of the input w ξ l (t).

Model inversion

Concerning the sample PSD matrix, obtained from the observed iEEG signals, it is computed using the SPM spectral toolbox with a 12-pole AR process [START_REF] Moran | Dynamic causal models of steady-state responses[END_REF].

The sample PSD, denoted by

G (ν) =   Gjj (ν) Gjk (ν) Gkj (ν) Gkk (ν)   ∈ R 2×2 at each
frequency ν, can also be considered as a noisy version of G (ν, θ):

G (ν) = G (ν, θ) + E (ν) (9)
where

E (ν) =   E jj (ν) E jk (ν) E kj (ν) E kk (ν)
  is the noise matrix. Let us consider N frequency bins (ν i , i = 1, 2, ..., N ) of G (ν) stacked in a column vector g, then following Eq. ( 9), g can be written:

g = g (θ) + ε ( 10
) where g = [g jj , gjk , gkj , gkk ] T ∈ R 4N ×1 gll = [ Gll (ν 1 ), ...., Gll (ν N )] T ∈ R N ×1 , {l, l } ∈ {j, k} g(θ) = [g jj (θ), g jk (θ), g kj (θ), g kk (θ)] T ∈ R 4N ×1 g ll (θ) = [G ll (ν 1 , θ), ...., G ll (ν N , θ)] T ∈ R N ×1 ε = [ε jj , ε jk , ε kj , ε kk ] T ∈ R 4N ×1 ε ll = [E ll (ν 1 ), ...., E ll (ν N )] T ∈ R N ×1
ε is the vector of prediction error with ε ∼ N (0, Σ λ ), where the posterior covariance Σ λ is related to a hyper-parameter vector λ under Gaussian assumptions [START_REF] Friston | Variational free energy and the Laplace approximation[END_REF]. According to Eq. ( 10), the vector of parameters θ for a given model structure Mm can be obtained by a Bayesian model inversion step using the model evidence p (g|Mm) [START_REF] Friston | Variational free energy and the Laplace approximation[END_REF][START_REF] Daunizeau | Variational Bayesian identification and prediction of stochastic nonlinear dynamic causal models[END_REF]. The logarithm of model evidence, ln p (g|Mm), can be expressed as:

ln p (g|Mm) =F m + KL (q (θ) , p (θ|g, Mm)) ( 11
)
where KL is the Kullback-Leibler divergence, p (θ|g, Mm) is the posterior distribution of the vector θ, q (θ) ∼ N θ, Σ θ and F m is the well-called free energy (of the model structure Mm). Since KL is non negative, then ∀q (θ) :

ln p (g |Mm ) ≥ F m , the equality being verified if KL (q (θ) , p (θ |g, Mm )) = 0,
i.e. if q (θ) = p (θ |g, Mm ) (almost everywhere), a maximum value of F m can be always seen as a lower bound of ln p (g|Mm). The definition of free energy is given as follows [START_REF] Friston | Variational free energy and the Laplace approximation[END_REF]:

F m = q(θ)lnp(g, θ|Mm) dθ -q(θ)lnq(θ) dθ (12)
Maximizing the free energy, F m , is usually done in an iterative way using a variational scheme [START_REF] Friston | Variational free energy and the Laplace approximation[END_REF], named as the variational EM algorithm 1 .

Model and family selection

Once the above model inversion step is achieved, a model selection step is performed. This step consists in selecting, in the predefined model space, the best model structure underlying the observed data. Recall that an optimal model structure is the one defined in terms of the maximum free energy compared to the free energies of the other considered model structures. Note that causality identification (i.e. effective connectivity identification) in cPBM can also be based on the family identification. Now, this causality identification depends on the vectors of coupling parameters K kj and K jk . As an example, the family is computed. This posterior probability is equal to the sum of model evidence p (g|Mm) times the prior p (Mm|Ff ) [START_REF] Stephan | Bayesian model selection for group studies[END_REF][START_REF] Penny | Comparing families of dynamic causal models[END_REF]:

F2 (i.e. K jk = K pp jk , K f p jk T = 0
p (g|Ff ) = Mm∈Ff p (g|Mm) p (Mm|Ff ) ( 13
)
where the prior p (Mm|Ff ) is given as follows (see Fig. 2 (c)):

p (M1|F1) = 1, M1 ∈ F1 p (Mm|F2) = 1 3 , ∀Mm ∈ F2 = {M2, M5, M6} p (Mm|F3) = 1 3 , ∀Mm ∈ F3 = {M3, M9, M11} p (Mm|F4) = 1 9 , ∀Mm ∈ F4 = {M4, M7, M8, M10, M12, M13, M14, M15, M16} (14)
1 Although we have referred to the variational scheme as a variational EM algorithm, technically it is known as Variational Laplace. Variational Laplace is a generic form of variational Bayes under the Laplace assumption (i.e. under the assumption that the prior and posterior are Gaussian). This enables one to estimate the posterior expectation and covariance of both the parameters θ and hyper-parameter vector λ. Please see [START_REF] Friston | Variational free energy and the Laplace approximation[END_REF] for details.

Similar to the logarithm of model evidence, the family free energy, F F f index which can be defined as lnp (g|Ff ), is given by:

F F f = lnp (g|Ff ) = ln Mm∈Ff p (g|Mm) p (Mm|Ff ) ln Mm∈Ff exp(F m ) + lnp (Mm|Ff ) (15)
where p (g|Mm) exp(F m ) (see Eq. ( 11)). Likewise, the best family in cPBM is identified as the one with the maximum family free energy. The pseudo-code summarizing spectral DCM algorithm is given in Algorithm 1. Note that, the variational EM algorithm stops when the convergence condition is met for eight successive iterations. The convergence condition is fulfilled when the difference in the estimation of F m exhibits a value smaller than or equal to a predefined threshold [START_REF] Brodie | Staged approach to epilepsy management[END_REF] or a maximum number of iterations (60) is reached.

Sub-estimation problem

As shown in Fig. 3 (a), the variational EM algorithm in DCM is highly sensitive to the initialization, so it may lead to a sub-estimation of the free energy and consequently may lead to a mis-estimation of the optimal model structure/family. Repeat until convergence or a maximum number of iterations E-step:

( θ, Σ θ ) ← argmax f ixed λ, Σ λ Fm( θ, λ, Σ θ , Σ λ ) M-step: (λ, Σ λ ) ← argmax f ixed θ, Σ θ Fm( θ, λ, Σ θ , Σ λ ) End Fm = Fm( θ, λ, Σ θ , Σ λ ) (d)

L-DCM

To deal with the sub-estimation problem in DCM, the L-DCM algorithm employing a local adjustment strategy was proposed in [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF], and the scheme of L-DCM to maximize the free energy is shown in Fig. 3 (b). As indicated in [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF],

the sub-estimation, when two populations are considered, relies on a threshold,

T h ll = ν b l ν=ν a l γ 2 gll (ν)
2 and a prediction error,

Er ll = ν b l ν=ν a l g ll ν, θ -gll (ν) 2 ,
l ∈ {j, k}, for each population, where gll (ν) and g ll ν, θ stand for the sample PSD computed using the SPM spectral toolbox with a 12-pole AR process [START_REF] Moran | Dynamic causal models of steady-state responses[END_REF] and the estimation of the predicted PSD obtained by substituting the model parameters estimate θ using the variational EM algorithm, respectively.

According to [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF], we set γ = 0.1, ν a l = ν m l -10 and ν b l = ν m l + 10 with ν m l standing for the frequency bin corresponding to the maximum value of gll (ν). A good estimation of the PSD g ll ν, θ is obtained when Er ll ≤ T h ll , l ∈ {j, k}. However, in case this inequality is not fulfilled, three possible situations can be figured out: situation 1 ○: Er jj ≤ T h jj , Er kk > T h kk ; situation 2 ○:

Er jj > T h jj , Er kk ≤ T h kk ; situation 3 ○: Er jj > T h jj , Er kk > T h kk . The L-DCM algorithm proposed in [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF] refined the variational EM estimation results by performing few supplementary well-conditioned variational EM sweeps with a local adjustment strategy for the initialization. The initial guess for the (r + 1)-th EM sweep θ the model structure (i.e. Mm), ∆θ is given as follows:

(r+1) 0 = vec θ (r+1) 0,j , θ (r+1) 0,k T is then locally adjusted, θ (r+1) 0 = θ(r) + ∆θ, where θ(r) = vec θ(r) j , θ (r)
∆θ =            vec [0, ∆θ k ] T , ∀Mm ∈ {F1, F2} under 1 ○ vec [∆θ j , 0] T , ∀Mm ∈ {F1, F3} under 2 ○ vec [∆θ j , ∆θ k ] T , otherwise (16)
where ∆θ j ∼ N 0, σ 2 I and ∆θ k ∼ N 0, σ 2 I , I is the identity matrix and σ 2 is a parameter that was empirically set to 0.02. The pseudo-code of L-DCM algorithm can be seen in Algorithm 2.

Algorithm 2 L-DCM algorithm

= vec [θ0,j, θ 0,k ] T , λ0, Σ θ 0 , Σ λ 0 (c) Perform the following EM steps:

Repeat until convergence or a maximum number of iterations E-step: Another way to deal with the sub-estimation issue in spectral DCM relies essentially on a deterministic annealing scheme [START_REF] Ueda | Deterministic annealing EM algorithm[END_REF][START_REF] Katahira | Deterministic annealing variant of variational Bayes method[END_REF], leading to a new algorithm called D-DCM hereafter. More precisely, the deterministic annealing scheme, inspired from the thermodynamic theory, is based on tracking the optimum of the objective function under consideration from a high temperature wherein the objective function is smoothed (i.e. it may have one global optimum, see Fig. 3 (c) β = 0.4) to low temperature wherein the objective function returns to its initial form (see Fig. 3 (c) β = 1). By introducing a temperature related parameter, denoted by β (temperature T = 1 β ), the objective function (e.g. the free energy) F m in Eq. ( 12) can be rewritten as:

( θ, Σ θ ) ← argmax f ixed λ, Σ λ Fm( θ, λ, Σ θ , Σ λ ) M-step: (λ, Σ λ ) ← argmax f ixed θ, Σ θ Fm( θ, λ, Σ θ , Σ λ ) End F (r) m = Fm( θ, λ, Σ θ , Σ λ ) (d) If
F m = q(θ)lnp(g, θ|Mm) dθ - 1 β q(θ)lnq(θ) dθ (17)
The variational EM algorithm [START_REF] Friston | Variational free energy and the Laplace approximation[END_REF] is used to perform such maximization in spectral DCM framework. In this way, one guarantees a good initialization

( θ, Σ θ ) ← argmax f ixed λ, Σ λ Fm( θ, λ, Σ θ , Σ λ , β) M-step: (λ, Σ λ ) ← argmax f ixed θ, Σ θ Fm( θ, λ, Σ θ , Σ λ , β) End Fm = Fm( θ, λ, Σ θ , Σ λ , β) (e) βnew = β + const (f) If βnew <= 1
Repeat from step (d) with θ0 = θ from βnew as a new initial value; Otherwise stop (g) Turn to next model structure Mm, repeat from step (a) (h) Model/Family selection to identify the optimal model structure/family 4. Results

Simulations

In a previous study, we have performed a comparative evaluation between the L-DCM, D-DCM and DCM algorithms in inferring effective connectivity among epileptic neural populations at the seizure onset using PBM [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF]. Hereafter, the simulated iEEG data are generated using the cPBM. Two connectivity configurations, the unidirectional and the bidirectional influence, between the two neural populations are considered. Several model structures can fulfill unidirectional and bidirectional coupling configurations between two populations (see Fig. 2 (b)). Therefore, model structures sharing the same coupling pattern (unidirectional/bidirectional) are considered to belong to the same family as shown in Fig. 2 (c). Furthermore, two families are also considered: (i) the family F2 composed of model structures M2, M5 and M6 reflecting unidirectional coupling from pop j to pop k (see section 4.1.1) and (ii) the family F4 composed of model structures M4, M7, M8, M10, M12, M13, M14, M15 and M16 reflecting bidirectional coupling between the two populations. For simplicity, only M8 and M16 were used to generate data. These models were taken to be representative of F4 (see section 4.1.2).

Unidirectional influence

As mentioned previously, the family F2 = {M2, M5, M6} is considered here. Figs. 4 (a-i), (b-i) and (c-i) show the three different model structures (corresponding to M2, M5 and M6 respectively) and their related vector θ. Note that, for each of model structure in F2, we simulated 100 trials of 10-second length signals by integrating delay differential equations [START_REF] David | Dynamic causal modeling of evoked responses in EEG and MEG[END_REF] and using Runge-Kutta 4-th order method to solve Eq. ( 3) at a sampling frequency of 256 Hz and the values of the parameters in θ are chosen in such a way that the model output, y(t), is equivalent to an epileptic iEEG signal with a peak around 25 Hz in its PSD observed at the seizure onset [START_REF] Wendling | From intracerebral EEG signals to brain connectivity: identification of epileptogenic networks in partial epilepsy[END_REF][START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF]. Fig. 4 (a-c) depicts the temporal signature of the simulated iEEG and the sample auto/cross PSDs for M2, M5 and M6, respectively. 0,0,0,0,0,0,0,0 0,0,0.3, 0.3,0, 0.2,0,0,ln 0.7 0.7 ,ln

j k = = - -   j pop k pop M2    
0,0,0,0,0,0,0,0 0,0,0.3, 0.3,0,0.4,0,0, 0.25,

j k = = - -    (a-i) ground truth j pop k pop M2    
0,0,0,0,0,0,0,0 0,0,0.3, 0.3,0,0.4,0,0, 0.25,

j k = = - -    (a-i) ground truth    
0,0,0,0,0,0,0,0 0,0,0.2, 0.3,0,0.1,0,0, , ln 0.7

j k = = -    j pop k pop M5 (b-i) ground truth    
0,0,0,0,0,0,0,0 0,0,0.2, 0.3,0,0.1,0,0, , ln 0.7 influence, the black arrows correspond to excitatory connections, whereas the red arrows correspond to extrinsic connections that target inhibitory populations, the corresponding vector θ = vec [θ j , θ k ] T , where θ j = θω e,j , θω s,j , θω f,j , θ G e,j , θ G s,j , θ G f,j , θ α p j , θ α f j T and

θ k = θω e,k , θω s,k , θω f,k , θ G e,k , θ G s,k , θ G f,k , θ α p k , θ α f k , θ K pp kj , θ K f p kj T .

Model selection

Table 3 shows the identification rates2 computed over 100 trials for the three methods in identifying model structures M2, M5 and M6 when these models were successively considered as a ground truth. For each algorithm and each simulated model structure, the identification rates of only the four most plausible model structures among the sixteen structures depicted in Fig. 2 (b) are reported in this table. Regarding the ground truth M2, both D-DCM and L-DCM algorithms are able to identify the true model structure (e.g. M2) with identification rates of 48% and 38% respectively. This fact does not hold true for the DCM algorithm where the M16 model structure with an identification rate of 59% is rather retained than the M2 model structure. As far as ground truth M5 is concerned, all algorithms are able to identify the true model structure (e.g. M5). However, the DCM algorithm shows a slight superiority (76%) especially compared to the D-DCM algorithm (71%). When the true model structure is M6, D-DCM, DCM and L-DCM algorithms succeed in identifying M6 with identification rates of 99%, 97% and 85%, respectively.

To sum up, false detection in the context of model identification is generally highly related to a spurious estimation of one or more coupling links between neural populations under study (more details are presented in section Support materials S1). For instance, when M2 is the ground truth, estimating M6 instead of M2 for DCM, L-DCM and D-DCM approaches is mainly due to the fact that estimation errors during the optimization process could lead to some small non-null value of the coupling parameter K f p jk (the red line in Fig. 2 (a)). To deal with such problem, spurious connections can be removed through a thresholding step (more details are presented in Appendix: B). On the other hand, in the case of false detection, one can interestingly note that, even if a wrong model structure is identified (e.g. M6 instead of M2 in our example), this structure generally belongs to the same family as the right model structure.

This shows that the algorithms are still able to quantify the direction of coupling between populations under study (e.g. M6 and M2 belong to the same family F2, see Fig. 2 (c)) which is the main goal of effective connectivity. Therefore, resorting to family identification analysis would provide a new insight regarding the estimation of effective connectivity between the neural populations.

The performance of the three algorithms is also evaluated using the coefficients of determination R 2 L-DCM , R 2 D-DCM and R 2 DCM , which stand for the measures of the goodness of fit, computed between the sample and estimated PSDs, as depicted in Figure 4 (a-c). This figure shows comparable behaviors of the three algorithms with slight superiority of the L-DCM one, since the L-DCM algorithm applies a local adjustment strategy based on the estimation errors [START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF].

Family selection

As just mentioned, an analysis of family identification could be preferable to circumvent the false detection encountered in the model identification due to spurious estimation of coupling parameters. This idea comes from the observation that some wrong identified model structures belong to the same family as the target model structure. Since the main goal of effective connectivity is to point out the influence of one neural population on another one, it could be interesting to quantify, in addition to model identification, the ability of the algorithms to identify the correct family. Table 4 shows both the family identification rate and the execution time (in minutes) over the 100 trials for L-DCM, D-DCM and DCM methods for the previous ground truths M2, M5 and M6, respectively. According to this table, the D-DCM algorithm shows higher performance in identifying the true family (e.g. F2) than the DCM and L-DCM algorithms for M2 and M6. Indeed, for M2, the D-DCM algorithm provides an identification rate of 83% compared to 73% and 18% for the L-DCM and DCM approaches respectively.

Furthermore, the D-DCM algorithm provides, for M6, an identification rate of 99% compared to 97% and 87% for the DCM and L-DCM approaches respectively.

Regarding ground truth M5, all algorithms succeed in identifying the true family (e.g. F2 family) with a slight preference for DCM. As far as the computation time is concerned, both L-DCM and D-DCM algorithms show, as expected, higher computation time compared to DCM in Table 4. In fact, this higher computation time is essentially due to the additional EM sweeps in the L-DCM (see Fig. 3 (b)) and to the adopted deterministic annealing scheme in D-DCM (see Fig. 3 (c)).

This result is assessed also in terms of the family free energy as depicted in Fig. 5. According to Fig. 5 (a) for ground truth M2, both L-DCM and D-DCM approaches show higher family free energies compared to the DCM algorithm. More particularly, for the true family (e.g. F2), the D-DCM and L-DCM algorithms provide a higher family evidence than the DCM one. This figure together with Table 4 highlight the fact that the DCM algorithm suffers from higher false family identification rate. Regarding the correct model M5 and M6 in Figs. (b-c), all methods are able to identify the true family F2 with comparable family identification rates.

Bidirectional influence

Family F4 comprises all model structures reflecting bidirectional coupling between pop j and pop k , as shown in Fig. 2 (c). The performance of L-DCM, D-DCM and DCM in identifying F4 is evaluated in three different ground truths presented in Fig. 6 (a-i), (b-i) and (c-i), where the model structure M8 and two variants of the model structure M16 were used. More precisely, symmetric extrinsic couplings (e.g. K pp jk = K pp kj , K f p jk = K f p kj , see Fig. 6 (b-i) for the values of the parameters) between the two populations were used giving rise to M16 s and asymmetric extrinsic couplings (e.g.K pp jk = K pp kj , K f p jk = K f p kj , see Fig. 6 (c-i) for the values of the parameters) were used giving rise to M16 as . For each ground truth, 100 trials of 10-second length signals at a sampling frequency of 256 Hz were simulated. The vector of parameters θ in model structure M8 and M16 s is tuned such that the model output, y(t), is a fast sinusoidal activity with only one peak around 25 Hz (Fig. 6 (a-b)) similar to the epileptic fast activity [START_REF] Wendling | From intracerebral EEG signals to brain connectivity: identification of epileptogenic networks in partial epilepsy[END_REF][START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF]. Regarding the symmetric nature of M16 s , the vectors of parameters θ j and θ k associated to pop j and pop k , respectively, are effectively equal. However, these two vectors are not identical in the case of ground truth M16 as . In Fig. 6 (c) we represent, under M16 as configuration, the epileptic-like activity y(t), the auto (cross) PSDs of the two populations displaying two peaks (around 20 Hz and 28 Hz), similar to possible activity at the seizure onset.

Model selection

Table 5 gives the four most plausible model structures and the model identification rate over 100 Monte Carlo trials for the bidirectional influence under the three ground truths M8, M16 s and M16 as , respectively. According to this table, all techniques globally succeed in identifying the correct model structure. Furthermore, L-DCM outperforms D-DCM and DCM in the ground truth M8 (85% identification rate vs 80% and 79% respectively) whereas the D-DCM approach outperforms DCM and L-DCM in the ground truth M16 s (96% identification rate vs 94% and 80% respectively). Finally, all methods can identify M16 with a 100% identification rate in the ground truth M16 as . More

j pop k pop M16s    
0,0,0,-0.4,0,0,0,0,ln0.5, 0,0,0.1, 0.4,0.1,0.2, 0 0,0,0, influence, the black arrows correspond to excitatory connections, whereas the red arrows correspond to extrinsic connections that target inhibitory populations, the corresponding vector θ = vec [θ j , θ k ] T , where θ j = θω e,j , θω s,j , θω f,j , θ G e,j , θ 6 (a-c).

G s,j , θ G f,j , θ α p j , θ α f j , θ K pp jk , θ K f p jk T and θ k = θω e,k , θω s,k , θω f,k , θ G e,k , θ G s,k , θ G f,k , θ α p k , θ α f k , θ K pp kj , θ K f p kj T .

Family selection

Regarding the family selection and the averaged execution time for each family over 100 trials, they are displayed in Table 6. All techniques perfectly identify the correct family. Furthermore, Fig. 7 depicts the boxplots of family free energy over the 100 trials using the three methods. Clearly, the maximal median of the family free energy is globally obtained for F4 whatever the method and the ground truth. Furthermore, D-DCM show relatively higher boxplots of family free energy than L-DCM and DCM for all ground truths and families. Let

Regarding the coefficients of determination between the sample and estimated

PSDs and considering the estimated PSD of pop j as given in Fig. 8 (c-i) (d-i),

under the cPBM framework, a higher estimation quality is obtained for the ). This is due to the fact that, under PBM framework, no method is able to reproduce a PSD with two peaks (around 22 Hz and 44 Hz) corresponding to the sample PSD computed from the observed iEEG signals.

On the contrary, the cPBM allows to recover two peaks in the PSD whatever the method, this PSD being close to the sample one. We can also observe that, in the case of cPBM, the D-DCM method provides the best performance in the population pop j , as shown in Fig. 8 (d-i). Regarding the estimated PSD of pop k , we also note from Fig. 8 (c 7, all methods show, as expected, an execution time that is considerably smaller under the cPBM compared to the PBM. Around 65%, 61% and 41% of averaged execution time is saved for L-DCM, D-DCM and DCM, under the assumption of cPBM compared to the case of PBM. This is due to the fact that, contrary to PBM, the equilibrium point for (θ) in Eq. ( 5) in cPBM to derive the estimated PSD during the optimization process does not need to be computed. This is since both the sigmoid function and the population inputs are centered in cPBM. Consequently, cPBM provides a faster model inversion than PBM even if the number of model parameters and model structures are higher in cPBM compared to PBM.

Discussion

Conclusions

We attempted to contribute to the analysis of brain effective connectivity in the context of drug-resistant epilepsy in the framework of neural mass models.

This work focused on the model structures used to simulate intracranial epileptic activities similar to the ones observed at the seizure onset, and more particularly on the fast activity around 25Hz (sinusoidal dynamics) [START_REF] Wendling | From intracerebral EEG signals to brain connectivity: identification of epileptogenic networks in partial epilepsy[END_REF].

Regarding the simulated iEEG signals using cPBM, the D-DCM algorithm showed some superiority in terms of identification when an extrinsic coupling from the pyramidal cells to the slow inhibitory interneurons exists (case of unidirectional influence) between two populations. In the case of bidirectional coupling, the three techniques we presented, performed well. The D-DCM and L-DCM algorithms outperformed DCM in terms of maximizing the free energy in cPBM, with a slight superiority in favour of D-DCM over L-DCM.

Regarding the real epileptic signals, the analysis was conducted on a couple of iEEG channels. Results showed that the potential causal influence was properly identified using either PBM or cPBM. Compared to PBM, cPBM offered a lower computational complexity since computing the equilibrium point was not any more required. This leaded to a faster model inversion even if the number of model parameters and model structures was relatively high. Moreover, the cPBM provided a more expressive model of PSDs and a better estimation accuracy compared to PBM as already stated in our obtained results on real epileptic iEEG signals.

Limitations

Clearly, there are still some issues that we have not addressed in this paper. Firstly, the D-DCM approach can be improved by incorporating a local adjustment initialization strategy as the one in L-DCM. Secondly, our study only concerned the causal coupling between two neural populations and it would be interesting to extend the L-DCM and D-DCM methods to the case of several coupled neural populations or larger network implementations since potentially more than two distant populations can be involved in the seizure initiation. Let us note that the scalability of the cPBM model for larger network implementations may be probably solved by the Baysian model reduction [START_REF] Friston | Bayesian model reduction[END_REF]. In addition to the direct coupling between several coupled neural populations, identifying indirect causal coupling could also improve effective connectivity inference. Finally, since epileptic seizures are patient-dependent, the performance of the L-DCM and D-DCM algorithms should be assessed by an extensive study using a larger dataset. A comparison to other methods, such as Granger Causality [START_REF] Seth | A MATLAB toolbox for Granger causal connectivity analysis[END_REF], transfer entropy [START_REF] Yang | A new strategy for model order identification and its application to transfer entropy for EEG signals analysis[END_REF] and DCM approaches could also be investigated. However, these different types of approaches do not address the same causality concept and are still disputed. expected. This is mainly due to the fact that the difference between the median of the maximized free energy of M6, M10 or M16 models and the one of the M2 one is significantly small (see Fig. 11). Now, the DCM shows higher false model identification rate compared to the L-DCM and D-DCM methods. This can be interpreted by the highest median of the maximized free energy obtained for the M16 model. Regarding the correct model M5, as the medians of the maximized free energy of models M5, M6, M7 and M8 are to some extent comparable (see Fig. 12), with slight superiority for M5, some false model identification can occur for the three algorithms. As far as the ground truth M6 is concerned, Table 3 clearly shows very high performance in terms of detection rate, especially for D-DCM (99%) and DCM (97%). These results can also be interpreted in terms of the median of the maximized free energy depicted in Fig. 13. According to this figure, the median of the maximized free energy obtained for identifying the correct model, M6, for D-DCM and DCM algorithms, is significantly higher than the ones related to the other model structures. Consequently, a very low false detection rate for these algorithms is expected.

As relatively high false detection rate is to some extent problematic, investigating the model evidence, p(g|Mm), would help to complete the whole scene regarding the ability of each algorithm in identifying the model structure.

This model evidence is nothing else than the model posterior probability under flat prior and Mm [START_REF] Friston | Bayesian model reduction[END_REF]. This quantity is computed via the softmax function of the median of the maximized free energy computed over the 100 conducted trials, that means: p(g|Mm) = exp( Fm)

Figure 1 .

 1 Figure1. The population pop l in cPBM[START_REF] Ursino | The generation of rhythms within a cortical region: analysis of a neural mass model[END_REF] with (a) interactions between four neuronal subpopulations (e.g. pyramidal neurons (P p ), excitatory interneurons (P e ), slow inhibitory interneurons (P s ) and fast inhibitory interneurons (P f ), respectively). (b) block diagram for each hidden state as described in Eq. (1), where excitatory and inhibitory intrinsic connections {C ep ll , C pe ll , C sp ll , C f p ll } and {C ps ll , C f s ll , C pf ll , C f f ll }, respectively. S (•) is the sigmoid function, u p l (t) and u f l (t) are the inputs, y l (t) is the output.

Figure 2 .

 2 Figure 2. The effective connectivity between popj and pop k in cPBM. (a) two potentially coupled populations; (b) the model space including the 16 model structures derived from possible extrinsic coupling(s) between the two populations; (c) the family space including four families and corresponding model structures (K jk = K pp jk , K f p

  ) includes three model structures, i.e. M2, M5, and M6 (see Fig.2 (c)), in cPBM with one or two element(s) of the vectorK kj = K pp kj , K f p kj Tthat are non-zero and positive. To decide the most plausible family in cPBM, the family posterior probability of each family Ff , p (g|Ff )

Algorithm 1 (

 1 Spectral) DCM algorithm (a) Considering the current model structure Mm based on θ (b) Initialization θ0 = vec [θ0,j, θ 0,k ] T , λ0, Σ θ 0 , Σ λ 0 (c) Perform the following EM steps:

3 Fm 3 FmFigure 3 .

 333 Figure 3. Scheme of three algorithms to maximize the free energy F m , under considering Mm with (a) DCM, (b) L-DCM and (c) D-DCM (β = 0.4, β = 0.6, β = 0.8 and β = 1) respectively.

  a random vector dropping from zero-mean normal distribution with a predefined variance. By taking into account (i) the encountered situation (i.e. 1 ○ or 2 ○ or 3 ○) and (ii)

( a )

 a Considering the current model structure Mm based on θ (b) Initialization r = 1, θ

=

  θ(r) + ∆θ as a new initial value and r = r + 1; Otherwise stop (e) Fm = max{F (r) m }, r = 1, 2, ..., 6 (f) Turn to next model structure Mm, repeat from step (a) (g) Model/Family selection to identify the optimal model structure/family 3.4. D-DCM

  through decreasing temperature values (e.g. β varying from very small values up to 1) as indicated in Fig.3(c). The pseudo-code of the D-DCM algorithm is given in Algorithm 3. Algorithm 3 D-DCM algorithm (a) Considering the current model structure Mm based on θ (b) Set β ← βint (0 < βint < 1) (c) Initialization θ0 = vec [θ0,j, θ 0,k ] T , λ0, Σ θ 0 , Σ λ 0 (d) Perform the following EM steps:Repeat until convergence or a maximum number of iterations E-step:

Figure 4 .

 4 Figure 4. (a-i), (b-i), (c-i): Three ground truths M2, M5 and M6 for unidirectional

  Simulated time signals displaying epileptic activities (i.e. (a-ii), (b-ii), (c-ii)), corresponding auto PSD amplitudes (i.e. (a-iii, a-iv),(b-iii, b-iv),(c-iii, c-iv)) and cross PSD amplitudes (i.e. (a-v),(b-v),(cv)) for the three ground truths M2, M5 and M6, respectively. The sample PSD, denoted by gn (ν), n ∈ {jj, kk, kj}, is computed from the simulated EEG signals based on the non-linear model Eq. (1) and the SPM spectral toolbox [13], gn(ν, θ) is the predicted PSD with the correct θ based on the linearized model Eq. (5), gn ν, θ PSDs with θ (under the optimized model structure) using L-DCM, D-DCM and DCM respectively. The coefficients of determination R 2 L-DCM , R 2 D-DCM and R 2 DCM stand for the measures of the goodness of fit between sample and estimated PSDs averaged over 100 trials, where R 2 κ = 1 -ν (gn(ν)-gn(ν, θ) κ ) 2 var(gn(ν)) , κ ∈ {L -DCM, D -DCM, DCM }.

Figure 5 .

 5 Figure 5. Boxplots of the family free energy (F F f , f = 1, 2, 3, 4) values over 100 trials using L-DCM, D-DCM and DCM for the four families (unidirectional influence) in cPBM. (a-c): the three ground truths (see Figure 4 (a-i), (b-i), (c-i)) corresponding to M2, M5 and M6, respectively.

Figure 6 .

 6 Figure 6. (a-i), (b-i), (c-i): Three ground truths M8, M16s and M16as for bidirectional

  Simulated time signals displaying epileptic activities (i.e. (a-ii), (b-ii), (c-ii)), corresponding auto PSD amplitudes (i.e. (a-iii, a-iv),(b-iii, b-iv),(c-iii, c-iv)) and cross PSD amplitudes (i.e. (av),(b-v),(c-v)) for the three ground truths M8, M16s and M16as, respectively. The sample PSD, denoted by gn (ν), n ∈ {jj, kk, kj}, is computed from the simulated EEG signals based on the non-linear model Eq. (1) and the SPM spectral toolbox [13], gn(ν, θ) is the predicted PSD with the correct θ based on the linearized model Eq. (5), gn ν, θ L-DCM , gn ν, θ D-DCM and gn ν, θ DCM represent the estimated PSDs with θ (under the optimized model structure) using L-DCM, D-DCM and DCM respectively. The coefficients of determination R 2 L-DCM , R 2 D-DCM and R 2 DCM stand for the measures of the goodness of fit between sample and estimated PSDs averaged over 100 trials, where R 2 κ = 1 -ν (gn(ν)-gn(ν, θ) κ ) 2 var(gn(ν)) , κ ∈ {L -DCM, D -DCM, DCM }.

  in section Support materials S2. Moreover, in terms of PSD estimation, all methods show again good estimation accuracy which is clearly assessed by the obtained coefficients of determination depicted in Fig.

Figure 7 .

 7 Figure 7. Boxplots of the family free energy (F F f , f = 1, 2, 3, 4) values over 100 trials using L-DCM, D-DCM and DCM for the four families (bidirectional influence) in cPBM.(a-c): the three ground truths (see Figure6 (a-i), (b-i), (c-i)) corresponding to M8, M16s and M16as, respectively.

4. 2 .

 2 Real iEEG signalsOur purpose is to detect brain regions responsible for the seizure onset and effective connectivity consists in inferring the directional flow between brain regions. In this study, iEEG signals were recorded in an epileptic patient in the Centre Hospitalier Universitaire in Rennes. Two particular channels were selected by the clinician: the Pp4 channel, whose activity was confirmed to be part of the seizure onset zone and the Dp1 channel, corresponding to a less epileptogenic region triggered in the seizure. It is worth mentioning that the bipolar signal derived from channel Pp4 (obtained as the difference between Pp4 and Pp5) was associated with the output y j (t) of population pop j whereas the bipolar signal derived from channel Dp1 (obtained as the difference between Dp1 and Dp2) was assigned with the output y k (t) of population pop k . Effective connectivity based on PBM, where the simulations were already investigated in[START_REF] Xiang | A local adjustment strategy for the initialization of dynamic causal modelling to infer effective connectivity in brain epileptic structures[END_REF], and cPBM between these two bipolar signals was evaluated for the first part of the ictal phase (22 ∼ 46s) where channel Pp4 showed a fast onset activity (see Fig.8 (a)). The segment for each scenario was partitioned using 2s length sliding windows without any overlap. Each windowed signal segment of 2s length was normalized and used to calculate the sample PSD, G(ν). A unidirectional connectivity from Pp4 to Dp1 in the ictal phase could be considered as ground truth according to the clinician. Therefore, the model structure M2 in PBM or family F2 (i.e. M2, M5, M6) in cPBM stands for the ground truth based on clinical expertise, as indicated in Fig.8 (b).Regarding the model/family selection, under the cPBM model, L-DCM, D-DCM and DCM methods identify M5 as the optimal model structure underlying the observed iEEG signals as shown in Support materials S3. In terms of the maximized free energy and the family free energy depicted in Fig.9 (a-b), all methods point to a causal effect from pop j to pop k (M2 for Fig.9(a) and F2 for Fig.9 (b)) whatever the NMM (e.g. PBM or cPBM).

Figure 8 .Figure 9 .

 89 Figure 8. (a): Two iEEG signals with a seizure onset up to 24s in channels y j (t) and y k (t) from 22s to 46s in the ictal phase. Each channel corresponds to a bipolar iEEG signal. (b): Prior information provided by the clinical expertise. Results corresponding to PBM (c-i, c-ii, c-iii) and cPBM (d-i, d-ii, d-iii) in terms of the sample and estimated auto/cross PSD amplitudes together with their coefficients of determination R 2 .

  -ii)(d-ii) that all methods show good estimation accuracy under cPBM (R 2 L-DCM = 0.9977, R 2 D-DCM =0.9975 and R 2 DCM =0.9967), compared to PBM (R 2 L-DCM = 0.9623, R 2 D-DCM =0.9565 and R 2 D-DCM =0.9368). Regarding the execution time shown in Table

( 2 ) 12 and 13 .

 21213 The D-DCM algorithm shows a higher model identification rate compared to the DCM and the L-DCM ones in the cases where M2, M6 and M16s are considered as ground truths and also whatever the threshold value. Besides, regarding the case where M8 is the ground truth, the superiority of L-DCM over DCM and D-DCM are also to be noticed whatever the value of φ.Support materialsS1: Bayesian model comparison for unidirectional influenceThese results for Table3can be interpreted in terms of the median (less prone to outliers than the mean) of the maximized free energy depicted in Figs. 11, For instance, regarding the L-DCM and D-DCM algorithms, a relatively higher false model identification rate occurred when the ground truth is M2, as

16 i=1

 16 exp( Fi) where Fi denotes the median of the maximized free energy associated to the model structure Mi. It is clear from Figs. 11, 12 and 13 that the model evidence measure allows for a global quantification of the ability of the different algorithms in identifying the correct model structure except in ground truth M2 where M16 is identified using DCM.

Figure 11 .

 11 Figure 11. Bayesian model comparison in terms of the median of maximized free energy over 100 trials and the model posterior probability using (a-b) L-DCM, (c-d) D-DCM and (e-f ) DCM for the ground truth M2 in the context of unidirectional influence.

Figure 12 .

 12 Figure 12. Bayesian model comparison in terms of the median of maximized free energy over 100 trials and the model posterior probability using (a-b) L-DCM, (c-d) D-DCM and (e-f ) DCM for the ground truth M5 in the context of unidirectional influence.

Figure 14 .

 14 Figure 14. Bayesian model comparison in terms of the median of maximized free energy over 100 trials and the model posterior probability using (a-b) L-DCM, (c-d) D-DCM and (e-f ) DCM for the ground truth M8 in the context of bidirectional influence.

  0, P i is considered as constant.

	Parameter (Pi)	Description (Unit)	Value (Vi)	Prior (p (θi))
		Synaptic time parameter (s)
	ω e,l ω s,l ω f,l	Excitatory Dendritic slow inhibitory Somatic fast inhibitory Intrinsic connection within population 1/75 1/30 1/100	N (0, 1/32) N (0, 1/32) N (0, 1/32)
	C ep ll C pe ll C sp ll C ps ll	Pp → Pe Pe → Pp Pp → Ps Ps → Pp	54 54 54 67.5	N (0, 0) N (0, 0) N (0, 0) N (0, 0)

Table 3 .

 3 The four most plausible model structures and model identification using simulated iEEG signals in the context of unidirectional influence in cPBM, the values in bold indicate the best performance under the true model structure.

	Ground truth Method	The four most plausible model structures (Model identification count over 100 trials)
		L-DCM	M2 (38/100)	M6 (35/100)	M16 (13/100)	M10 (13/100)
	M2	D-DCM	M2 (48/100)	M6 (35/100)	M10 (7/100)	M16 (7/100)
		DCM	M16 (59/100)	M2 (4/100)	M6 (14/100)	M10 (22/100)
		L-DCM	M5 (57/100)	M6 (14/100)	M7 (19/100)	M8 (8/100)
	M5	D-DCM	M5 (71/100)	M6 (9/100)	M7 (17/100)	M8 (0/100)
		DCM	M5 (76/100)	M6 (6/100)	M7 (14/100)	M8 (4/100)
		L-DCM	M6 (87/100)	M2 (0/100)	M8 (12/100)	M10 (0/100)
	M6	D-DCM	M6 (99/100)	M2 (0/100)	M10 (0/100)	M12 (0/100)
		DCM	M6 (97/100)	M2 (0/100)	M10 (0/100)	M16 (0/100)

Table 4 .

 4 Family selection over 100 trials using simulated iEEG signals in the context of unidirectional influence in cPBM, the values in bold indicate the best performance.

			Family identification count over 100 trials
	Ground truth Method	(Averaged run time over 100 trials (mn))
			F1	F2	F3	F4
		L-DCM	0/100 (0.4)	73/100 (1.9)	0/100 (1.0)	27/100 (3.4)
	M2	D-DCM	0/100 (0.2)	83/100 (1.0)	0/100 (0.9)	17/100 (2.9)
		DCM	0/100 (0.1)	18/100 (0.6)	0/100 (0.6)	82/100 (2.0)
		L-DCM	0/100 (0.1)	71/100 (0.9)	0/100 (0.4)	29/100 (2.1)
	M5	D-DCM	0/100 (0.1)	80/100 (0.5)	0/100 (0.6)	20/100 (2.2)
		DCM	0/100 (0.1)	82/100 (0.3)	0/100 (0.3)	18/100 (1.7)
		L-DCM	0/100 (0.2)	87/100 (0.8)	0/100 (0.7)	13/100 (2.3)
	M6	D-DCM	0/100 (0.1)	99/100 (0.5)	0/100 (0.8)	1/100 (2.5)
		DCM	0/100 (0.1)	97/100 (0.4)	0/100 (0.5)	3/100 (1.8)

Table 5 .

 5 The four most plausible model structures and model identification using simulated iEEG signals in the context of bidirectional influence in cPBM. M16 s : symmetrical structure. M16 as : asymmetrical structure, the values in bold indicate the best performance under the true model structure.

	Ground truth Method	The four most plausible model structures (Model identification count over 100 trials)
		L-DCM	M8 (85/100)	M16 (14/100)	M12 (1/100)	M2 (0/100)
	M8	D-DCM	M8 (80/100)	M16 (19/100)	M12 (1/100)	M6 (0/100)
		DCM	M8 (79/100)	M16 (21/100)	M12 (0/100)	M6 (0/100)
		L-DCM	M16 (80/100)	M8 (11/100)	M12 (9/100)	M6 (0/100)
	M16 s	D-DCM	M16 (96/100)	M8 (2/100)	M12 (2/100)	M6 (0/100)
		DCM	M16 (94/100)	M12 (1/100)	M8 (5/100)	M6 (0/100)
	M16 as

Table 6 .

 6 Family selection over 100 trials using simulated iEEG signals in the context of bidirectional influence in cPBM, the values in bold indicate the best performance.

			Family identification count over 100 trials
	Ground truth Method	(Averaged run time over 100 trials (mn))
			F1	F2	F3	F4
		L-DCM	0/100 (0.2)	0/100 (0.7)	0/100 (0.8)	100/100 (2.4)
	M8	D-DCM	0/100 (0.2)	0/100 (0.7)	0/100 (0.8)	100/100 (2.6)
		DCM	0/100 (0.1)	0/100 (0.4)	0/100 (0.5)	100/100 (1.8)
		L-DCM	0/100 (0.3)	0/100 (0.8)	0/100 (0.8)	100/100 (2.8)
	M16 s	D-DCM	0/100 (0.2)	0/100 (0.7)	0/100 (0.7)	100/100 (2.4)
		DCM	0/100 (0.1)	0/100 (0.5)	0/100 (0.5)	100/100 (1.7)
		L-DCM	0/100 (0.1)	0/100 (0.7)	0/100 (0.8)	100/100 (2.4)
	M16 as	D-DCM	0/100 (0.1)	0/100 (0.6)	0/100 (0.6)	100/100 (2.2)
		DCM	0/100 (0.1)	0/100 (0.4)	0/100 (0.5)	100/100 (1.6)

Table 7 .

 7 Averaged execution time for all model structures over all windows (mn) for PBM and cPBM, NMM: neural mass model.

	NMM	Method L-DCM D-DCM DCM
		PBM	10.49	12.25	5.73
		cPBM	3.70	4.75	3.36

Table 9 .

 9 The identification rates of the second plausible model structure (according to Fig.10) using simulated iEEG signals under six ground truths over 100 trials. The values in bold indicate the best performance.

	Ground truth	Method Mm	Model identification count over 100 trials φ = 0 φ = 0.05 φ = 0.1 φ = 0.15
		L-DCM	M6 35/100	33/100	24/100	11/100
	M2	D-DCM M6 35/100 30/100 23/100 10/100
		DCM	M2	4/100	17/100	31/100	37/100
		L-DCM	M6	14/100	11/100	6/100	5/100
	M5	D-DCM M6	9/100	8/100	7/100	4/100
		DCM	M6	6/100	6/100	6/100	5/100
		L-DCM	M2	0/100	0/100	0/100	0/100
	M6	D-DCM M2	0/100	0/100	0/100	0/100
		DCM	M2	0/100	0/100	0/100	0/100
		L-DCM M16 14/100 14/100 14/100 13/100
	M8	D-DCM M16 19/100	19/100	18/100	17/100
		DCM	M16 21/100	21/100	21/100	19/100
		L-DCM	M8	11/100	11/100	11/100	11/100
	M16 s	D-DCM M8	2/100	2/100	2/100	2/100
		DCM	M12	1/100	1/100	1/100	1/100
		L-DCM M10	0/100	0/100	0/100	0/100
	M16 as	D-DCM M10	0/100	0/100	0/100	0/100
		DCM	M10	0/100	0/100	0/100	0/100
	whatever the threshold value.

The identification rate of Mm refers to the degrees of preference of Mm computed over 100 trials. The identifiability of different models was assessed using the Bayesian model comparison at the group level where the fixed effects analysis was conducted[START_REF] Stephan | Bayesian model selection for group studies[END_REF]. More precisely, for each model in a set of plausible sixteen model structures, the median (less prone to outliers than the mean) of the maximized free energy was computed over 100 trials and then used to define the confusion matrices (more details are given in Appendix: A).

Acknowledgments

This work has been supported by Natural Science Foundation of China under Grant (62001240, 31400842, 61876037), the National Key Research and Development Program of China (2017YFB1303200), the National Key R&D Program of China (2017YFC0107900), the Short-term Recruitment Program of Foreign Experts (WQ20163200398), the CRIBs, the Leading-edge Technology and Basic Research Program of Jiangsu Province (BK20192004) and the China

Declaration of Ethical Approval

All experimental procedures have been conducted at the CHU (Centre Hospitalier Universitaire), Rennes, France, following the ethical and regulatory standards. The patient has signed a consent and has been informed that his iEEG data would be used for clinical research and might serve for publication.

Declaration of Competing Interest

The authors declare that they have no conflict of interest. Scholarship Council. The authors thank Jean-Jacques Bellanger, Isabelle Merlet, and Fabrice Wendling for their fruitful discussions and Rennes CHU for providing the database.

Appendix A: confusion matrices

Bayesian model comparison at the group level [START_REF] Penny | Comparing families of dynamic causal models[END_REF] can be also performed using the confusion matrix where the evidence of the sixteen model structures that has been used to invert the different datasets generated using the six model structures (e.g. M2, M5, M6, M8, M16 s and M16 as ) is reported. To this end, three confusion matrices, one per method, were generated and illustrated in Fig. 10. Each entry of these matrices is the median of the maximized free energy computed over 100 trials. From this figure, we observed that the three methods can select the correct model structure except that the DCM method which selected M16 whereas the ground truth was M2. The superscript number indicates the four most plausible model structures for each method and each ground truth as given in Tables 3 and5.

Appendix B: thresholding for removing spurious connections

To deal with potential false detections induced by the small non-null value(s) of the coupling parameter(s), a thresholding based post-processing step is adopted.

This means that, if one/some of the estimated parameter(s) exp( θK pp kj ), exp( θK pp kj ), exp( θK pp kj ) and exp( θK pp kj ) is/are less than a pre-defined threshold parameter φ, it/they is/are set to 0. This can potentially lead to define another model structure and consequently modify its identification rate as shown in Tables 8 and9, where several threshold values (e.g. φ ∈ {0, 0.05, 0.1, 0.15}) are considered. From the above two tables, we can observe that: (1) in the case of unidirectional flow, the model identification rate of the correct model structure is increased when the value of the threshold φ is increased as well regardless the considered method.

It results in a decreased false detection rate. Regarding the case of bidirectional flow, no significant improvement in the model identification rate can be observed

S2: Bayesian model comparison for bidirectional influence

The results for Table 5 can be furthermore interpretable in Figs. vs 5% and 11% respectively (resp. 2% vs 1% and 9% respectively) as displayed in Table 5. Furthermore, as depicted in Fig. 16 where model structure M16 as stands for the ground truth, the gap on the median of the maximized free energy between the optimal model structure M16 and the model structure M10 (the model structure giving the second more important free energy among the 16 model structures in this situation) is significant (higher than 270) whatever the method. The optimal model structure (M16) has a greater model posterior probability than any other model structure, so that 100% model identification rate is obtained by all techniques in this case (see Table 5).


Pierre Cladé

Saïda Guellati-Khélifa

François Nez

François Biraben

Large Momentum Beamsplitter using Bloch Oscillations

Keywords: numbers: 37.25.+k, 03.75.Dg, 67.85.-d

The sensitivity of an inertial sensor based on an atomic interfermometer is proportional to the velocity separation of atoms in the two arms of the interferometer. In this paper we describe how Bloch oscillations can be used to increase this separation and to create a large momentum transfer (LMT) beamsplitter. We experimentally demonstrate a separation of 10 recoil velocities. Light shifts during the acceleration introduce phase fluctuations which can reduce the contrast of the interferometer. We precisely calculate this effect and demonstrate that it can be significantly reduced by using a suitable combination of LMT pulses. We finally show that this method seems to be very promising to realize LMT beamsplitter with several 10s of recoil and a very good efficiency.

Over the past two decades, the impressive advances on the control and the manipulation of atomic de Broglie waves by coherent light pulses had led to the emergence of highly sensitive atomic interferometers [1,2,3,4,5]. Atom interferometry is based on coherent splitting and recombination of atomic wave packets. The sensitivity of the atomic interferometer, which is proportional to the spatial separation of the wave packets, can be enhanced by increasing the interaction time. This time is limited in some experiments: this is the case for gravity measurement where the free fall of atoms is a limitation as well as recoil velocity measurement. In our recent measurement of the recoil velocity [6], this time was limited by the size of the vacuum cell. Indeed atoms accelerated to velocities close to 1600 recoils (≃ 10 m/s for 87 Rb atoms) move of several centimeter during the measurement. Furthermore, vibrations or phase noise of lasers which are more important for longer interaction time are also limitations for the duration of the measurement.

The velocity splitting is usually achieved by using coherent momentum exchange between atoms and light (Raman or Bragg transition) [START_REF] Bergman | Atom interferometry[END_REF]: with such a process atoms are prepared in a superposition of two states with a well control relative velocity. In the case of Raman transitions and first order Bragg transitions (the most widely used configurations), this velocity difference δv is fixed to 2v r where v r is the recoil velocity. To circumvent the limitation due to the time measurement, one can enhance the interferometer sensitivity by increasing the velocity separation between its two arms. Several attempts have been investigated to increase this separation [START_REF] Gupta | [END_REF]9,10,11,12,13].

In this paper, we investigate a method to realize the momentum transfer using the so called "Bloch oscillation" (BO) method [14,15]. In this scheme, the large momentum transfer beam splitter (LMTBS) is obtained in a two steps process: the first one is a regular beamsplitter which transfers 2 recoil velocities while the second one is a large momentum transfer pulse (LMT pulse), that uses BO to coherently accelerate atoms of one of the two wave packets. When the beam splitter is used to recombine wave packets, the order of the two steps is reversed. This method was initially suggested in Ref. [11], where a proof of principle experiment is described.

Compare to other methods, BO allow to transfer a large number of recoils with a high efficiency. This method seems therefore to be the most suitable for high precision measurements. In this letter we focus on one of the main systematic effect which is due to level shift of atoms in the lattice and produces a large phase shift in the interferometer. This phase shift depends on the position of atoms in the laser beam and the motion of atoms induces fluctuations which reduce the contrast of the interferometer. It is especially large in our experiment in which, compared to Ref. [11], we use a non degenerate gas. We calculate this effect and demonstrate that it can be significantly reduced by using a suitable combination of LMT pulses.

The principle of Bloch oscillations [15] consists in shining on the atoms two counterpropagating laser beams of similar frequencies (ν and ν +δν). The interference of the two lasers results in an optical lattice of depth U 0 whose velocity is proportional to δν. When this lattice is accelerated, and under specific conditions, atoms are coherently accelerated. Because of the periodic potential, the eigenenergies of the system exhibit a band structure [16], each eigen state being described by a quasimomentum q defined modulo 2 k (where k = 2π/λ is the laser wave vector) and a band index n. The figure 1 represents this band structure, where we have unfold the first Brillouin zone. There is a gap between bands either at the center (q = 0) or at the edge of the first Brillouin zone (q = ± k) when atoms are resonantly coupled. When an atom in the first band is subjected to a constant and uniform force, its quasimomentum increases lineary with time. If the force is weak enough to avoid adiabatic transition, the atom stays in the first band and therefore has a periodic motion (Bloch oscillations [START_REF] Zener | Proc. Soc. Lond. A[END_REF]). The oscillations can occur also in higher band. However, the bandgap decreases with higher band index. Therefore, if an atom passes through the crossing at a given speed, the probability to make an adiabatic transition (to stay in the same band) will be higher for low value of the band index.

The principle of the large momentum transfer beam splitter consists in creating a superposition of two wavepackets separated by 2 recoil velocities and in loading them in the optical lattice such that one wavepacket is in the first band (A, see Fig. 1) and the second in the third band (B). A constant acceleration is then applied to the lattice. This acceleration (which acts like a force in the frame of the lattice) brings atoms to the edge of the first Brillouin zone. The acceleration of the lattice is small enough such that atoms in the first band have a large probability to make an adiabatic transition but high enough so that atoms in the third band change band. Each oscillation increases the momentum of the atom by 2 k. On the other hand, atoms in the third band (atoms that change band) are not accelerated.

For the efficiency of the LMT beamsplitter, the acceleration of the lattice should be weak enough to adiabatically accelerate atoms in the first band and strong enough to induce non-adiabatic transition for atoms in the third band. Reciprocally, for a given acceleration, the amplitude of the lattice should be in a given range. The figure 2 depicts the probability for an atom to stay in its band as a function of the lattice amplitude for the first (solid line) and third (dashed line) band. There is clearly an intermediate regime where the probability η 11 for an atom to stay in the first band is high whereas the probability η 34 for an atom to leave the third band (and reach the fourth one) is also high. For an acceleration of 2 recoils in 200 µs, the total probability (η = η 11 η 34 , diamond) presents a maximum around U 0 = 8E r . The value of the maximum (97%) depends on the duration of the acceleration and increases with this parameter.

We have plotted on the right side of Fig. 2 the total efficiency η as a function of the initial momentum p 0 . This efficiency is computed including the loading and unloading of atoms in the lattice: it is initially ramped up during a time t adiab , then accelerated during T acc and ramped down during t adiab . By ramping adiabatically up and down the lattice, atoms from plane wave states are transferd to Bloch states and vice versa. As this process is not fully adiabatic, the efficiency of the LMT is reduced. At the center and the edge of the frist Brillouin zone (q 0 = 0 and q 0 = 1), the efficiency is strongly reduced because atoms cannot be loaded adiabatically in the lattice (those points are initally degenerate). For the choosen parameters (t adiab = 150 µs, T acc = 200 µs, N = 2 oscillations, U 0 = 8E r ), we see that the efficiency is larger than 95% on a large zone. An important issue is to maximize the width in initial momentum where the process is very efficient. Indeed atoms used in the interferometer have an initial velocity distribution selected by the Raman beam. The wider is the initial velocity distribution loaded into the LMT pulse, the higher is the number of atoms that contributes to the interferometer and so is the signal to noise ratio. We have optimized the efficiency by varying the amplitude and the temporal parameters keeping the total time 2t adiab + T acc = 500 µs constant.

One of the main drawback of the LMTBS based on Bloch oscillations is the light shift of atoms in the lattice. In the case of a blue detuned lattice, atoms in the first band are in a dark region and are almost not shifted while non-accelerated atoms in excited band see an average shift coresponding to the mean value of the potential of the lattice. For typical parameters, this light shift, much larger than 2π, must be cancelled in order to run the interferometer. This cancellation occurs in the Mach-Zender configuration described on Fig. 3. The configuration used for the Raman pulses is similar to a regular interferometer with four π/2-pulses and the LMT pulses are added inside each pair of π/2 pulses used either for selection or measurement. With this scheme, the LMT pulses are applied symetrically on each arm of the interferometer, i.e. one arm of the interferometer is initially in the first band and then in an excited band and vice versa for the other arm. Therefore, the phase shift accumulated on each arm is the same and there is no systematic effect if the laser intensity seen by atoms is constant. However, this is not the case because of temporal fluctuations of the laser intensity (leading to a phase noise in the interference pattern) or motion of atoms through the spatial profile of the laser beam (leading to a reduction of the fringes contrast). This last effect is not negligeable in our experiment: indeed with our experimental parameters, the typical relative variation of the laser intensity due to the motion of atoms in the gaussian profile of the laser is 6 s -1 (meaning a variation of 6% in 10 ms). We numerically estimate that the sensitivity of the interferometer to such a linear temporal variation of the intensity is about 500 rad.s, leading to a typical phase shift of ∆φ = 3 rad. Consequelty, there is a complete supression of fringes.

In order to run the inteferometer, we used a sequence (see Fig. 4) with 8 LMT pulses instead 4. After (or before) each π/2 pulses, we accelerate successively each arm of the interferometer with two LMT pulses of opposite directions. With this scheme, the second LMT pulse compensates 90% of the light shift [19]. The sensitivity of the interferometer to a variation of the intensity is then reduced by a factor of 10 compared to the 4 LMT pulses interferometer. The resulting typical phase shift is ∆φ = 0.3 rad, which does not affect significantly the contrast of fringes. In the experiment, there is a correlation between the velocity and the position of atoms due to the time of flight before the begining of the interferometer. This results in a biased variation of the intensity which leads to a systematic effect in the interferometer. This effect can be canceled by inverting the order of the two LMT pulses for each beamsplitter.

The experimental setup is similar to the one described in [START_REF] Cladé | [END_REF]. We use 87 Rb atoms cooled using a magnetooptical trap followed by an optical molasses. The Raman beams, in a vertical configuration, are produced by two phase locked laser diodes and an acousto optic modulator (AOM) is used to control the shape of the π/2 pulses. The Raman beams are brought to the science chamber using an optical fiber. In order to reduce phase noise, a single polarization maintaining fiber is used. One of the two beams, selected by its polarization, is retroreflected in order to realize a counterpropagating transition.

The Bloch beams are produced by a Ti:Sa laser. A new scheme has been implemented in order to reduce the phase of the Bloch beams, which affect the interferometer. As for the Raman beams, a single fiber and a retroreflecting mirror are used to create the optical lattice. But in contrary to the Raman beams, both beams used for Bloch are retroreflected. By sending to the AOM two frequencies (ν and ν +δν) we create four lattices, one at the velocity of the atoms and three other lattices (one with the opposite velocity and two at rest). In the experiment, we wait about 10 ms between the end of the molasses and the beginning of the interferometer. Because of the gravity, the velocity of the atoms is then about 20 recoils and the three other lattices are sufficiently out of resonance not to disturb the interferometer. The amplitude and frequencies sent to the AOM are controled using an RF chain based on an arbitrary waveform generator (similar to the one described in [START_REF] Cladé | [END_REF]). In order to measure the velocity change of atoms between the first part and second part of the interferometer, the frequency of both the Raman and Bloch beams are changed in an identical way. The frequency of the central peak corresponds then to the Doppler effect induced by the variation of velocity.

The left side of Fig. 5 shows the proportion of atoms detected in the internal state F = 1 as a function of the frequency difference between selection and measurement (see the temporal sequence of Fig. 4). The fringes sepration is δ sep = 50 Hz for this sequence. It corresponds to an effective time T eff = 1/δ sep of about 20 ms. We can compare this periodicity with the periodicity of fringes of an interferometer without LMT Bloch pulses. With the delay T Ramsey of 5 ms between the begining of the two Ramsey pulses, the periodicity would be 200 Hz. There On the right side of Fig. 5 we have plotted T eff for three different durations of the selection: keeping the same π/2-Bloch-Bloch sequence, we change T Ramsey . In the inset of the figure, we have plotted the gain in resolution (ratio between the effective time and T Ramsey ) as a function of T Ramsey . At the limit where T Ramsey is long, the duration of the LMT Bloch pulses can be neglected and the gain is exactly 2N + 1, i.e. 5 for our parameters. At shorter time, this gain is smaller. The predicted effective time is plotted on Fig 5 (solid line). The good adjustement between the predicted time and our experimental data is a strong indication that the oscillations observed are due to the LMT pulses.

The amplitude of the signal is roughly 10 times smaller than the amplitude without LMT pulses. Different sources seems to contributes to these losses of amplitude. Our model does not take into account the non resonnant lattices due to the retroreflecting configuration. Also, as explained earlier in the paper (Fig. 2), we have to take into acount the initial distribution in quasi-momentum. We calculate that the total efficiency for 8 pulses (in the optimal configuration) is about 50%. This is the maximal efficiency and higher losses are actually expected because the intensity of the lattice is not well controlled and also is not constant over the atomic cloud. Spontaneous emission, which for our parameters is of the order of 4% for each of the LMT pulse, leads also to an additional reduction of 30% of the signal. The use a laser beam with a wider waist, a higher intensity (and a larger detuning) should allow us to reduce significantly those effects.

In order to further increase the velocity separation, one can imaging to have a two step process: the first step is identical to the one described earlier, with two Bloch oscillations. For the second step, we can increase adiabatically the depth of the lattice and then use a larger acceleration. A long time (and a relatively weaker lattice) is needed in order to realize the first separation with a good efficiency. However, for atom in the higher band, the Landau-Zener tunneling is larger and one can have a good efficiency with a shorter time and a deeper lattice. With an initial separation obtained with 2 oscillations in 250 µs with U 0 = 6E r , a second separation with 8 oscillations in 150 µs (with U 0 = 24E r ) and a total duration (taking into acount the adiabatic ramps) of 1 ms, the calculated efficiency is close to the 4 recoils LMT pulse (the losses are mainly due to the initial splitting). This scheme could then be used on both arms of the interferometer to realize a separation of 42 recoils.

In this paper, we have theoretically and experimentaly studied the implementation a large momentum beamsplitter in an atomic interferometer. It is possible to realize a separation of 4 recoils between the 2 arms of the interfermeter within 500 µs. The main issue to take care of is the light shift which degrades strongly the contrast of the interferometer. Our proposal allows us to reduce this effect by a factor of 10 and finally permits us to realize an interferometer with a separation of 10 recoils between the two arms. A gain of 4 in the resolution compared to a usual interferometer has been observed. This method seems to be very promizing for the realization of a beam splitter with separation of several tens of recoil velocity. The effect due to light shifts has been reduced by accelerating successivelly both arm of the interferometer. A way to suppress it systematically would be to accelerate both arms on the interfereometer simultaneously (instead of doing it successively) by applying two counter-propagating accelerated lattices. In this case light shift would be completely supressed. In the futur, we plan to combine this kind of LMT beam splitter with Bloch oscillations between the two pairs of π/2 pulses [6] to realize a new measurement of the recoil velocity and therefore improve the determination of the fine structure constant.

FIG. 1 :

 1 FIG. 1: Bande structure of the optical lattice. Trajectories of the accelerated (A) and non accelerated atoms (B).

FIG. 2 :

 2 FIG. 2: Left: transfer probability as a function of the maximal optical depth U0 of the lattice. U0 is in unit of 8Er. The acceleration is of 2 recoils in 200 µs. Solid line (blue): transfer probability for the first band η11; dashed line (red): for the third band (η33 ≈ 1 -η34); diamond: Efficiency of the LMT pulse, η = η11η34. Right: efficiency of the LMT pulse as a function of the initial momentum p0.

FIG. 3 :

 3 FIG.3: Temporal sequence of the interferometer. Bottom: intensity of the Raman (dashed line) and Bloch (solid line) beam; Top: trajectories of atoms in the two arms; dashed and solid lines correspond to the two internal states of the atom.

FIG. 4 :

 4 FIG. 4: Temporal sequence of the 4 LMT interferometer. Bottom: intensity of the Raman (dashed line) and Bloch (solid line) beam; Top: trajectories of atoms; dashed and solid line correspond to the two internal states of the atom.

FIG. 5 :

 5 FIG. 5: Left: Interferences fringes observed using the scheme depicted Fig. 4. Right: Effective time T eff as a function of the delay TRamsey between the two π/2-pulses; solid curve: numerical calculation; Diamond: experimental results