Emilie Blanc
email: eblanc@lma.cnrs-mrs.fr

Guillaume Chiavassa
email: gchiavassa@centrale-marseille.fr

Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

A time-domain numerical method for Biot-JKD poroelastic waves in 2D heterogeneous media

Keywords: numbers: 43.20.-Gp porous media, elastic waves, Biot-JKD model, fractional derivatives, numerical modeling, diffusive representation, finite difference methods, Cartesian grid

An explicit finite-difference scheme is presented for solving the two-dimensional Biot equations of poroelasticity across the full range of frequencies. The key difficulty is to discretize the Johnson-Koplik-Dashen (JKD) model which describes the viscous dissipations in the pores. Indeed, the timedomain version of Biot-JKD model involves order 1/2 shifted fractional derivatives which amounts to a time convolution product. To avoid storing the past values of the solution, a diffusive representation of fractional derivatives is used: the convolution kernel is replaced by a finite number of memory variables that satisfy local-in-time ordinary differential equations. The coefficients of the diffusive representation follow from an optimization procedure of the dispersion relation. Then, various methods of scientific computing are applied: the propagative part of the equations is discretized using a fourth-order ADER scheme, whereas the diffusive part is solved exactly. An immersed interface method is implemented to discretize the geometry on a Cartesian grid, and also to enforce the jump conditions at interfaces. Numerical experiments are proposed in various realistic configurations.

I. INTRODUCTION

Two frequency regimes have to be distinguished when dealing with poroelastic waves. In the low-frequency range (LF), [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I: Low-frequency range[END_REF] the viscous efforts are proportional to the relative velocity of the fluid and the solid matrix. In the high-frequency range (HF), modeling the dissipation is a more delicate task: Biot first presented an expression for particular pore geometries. [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II: High-frequency range[END_REF] In 1987, Johnson-Koplik-Dashen (JKD) [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF] published a general expression for the dissipation in the case of random pores. The viscous efforts depend in this model on the square root of the frequency. In the time domain, time fractional derivatives are then introduced, which involves convolution products with singular kernels. Transient analytical solutions of Biot-JKD model have been derived in academic geometries and homogeneous media. [START_REF] Fellah | Ultrasonic wave propagation in human cancellous bone: application of Biot theory[END_REF] When dealing with more complex geometries, numerical methods are required.

Many time-domain simulation methods have been developed in the LF regime: see Ref. 5 and the introduction to Ref. 6 for general reviews. In the HF regime, the fractional derivatives greatly complicate the numerical modeling of the Biot-JKD equations. The past values of the solution are indeed required in order to evaluate the convolution products, which means that the time evolution of the solution must be stored. It greatly increases the memory requirements and makes large-scale simulations impossible. To our knowledge, two approaches have been proposed so far in the literature. The first approach consists in discretizing the convolution products, [START_REF] Masson | Finite-difference modeling of Biot's poroelastic equations across all frequencies[END_REF] and the second one is based on the use of a diffusive representation of the fractional derivative. [START_REF] Lu | Wave field simulation for heterogeneous porous media with singular memory drag force[END_REF] In the latter approach, the convolution product is replaced by a continuum of diffusive variables -or memory variables -satisfying local differential equations. [START_REF] Haddar | Efficient solution of a wave equation with fractional-order dissipative terms[END_REF] This continuum is then discretized using appropriate quadrature formulas, resulting in the Biot-DA (diffusive approximation) model.

In Ref. 10, the diffusive approach was followed in a one-dimensional homogeneous configuration. Compared with Ref. 8, important improvements were introduced: good representation of the viscous dissipation on the whole range of frequencies; optimization of the model on the range of interest; estimation of the computational effort in terms of the required accuracy. The goal of the present work is to extend the algorithm to two-dimensional heterogeneous porous media with interfaces.

This paper is organized as follows. The original Biot-JKD model is briefly outlined in section II and the principles underlying the diffusive representation of fractional derivatives are described. The decrease of energy and the dispersion analysis are addressed. In section III, the method used to discretize the diffusive model is presented: the diffusive approximation thus obtained is easily treatable by computers. Following a similar approach than in viscoelasticity, [START_REF] Gautier | Propagation of acoustic waves in a one-dimensional macroscopically inhomogeneous poroelastic material[END_REF] the coefficients of the model are determined using an optimization procedure in the frequency range of interest, giving an optimum number of additional computational arrays. The numerical modeling is addressed in section IV, where the equations of evolution are split into two parts: a propagative part, which is discretized using a fourth-order scheme for hyperbolic equations, and a diffusive part, which is solved exactly. An immersed interface method accounts for the jump conditions and for the geometry of the interfaces, preventing from the usual limitations of finite differences on a Cartesian grid. Numerical experiments performed with realistic values of the physical parameters are presented in section V. In section VI, a conclusion is drawn and some futures lines of research are given.

II. PHYSICAL MODELING

A. Biot model

The Biot model describes the propagation of mechanical waves in a macroscopic porous medium consisting of a solid matrix saturated with a fluid circulating freely through the pores. [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I: Low-frequency range[END_REF][START_REF] Bourbié | Acoustics of Porous Media[END_REF][START_REF] Carcione | Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media[END_REF] It is assumed that

• the wavelengths are large in comparison with the diameter of the pores;

• the amplitude of the perturbations is small;

• the elastic and isotropic matrix is completely saturated with a single fluid phase;

• the thermo-mechanical effects are neglected. This model involves 10 physical parameters: the density ρ f and the dynamic viscosity η of the fluid; the density ρ s and the shear modulus µ of the elastic skeleton; the porosity 0 < φ < 1, the tortuosity a ≥ 1, the absolute permeability κ, the Lamé coefficient λ f and the two Biot's coefficients β and m of the saturated matrix. The following notations are introduced

ρ w = a φ ρ f , ρ = φ ρ f + (1 -φ) ρ s , χ = ρ ρ w -ρ 2 f > 0, λ 0 = λ f -m β 2 . ( 1
)
Taking u s and u f to denote the solid and fluid displacements, the unknowns in 2D are the elastic velocity

v s = ∂ us ∂ t , the filtration velocity w = ∂ W ∂ t = φ ∂ ∂ t (u f -u s )
, the elastic stress tensor σ, and the acoustic pressure p. The constitutive laws are

σ = (λ f tr ε -β m ξ) I + 2 µε, p = m (-β tr ε + ξ), (2)
where ε = 1 2 (∇u s + ∇u s T ) is the strain tensor and ξ = -∇.W is the rate of fluid change. The symmetry of σ implies compatibility conditions between spatial derivatives of ε, leading to the Beltrami-Michell equation: 14

∂ 2 σ xy ∂x ∂y = θ 0 ∂ 2 σ xx ∂x 2 + θ 1 ∂ 2 σ yy ∂x 2 + θ 2 ∂ 2 p ∂x 2 + θ 1 ∂ 2 σ xx ∂y 2 + θ 0 ∂ 2 σ yy ∂y 2 + θ 2 ∂ 2 p ∂y 2 , θ 0 = - λ 0 4(λ 0 + µ) , θ 1 = λ 0 + 2 µ 4(λ 0 + µ) , θ 2 = µ β 2(λ 0 + µ) . (3)
On the other hand, the conservation of momentum yields

       ρ ∂v s ∂t + ρ f ∂w ∂t = ∇σ, ρ s ∂v s ∂t + ρ w ∂w ∂t + η κ F * w = -∇p, (4)
where * is the convolution product in time. The second equation of ( 4) is a generalized Darcy law. The quantity F * w denotes the viscous dissipation induced by the relative motion between the fluid and the elastic skeleton.

B. High frequency dissipation: the JKD model

The frontier between the low-frequency range (LF) and the high-frequency range (HF) is reached when the viscous efforts are similar to the inertial effects. The transition frequency is given by [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I: Low-frequency range[END_REF][START_REF] Bourbié | Acoustics of Porous Media[END_REF]

f c = η φ 2 π a κ ρ f = ω c 2 π . (5)
In LF, the flow in the pores is of the Poiseuille type, and dissipation efforts in the second equation of (2) are given by

F (t) = δ(t) ⇐⇒ F (t) * w(x, y, t) = w(x, y, t), ( 6
)
where δ is the Dirac distribution. In HF, the width of the viscous boundary-layer is small in comparison with the size of the pores, and modeling the dissipation process is a more complex task. Here we adopt the widely-used model proposed by Johnson-Koplik-Dashen (JKD) in 1987, which is valid for random networks of pores with constant radii. [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF] The only additional parameter is the viscous characteristic length Λ. We take

P = 4 a κ φ Λ 2 , Ω = ω c P = η φ 2 Λ 2 4 a 2 κ 2 ρ f , ( 7
)
where P is the Pride number (typically P ≈ 1/2). Based on the Fourier transform in time, F (ω) = R F (t)e -iωt dt, the frequency correction given by the JKD model can be written

F (ω) = 1 + i ω 4 a 2 κ 2 ρ f η Λ 2 φ 2 1/2 , = 1 + i P ω ω c 1/2 , = 1 √ Ω (Ω + i ω) 1/2 . ( 8
)
This correction is the simplest function satisfying the LF and HF limits of the dynamic permeability. [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF] Therefore, the term F (t) * w(x, y, t) involved in the second equation of (2) is

F (t) * w(x, y, t) = F -1 1 √ Ω (Ω + i ω) 1/2 w(x, y, ω) , = 1 √ Ω (D + Ω) 1/2 w(x, y, t). (9)
The operator D 1/2 is a shifted order 1/2 time fractional derivative, generalizing the usual derivative characterized by

∂ w ∂ t = F -1 (i ω w(x, y, ω)).
The notation (D + Ω) 1/2 accounts for the shift Ω in (9).

C. The Biot-JKD equations of evolution

Based on (2), ( 4) and ( 9), the Biot-JKD equations can be written

                   ρ ∂ v s ∂ t + ρ f ∂ w ∂ t = ∇σ, ρ f ∂ v s ∂ t + ρ w ∂ w ∂ t + η κ 1 √ Ω (D + Ω) 1/2 w = -∇p, σ = (λ f tr ε -β m ξ)I + 2 µ ε, p = m (-β tr ε + ξ). (10)
We rearrange this system by separating ∂ vs ∂ t and ∂ w ∂ t in the first two equations of (10) and by using the definitions of ε and ξ. Taking

γ = η κ ρ χ 1 √ Ω , (11)
Biot-JKD poroelastic waves in 2D heterogeneous media 3 one obtains the following system of equations of evolution

                                                                           ∂ v sx ∂ t - ρ w χ ∂ σ xx ∂ x + ∂ σ xy ∂ y - ρ f χ ∂ p ∂ x = ρ f ρ γ (D + Ω) 1/2 w x + f vsx , ∂ v sy ∂ t - ρ w χ ∂ σ xy ∂ x + ∂ σ yy ∂ y - ρ f χ ∂ p ∂ y = ρ f ρ γ (D + Ω) 1/2 w y + f vsy , ∂ w x ∂ t + ρ f χ ∂ σ xx ∂ x + ∂ σ xy ∂ y + ρ χ ∂ p ∂ x = -γ (D + Ω) 1/2 w x + f wx , ∂ w y ∂ t + ρ f χ ∂ σ xy ∂ x + ∂ σ yy ∂ y + ρ χ ∂ p ∂ y = -γ (D + Ω) 1/2 w y + f wy , ∂ σ xx ∂ t -(λ f + 2µ) ∂ v sx ∂ x -m β ∂ w x ∂ x -λ f ∂ v sy ∂ y -m β ∂ w y ∂ y = f σxx , ∂ σ xy ∂ t -µ ∂ v sy ∂ x + ∂ v sx ∂ y = f σxy , ∂ σ yy ∂ t -λ f ∂ v sx ∂ x -m β ∂ w x ∂ x -(λ f + 2µ) ∂ v sy ∂ y -m β ∂ w y ∂ y = f σyy , ∂ p ∂ t + m β ∂ v sx ∂ x + ∂ w x ∂ x + β ∂ v sy ∂ y + ∂ w y ∂ y = f p . (12)
Terms f vsx , f vsy , f wx , f wy , f σxx , f σxy , f σyy and f p have been to simulate sources.

D. The diffusive representation

Taking

D Ω w(x, y, t) = ∂ w ∂ t + Ω w, (13)
the shifted fractional derivative (9) can be written 15

(D + Ω) 1/2 w(x, y, t) = 1 √ π t 0 e -Ω(t-τ ) √ t -τ D Ω w(x, y, τ ) dτ. (14)
The operator (D + Ω) 1/2 is not local in time and involves the entire time history of w. As we will see in section III, a different way of writing this derivative is more convenient for numerical evaluation. Based on Euler's Γ function, the diffusive representation of the totally monotone function

1 √ t is 9 1 √ t = 1 √ π ∞ 0 1 √ θ e -θt dθ. (15)
Substituting ( 15) into ( 14) gives

(D + Ω) 1/2 w(x, y, t) = 1 π t 0 ∞ 0 1 √ θ e -θ(t-τ ) e -Ω(t-τ ) D Ω w(x, y, τ ) dθ dτ, = 1 π ∞ 0 1 √ θ ψ(x, y, θ, t) dθ, (16)
where the diffusive variable is defined as

ψ(x, y, θ, t) = t 0 e -(θ+Ω)(t-τ ) D Ω w(x, y, τ ) dτ. (17)
Biot-JKD poroelastic waves in 2D heterogeneous media 4

For the sake of clarity, the dependence on Ω and w is omitted in ψ. From (17), it follows that the diffusive variable ψ satisfies the ordinary differential equation

   ∂ ψ ∂ t = -(θ + Ω) ψ + D Ω w, ψ(x, y, θ, 0) = 0. (18)
The diffusive representation therefore transforms a non-local problem ( 14) into a continuum of local problems (18). It should be emphasized at this point that no approximation have been made up to now. The computational advantages of the diffusive representation will be seen in sections III and V, where the discretization of ( 16) and ( 18) will yield a numerically tractable formulation.

E. Energy of Biot-JKD

Now, we express the energy of the Biot-JKD model (10). This result generalizes the analysis performed in the 1D case in the Ref. 10.

Proposition 1 Setting C the 4 × 4 poroelastic tensor such that σ = C ε -β p I in (2), let E = E 1 + E 2 + E 3 , with E 1 = 1 2 R 2 ρ v 2 s + ρ w w 2 + 2 ρ f v s .w dx dy, E 2 = 1 2 R 2 C ε : ε + 1 m p 2 dx dy, E 3 = 1 2 R 2 θ∈R + η κ 1 π 1 √ Ω θ 1 θ + 2 Ω (w -ψ) 2 dθ dx dy. (19)
Without any source terms, E is an energy which satisfies

dE dt = - R 2 θ∈R + η κ 1 π 1 √ Ω θ 1 θ + 2 Ω Ω w 2 + (θ + Ω) ψ 2 dθ dx dy ≤ 0. ( 20
)
Proposition 1 calls for the following comments:

• the Biot-JKD model is well-posed;

• when the viscosity of the saturating fluid is neglected (η = 0), the energy of the system is conserved;

• the terms in ( 19) have a clearly physical significance: E 1 is the kinetic energy, and E 2 is the potential energy. The term E 3 corresponds to the kinetic energy resulting from the relative motion of the fluid inside the pores.

F. Dispersion analysis

A plane wave de i(ωt-k.r) is injected in (12), where k = k e and d are the wavevector and the polarization, respectively; r is the position, ω = 2 π f is the angular frequency and f is the frequency. If d is collinear with k, the dispersion relation of compressional waves is obtained. Setting

             D 4 = m (λ 0 + 2 µ), D 2 (ω) = -((λ f + 2 µ) ρ w + m (ρ -2 ρ f β)) ω 2 + i ω η κ F (ω) (λ f + 2 µ), D 0 (ω) = χ ω 4 -i ω 3 η κ ρ F (ω), (21)
Biot-JKD poroelastic waves in 2D heterogeneous media 5

the dispersion relation of compressional waves takes the form

D e (k, ω) = D 4 k 4 + D 2 (ω) k 2 + D 0 (ω) = 0. ( 22
)
If d is orthogonal with k, the dispersion relation of the shear wave is obtained. Setting

             D 4 (ω) = ω 2 (ρ + φ ρ f (a -2)) -i ω φ 2 η κ F (ω), D 2 (ω) = -ω 2 φ ρ f (a -1) + i ω φ 2 η κ F (ω), D 0 (ω) = ω 2 φ ρ f a -i ω φ 2 η κ F (ω), (23)
the dispersion relation of shear wave takes the form

k = 1 √ µ D 4 D 0 -D 2 2 D 0 1/2 . ( 24
)
Expressions ( 21)-( 22)-( 23)-( 24) are valid in the case of both the Biot-LF and Biot-JKD models with the frequency correction defined by

F (ω) =      F LF (ω) = 1 Biot-LF, F JKD (ω) = 1 √ Ω (Ω + i ω) 1/2 Biot-JKD. ( 25
)
The solutions k pf , k ps of ( 22) and the solution k s of ( 24) give the phase velocities c pf = ω/ℜe(k pf ) of the fast compressional wave, c ps = ω/ℜe(k ps ) of the slow compressional wave, and c s = ω/ℜe(k s ) of the shear wave, with 0 < c ps < c pf and 0 < c s . The attenuations α pf = -ℑm(k pf ), α ps = -ℑm(k ps ) and α s = -ℑm(k s ) can also be deduced. Both the phase velocities and the attenuations of Biot-LF and Biot-JKD are strictly increasing functions of the frequency. The high frequency limits of phase velocities of compressional waves, c ∞ pf and c ∞ ps , obtained by diagonalizing the left-hand side of system (12), satisfy the relation

χ c 4 -((λ f + 2 µ) ρ w + m (ρ -2 ρ f β)) c 2 + m (λ 0 + 2 µ) = 0, (26)
and the high frequency limit of phase velocity of the shear wave, c ∞ s , is

c ∞ s = √ µ ρ - φ ρ f a -1/2 . ( 27
)
Figure 1 shows the dispersion curves corresponding to the Biot-LF and Biot-JKD models. The physical parameters are those used in the numerical experiments presented in section V (medium Ω 0 ). Note that the scales are radically different in the case of fast and slow waves. The following properties can be observed:

• when f < f c , the Biot-JKD and Biot-LF dispersion curves are very similar as might be expected, since

lim ω→0 F JKD (ω) = 1;
• the fast compressional wave and the shear wave are almost not affected by the frequency correction F (ω), while the slow compressional wave is greatly modified;

• when f ≪ f c , the slow compressional wave is almost static. When f > f c , the slow wave propagates but is greatly attenuated.

III. THE BIOT-DA (DIFFUSIVE APPROXIMATION) MODEL

The aim of this section is to approximate the Biot-JKD model, using a numerically tractable approach.

Biot-JKD poroelastic waves in 2D heterogeneous media 6

phase velocity of the attenuation of the fast compressional wave fast compressional wave Biot-JKD poroelastic waves in 2D heterogeneous media 7

A. The Biot-DA first-order system Using a quadrature formula on N points, with weights a ℓ and abscissa θ ℓ > 0, the diffusive representation ( 16) can be approximated by

(D + Ω) 1/2 w(x, y, t) = 1 π ∞ 0 1 √ θ ψ(x, y, t, θ) dθ, ≃ N ℓ=1 a ℓ ψ(x, y, t, θ ℓ ), ≡ N ℓ=1 a ℓ ψ ℓ (x, y, t). (28)
>From (18), the N diffusive variables ψ ℓ satisfy the ordinary differential equations

   ∂ ψ ℓ ∂ t = -(θ ℓ + Ω) ψ ℓ + D Ω w, ψ ℓ (x, y, 0) = 0. (29)
The fractional derivatives are replaced by their diffusive approximation (28) in the JKD model (12). Upon adding the equations (29) and performing some straightforward operations, the Biot-DA system is written as a first-order system in time and space

                                                                                                         ∂ v sx ∂ t - ρ w χ ∂ σ xx ∂ x + ∂ σ xy ∂ y - ρ f χ ∂ p ∂ x = ρ f ρ γ N ℓ=1 a ℓ ψ xℓ + f vsx , ∂ v sy ∂ t - ρ w χ ∂ σ xy ∂ x + ∂ σ yy ∂ y - ρ f χ ∂ p ∂ y = ρ f ρ γ N ℓ=1 a ℓ ψ yℓ + f vsy , ∂ w x ∂ t + ρ f χ ∂ σ xx ∂ x + ∂ σ xy ∂ y + ρ χ ∂ p ∂ x = -γ N ℓ=1 a ℓ ψ xℓ + f wx , ∂ w y ∂ t + ρ f χ ∂ σ xy ∂ x + ∂ σ yy ∂ y + ρ χ ∂ p ∂ y = -γ N ℓ=1 a ℓ ψ yℓ + f wy , ∂ σ xx ∂ t -(λ f + 2µ) ∂ v sx ∂ x -m β ∂ w x ∂ x -λ f ∂ v sy ∂ y -m β ∂ w y ∂ y = f σxx , ∂ σ xy ∂ t -µ ∂ v sy ∂ x + ∂ v sx ∂ y = f σxy , ∂ σ yy ∂ t -λ f ∂ v sx ∂ x -m β ∂ w x ∂ x -(λ f + 2µ) ∂ v sy ∂ y -m β ∂ w y ∂ y = f σyy , ∂ p ∂ t + m β ∂ v sx ∂ x + ∂ w x ∂ x + β ∂ v sy ∂ y + ∂ w y ∂ y = f p , ∂ ψ xj ∂ t + ρ f χ ∂ σ xx ∂ x + ∂ σ xy ∂ y + ρ χ ∂ p ∂ x = Ω w x -γ N ℓ=1 a ℓ ψ xℓ -(θ j + Ω) ψ xj + f wx , ∂ ψ yj ∂ t + ρ f χ ∂ σ xy ∂ x + ∂ σ yy ∂ y + ρ χ ∂ p ∂ y = Ω w y -γ N ℓ=1 a ℓ ψ yℓ -(θ j + Ω) ψ yj + f wy . (30)
where j = 1, ..., N . Taking the vector of unknowns U = (v sx , v sy , w x , w y , σ xx , σ xy , σ yy , p, ψ x1 , ψ y1 , . . . , ψ xN , ψ yN ) T (31)

and the source vector

F = (f vxs , f vys , f wx , f wy , f σxx , f σxy , f σyy , f p , f wx , f wy , . . . , f wx , f wy ) T , (32)
Biot-JKD poroelastic waves in 2D heterogeneous media 8

the system (30) can be written

∂ U ∂ t + A ∂ U ∂ x + B ∂ U ∂ y = -S U + F, ( 33
)
where A is the (2 N + 8) 2 propagation matrix along x, B is the (2 N + 8) 2 propagation matrix along y and S is the (2 N + 8) 2 dissipation matrix (A1, A2, A3). The size of the system increases linearly with the number N of diffusive variables.

B. Properties

Four properties of system (33) are specified:

• the eigenvalues of A (A1) and B (A2) are real: 0 with multiplicity 2 N + 2, ±c ∞ pf , ±c ∞ ps and ±c ∞ s . The system (33) is therefore hyperbolic;

• since the eigenvalues and eigenvectors do not depend on the diffusive coefficients, they are the same in both the Biot-DA and Biot-LF or Biot-JKD models. This is not so in the case of the method presented in the Ref. 7, where the propagation matrix is modified to account for the fractional derivative;

• the dispersion analysis is obtained in the case of the Biot-DA model by replacing F by

F AD (ω) = Ω + i ω √ Ω N ℓ=1 a ℓ θ ℓ + Ω + i ω (34)
in equations ( 21)- (22).

C. Determining the Biot-DA parameters

The coefficients a ℓ and θ ℓ in (28) have to be determined. For this purpose, we want to approach the original frequency correction F JKD (ω) (25) by F AD (ω) (34) in a given frequency range of interest. Let Q(ω) be the optimized quantity and Q ref (ω) be the desired quantity

       Q(ω) = F AD (ω) F JKD (ω) = N ℓ=1 a ℓ (Ω + i ω) 1/2 θ ℓ + Ω + i ω = N ℓ=1 a ℓ q ℓ (ω), Q ref (ω) = 1. (35)
We implement a linear optimization procedure [START_REF] Blanc | Biot-JKD model: simulation of 1D transient poroelastic waves with fractional derivatives[END_REF][START_REF] Emmerich | Incorporation of attenuation into time-domain computations of seismic wave fields[END_REF][START_REF] Groby | A time domain method for modeling viscoacoustic wave propagation[END_REF][START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF] in order to minimize the distance between Q(ω) and Q ref (ω) in the interval [ω min , ω max ] centered on ω 0 = 2 π f 0 , where f 0 is the central frequency of the source. The abscissas θ ℓ are fixed and distributed linearly on a logarithmic scale

θ ℓ = ω min ω max ω min ℓ-1 N -1 , ℓ = 1, ..., N. (36)
The weights a ℓ are obtained by solving the system

N ℓ=1 a ℓ q ℓ (ω k ) = 1, k = 1, ..., K, (37)
where the ωk are also distributed linearly on a logarithmic scale of K points

ωk = ω min ω max ω min k-1 K-1 , k = 1, ..., K. (38)
Since the q ℓ (ω) are complex functions, optimization is performed simultaneously on the real and imaginary parts

             N ℓ=1 a ℓ Re(q ℓ (ω k )) = 1, N ℓ=1 a ℓ Im(q ℓ (ω k )) = 0, k = 1, ..., K. (39)
For practical purposes, we use ω min = ω 0 /10 and ω max = 10 ω 0 , and we solve an overdetermined system (K = N ), which can be solved by writing normal equations. [START_REF] Flannery | Numerical Recipes in C: the Art of Scientific Computing[END_REF] The number of diffusive variables N can be determined in terms of the ratio f 0 /f c and of the chosen accuracy: see Ref. 10 for details.

IV. NUMERICAL MODELING A. Splitting

In order to integrate the Biot-DA system (33), a uniform grid is introduced, with mesh size ∆ x = ∆ y and time step ∆ t. The approximation of the exact solution U(

x i = i ∆ x, y j = j ∆ y, t n = n ∆ t) is denoted by U n ij .
If an unsplit integration of ( 33) is performed, a Von-Neumann analysis typically yields the stability condition

∆ t ≤ min Υ ∆ x c ∞ pf , 2 R(S) , ( 40
)
where R(S) is the spectral radius of S, and Υ > 0 depends on the numerical scheme. We have no theoretical estimate of R(S), but numerical studies have shown that this value is similar to that of the spectral radius in LF: η κ ρ χ , which can be very large. [START_REF] Chiavassa | Time domain numerical modeling of wave propagation in 2D heterogeneous porous media[END_REF] The time step can therefore be highly penalized in this case (40).

A more efficient strategy is adopted here, which consists in splitting the original system (33) into a propagative part and a diffusive part:

       ∂ U ∂ t + A ∂ U ∂ x + B ∂ U ∂ y = 0, ∂ U ∂ t = -S U. ( 41
)
For the sake of simplicity, the source term F has been omitted here. The discrete operators associated with the propagative part and the diffusive part are denoted by H a and H b , respectively. The second-order Strang splitting is then used to integrate (33) between t n ant t n+1 , giving the time-marching

• U (1) ij = H b ( ∆ t 2 ) U n ij , • U (2) ij = H a (∆ t) U (1) ij , • U n+1 ij = H b ( ∆ t 2 ) U (2) ij . (42)
The discrete operator H a associated with the propagative part (first equation of ( 41)) is an ADER 4 (Arbitrary DERivatives) scheme. [START_REF] Schwartzkopff | Fast high order ADER schemes for linear hyperbolic equations[END_REF] This scheme is fourth-order accurate in space and time, is dispersive of order 4 and dissipative of order 6, and has a stability limit Υ = 1. On Cartesian grids, ADER 4 amounts to a fourth-order Lax-Wendroff scheme.

Since the physical parameters do not vary with time, the diffusive part (second equation of 41)) can be solved exactly. This gives

H b ∆ t 2 U ij = e -∆ t 2 S U ij . (43)
The matrix e -∆ t 2 S is computed numerically using the (6, 6) Padé approximation in the "scaling and squaring method". [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF] The numerical integration of the diffusive step (43) is unconditionally stable. [START_REF] Blanc | Biot-JKD model: simulation of 1D transient poroelastic waves with fractional derivatives[END_REF] The full algorithm (42) is therefore stable under the optimum CFL (Courant-Friedrichs-Lewy) stability condition

CFL = c ∞ pf ∆ t ∆ x Υ, (44)
which is always independent of the Biot-DA model coefficients.

Biot-JKD poroelastic waves in 2D heterogeneous media 10 Let us consider two homogeneous poroelastic media Ω 0 and Ω 1 separated by a stationary interface Γ, as shown in figure 2. The governing equations (30) in each medium have to be completed by a set of jump conditions. The simple case of perfect bonding and perfect hydraulic contact along Γ is considered here, modeled by the jump conditions: 22

[v s ] = 0, [w.n] = 0, [σ.n] = 0, [p] = 0. ( 45
)
The discretization of the jump conditions requires special care, for three reasons. First, if the interfaces do not coincide with the uniform meshing, then geometrical errors will occur. Secondly, the jump conditions (45) will not be enforced numerically by the finite-difference scheme, and the numerical solution will therefore not converge towards the exact solution. Lastly, the smoothness of the solution required to solve (33) will not be satisfied across the interface, which will decrease the convergence rate of the ADER scheme.

To overcome these drawbacks without detracting the efficiency of Cartesian grid methods, an immersed interface method 6,23 is used. The latter studies can be consulted for a detailed description of this method. The basic principle is as follows: at the irregular nodes where the ADER scheme crosses an interface, modified values of the solution are used on the other side of the interface instead of the usual numerical values.

Calculating these modified values is a complex task involving high-order derivation of jump conditions (45), Beltrami-Michell relations (3) and singular value decompositions. Fortunately, all these time consuming procedures can be carried out during a preprocessing stage and only small matrix-vector multiplications need to be performed during the simulation. After optimizing the code, the extra CPU cost can be practically negligible, i.e. lower than 1% of that required by the time-marching procedure.

The immersed interface method uses only the propagative part (first equation of (41)) and the jump conditions (45), whose do not depend on the frequency regime. Consequently, no changes are required in the high-frequency range compared to the low-frequency range. The algorithm is the same as in the low-frequency range, see Ref. 6 for details.

V. NUMERICAL EXPERIMENTS

General configuration

The physical parameters used in all the numerical experiments, which are given in table I, correspond to Cold Lake sandstone (Ω 0 ) and Berea sandstone (Ω 1 ) saturated with water. The computational domain [-0.1, 0.1] 2 m is discretized with N x × N y grid points, N x = N y = 1600, which amounts to 29 points per slow compressional wavelength in Ω 0 and in Ω 1 for a source of central frequency f 0 = 200 kHz. As seen in 1D, 10 N = 6 diffusive variables are used to ensure an error of model smaller than 5.58 %. The time step is deduced from (44), taking CFL = 0.95. The numerical experiments are performed on an Intel Core i7 processor at 2.80 GHz.

Parameters

Ω0 Ω1

Saturating fluid ρ f (kg/m 3 ) 1040 1000 η (Pa.s) 1.5 10

Test 1: homogeneous medium

The unbounded homogeneous medium Ω 0 is excited by a point source. The only non-null component of the vector of source F (32) is f σxy = g(x, y) h(t), where g(x, y) is a truncated gaussian centered at point (0, 0), of radius R 0 = 3.79 10 -3 m and Σ = 1.90 10 -3 m:

g(x, y) =      1 π Σ 2 exp - x 2 + y 2 Σ 2 if 0 x 2 + y 2 R 2 0 , 0 otherwise, (46)
and h(t) is a Ricker signal of central frequency f 0 = 200 kHz and of time-shift t 0 = 2/f 0 = 10 -5 s:

h(t) =      2 π 2 f 2 0 t - 1 f 0 2 -1 exp -π 2 f 2 0 (t - 1 f 0 ) 2 if 0 t t 0 , 0 otherwise. (47)
We use a truncated gaussian for g(x, y) rather than a Dirac to avoid spurious numerical artifacts localized around the source point. This source generates cylindrical waves of all types: fast and slow compressional waves and shear wave, which are denoted P f , P s and S, respectively, in figure 3. Fast and slow compressional waves are observed as regards the pressure, while the additional shear wave is present in the σ yy component of the stress tensor. No special care is applied to simulate outgoing waves (with PML, for instance), since the simulation is stopped before the waves have reached the edges of the computational domain.

Test 2: diffraction of a plane wave by a plane interface In all the following tests, the source is a plane right-going fast compressional wave, whose wavevector k makes an angle θ = 0 with the horizontal x-axis. Its time evolution is the same Ricker signal as in the first test (47). We use periodic boundary conditions at the top and at the bottom of the domain. This test is conducted on a vertical plane interface at x = 0 m. The incident P f -wave (Ip f ) propagates in the medium Ω 0 . The incident wave crosses the interface normally, leading to a 1D configuration. Consequently, only the fast and slow compressional waves propagate. From a numerical point of view, however, the problem is fully bidimensional. The advantage of such a 1D configuration is that each diffracted wave has interacted with the interface and is consequently very sensitive to the discretization of the interface conditions (45). The figure 4 shows a snapshot of the pressure at t 1 = 1.82 10 -5 s, on the whole computational domain. The reflected fast and slow compressional waves, denoted Biot-JKD poroelastic waves in 2D heterogeneous media 12 In this case, we compute the exact solution of Biot-DA thanks to standard tools of Fourier analysis. Special care is taken to ensure that the initial data and the exact solution are valuable reference solutions : N f = 512 Fourier modes are used, with a frequency step ∆ f = 4000 Hz. The figure 5 shows the excellent agreement between the analytical and the numerical values of the pressure along the line y = 0 m.

Test 3: diffraction of a plane wave by a cylinder or a circular shell

In the first case, a cylindrical scatterer Ω 1 of radius 0.015 m is centered at point (0.01, 0). In the second case, we consider a poroelastic shell Ω 1 of external radius 0.03 m and internal radius 0.015 m, centered at point (0.025, 0). The Biot-JKD poroelastic waves in 2D heterogeneous media 13 Test 4: multiple ellipsoidal scatterers In the last test, the ability of the proposed numerical strategy to handle complex geometries with variable curvature is illustrated. 15 ellipsoidal scatterers, of major radius 0.008 m and of minor radius 0.005 m, of medium Ω 1 are randomly embedded in Ω 0 . Performing such a simulation with hundreds of scatterers has physical applications. Doing so provide a mean to check the validity of multiple scattering models. [START_REF] Luppé | Coherent waves in a multiply scattering poro-elastic medium obeying Biot's theory[END_REF][START_REF] Chekroun | Time-domain numerical simulations of multiple scattering to extract elastic effective wavenumbers[END_REF] Biot-JKD poroelastic waves in 2D heterogeneous media 14

VI. CONCLUSION

An explicit finite-difference method is presented here for simulating transient poroelastic waves in the full range of validity of poroelasticity. The Biot-JKD model, which involves order 1/2 fractional derivatives, was replaced here by an approximate Biot-DA model, which is much more tractable numerically. The Biot-DA coefficients are determined here using an optimization procedure, which depends on the frequency range of interest. The hyperbolic system of partial differential equations was discretized using efficient tools (Strang splitting, fourth-order ADER scheme, immersed interface method). This numerical methodology provides highly accurate numerical solutions, allowing a fine numerical investigation of realistic porous media. Some future lines of research are suggested:

• Continuously variable porous media. In the previous numerical experiments, the physical parameter characterizing the porous media are constant and discontinuous across the interfaces. The possibility of applying the presented numerical method to 2D porous media with continuously variable coefficients, [START_REF] Gautier | Propagation of acoustic waves in a one-dimensional macroscopically inhomogeneous poroelastic material[END_REF] where no analytical expressions are available, is currently in progress: see Ref. 10 in 1D.

• Multiple scattering. The numerical tools presented here make possible the modeling of multiple scattering in random media. Based on simulated data, the properties of the effective medium amount to the disordered under study can be deduced. [START_REF] Luppé | Coherent waves in a multiply scattering poro-elastic medium obeying Biot's theory[END_REF][START_REF] Chekroun | Time-domain numerical simulations of multiple scattering to extract elastic effective wavenumbers[END_REF] Optimization and parallelization of the method must be realized to take into account of hundreds of inclusions in our simulations.

• Anisotropy. The anisotropy of some porous media is sometimes a main feature, for instance in biomechanics to model trabecular and cortical bones. [START_REF] Bossy | Three-dimensional simulations of ultrasonic axial transmission velocity measurement on cortical bone models[END_REF][START_REF] Lemoine | High-resolution finite volume modeling of wave propagation in orthotropic poroelastic media[END_REF] The anisotropy of a medium changes only the 4 × 4 poroelastic tensor C (1). So the ADER 4 scheme and the immersed interface method have to be modified, but the approximation of the fractionnal derivative is still valid.

• Thermic boundary-layer. In cases where the saturating fluid is a gas, thermo-mechanical effects have to be taken into account. Extended versions of the Biot-JKD have been developed, [START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF] involving additional order 1/2 fractional derivatives. The numerical method developed in this paper should lend itself well to working with this model.

Figure 1 .

 1 Figure 1. Dispersion curves: comparison between Biot-LF and Biot-JKD.

B. Immersed interface methodFigure 2 .

 2 Figure 2. Interface Γ between two poroelastic media Ω0 and Ω1.

Figure 3 .

 3 Figure 3. Test 1. Fast and slow compressional waves, respectively P f and Ps, and shear wave S emitted by a source point at (0.01, 0). Pressure (left) and σyy component of the stress tensor (right) at t1 = 3.98 10 -5 s.

Figure 4 .

 4 Figure 4. Test 2. Snapshot of p at initial time (a) and at t1 = 1.82 10 -5 s (b).

Pressure

  Figure 5. Test 2. Pressure along the line y = 0 m; vertical line denote the interface. Comparison between the numerical values (circle) and the analytical values (solid line) of p at t1 = 1.82 10 -5 s.

Figure 6 .

 6 Figure 6. Test 3. Scatterer: cylinder. Snapshot of σyy at initial time (a) and at t1 = 3.63 10 -5 s (b).

Figure 7 .

 7 Figure 7. Test 3. Scatterer: circular shell. Snapshot of σyy at initial time (a) and at t1 = 3.63 10 -5 s (b).

Figure 8 .

 8 Figure 8. Test 4. Snapshot of σyy at initial time (a) and at t1 = 4.03 10 -5 s (b).

  Appendix A: MATRICES OF PROPAGATION AND DISSIPATIONThe matrices in (33) are λ f + 2 µ) 0 -m β 0

Table I .

 I Physical parameters used in numerical experiments.

			-3	10 -3
	Grain	ρs (kg/m 3 )	2650	2644
		µ (Pa)	2.93 10 9	7.04 10 9
	Matrix	φ	0.335	0.2
		a	2	2.4
		κ (m 2 )	10 -11	3.6 10 -13
		λ f (Pa)	6.14 10 9	1.06 10 10
		m (Pa)	6.49 10 9	9.70 10 9
		β	9.56 10 -1	7.20 10 -1
		Λ (m)	2.19 10 -5	5.88 10 -6
	Phase velocities c ∞ pf (m/s)	2384.17	3269.89
		c ∞ ps (m/s)	758.95	814.95
		c ∞ s (m/s)	1229.00	1776.16
		fc (Hz)	3.84 10 3	3.68 10 4


Bilger Alexandre

Essert Caroline

Duriez Christian

Cotin Stéphane

Brain-shift aware risk map for Deep Brain Stimulation Planning

In Deep Brain Stimulation surgery, the eciency of the procedure heavily relies on the accuracy of the placement of the stimulating electrode. Meanwhile, the eectiveness of the placement is dicult due to brain shifts occurring during and after the procedure. We propose an approach to overcome the limitations of current planning software that ignores brain shift. In particular, we consider the motion of vascular structures in order to reduce risks of dissecting a vessel during the procedure. Facing the diculty to produce an exact brain shift prediction, we propose to build a brain shift aware risk map which embeds the vascular motion risk. This risk map is extrapolated using simulation from clinical studies that provide statistics on the displacement of anatomical landmarks during the procedure. Risk maps can be directly integrated into automatic path planning algorithms to better predict optimal electrode trajectories. The method relies on a physics-based simulation that takes into account brain deformation, electrode placement, cerebrospinal uid, and vascular motion. The goal is to reproduce the spread of brain shift situations that are noted in clinical studies. Preliminary results show that it is possible to compute safe electrode trajectories even in case of brain shift and yet optimal regarding the placement within the targeted area.

Introduction

Over the past decade, Deep Brain Stimulation (DBS) has known a growing interest in neurology for the treatment of movement disorders such as Parkinson's disease or dystonia. This surgical procedure consists in implanting an electrode in deeply located structures of the brain, among which the SubThalamic Nucleus (STN). An accurate placement of the electrode is crucial to maximize outcomes and to prevent adverse eects. This placement is achieved in two main stages: rst, pre-operative medical images of the patient are combined with the use of a stereotactic frame (and sometimes an atlas of the brain) to determine the target coordinates and optimal trajectory for the electrode(s). Second, the patient is taken to the operating room where a hole is drilled in the skull to insert electrodes into the brain according to the planned trajectory (two in the case of a bilateral implantation). However, the opening of the skull changes the intracranial pressure which causes a combination of deformation and motion of the brain known as brain shift. The trajectory is in green, from an insertion point on the scalp to the target in red. The area in black is the pneumocephalus. The white circle represents a vessel. In (a) the vessel lies in its initial position. In (b) and (c), we see its new position according to the brain shift, and the whole shape of its movement during the progression of the brain shift (light dotted shape).

The brain shift model remaining unknown, the deformation can not be anticipated. Fig. 1 represents dierent possible intensities of brain shift that might occur with the same insertion point. When this brain shift is signicant, it can have an impact on the procedure and its outcome [START_REF] Miyagi | Brain shift: an error factor during implantation of deep brain stimulation electrodes[END_REF] [START_REF] Elias | Cortical and subcortical brain shift during stereotactic procedures[END_REF]. The three main known issues are the following. The brain shift results in a dierence of location of anatomical structures between the planning stage and the intra-operative set-up (blood vessels, ventricles, STN, Anterior and Posterior Commissures (AC/PC)). This dierence means that 1) the eective location for the electrode can be quite remote from the planned location (5 mm or more in rare cases) [START_REF] Miyagi | Brain shift: an error factor during implantation of deep brain stimulation electrodes[END_REF], and 2) there is a risk that an anatomical structure (such as vessels) shifts to the area of the burr-hole. Therefore, it becomes impossible to implant the electrode without risking intracranial hemorrhage. In their survey, Benabib et al. [START_REF] Benabid | Deep brain stimulation of the subthalamic nucleus for the treatment of parkinson's disease[END_REF] report that hemorrhages occurred in 8.4% of all the DBS patients, mostly at the entry point. Moreover, several weeks after the surgery a post-operative electrode migration can appear as the brain returns to its initial position when the subdural air introduced during surgery has resolved. As a result, 3) the electrode might no longer be in contact with the subthalamic area, and brain disorders might resume. Our main objective is to improve the trajectory planning part of the DBS procedure. In this paper, we propose to address problem #2 mentioned above, by presenting a preoperative approach that can reduce risks of dissecting a vessel due to brain shift.

In clinical routine, most neurosurgeons use commercial software such as Medtronic StealthStation or GE Advantage Windows. Such software provide a signicant assistance in determining a good trajectory for the electrode. Usually the software proposes a default solution based on standard angles relatively to AC/PC landmarks or the stereotactic frame, without considering the surrounding anatomy, and lets the neurosurgeon interactively move the virtual electrode to a location which is assumed to be optimal in terms of targeting and avoidance of certain structures (mainly ventricles and blood vessels). As this task can be tedious and time consuming (from 10 to 50 minutes depending on the diculty of the case), a few research groups have proposed to automate the search for an optimal trajectory [3] [11] [6] [START_REF] Bériault | Automatic trajectory planning of DBS neurosurgery from multi-modal MRI datasets[END_REF]. Brunenberg et al. [START_REF] Brunenberg | Automatic trajectory planning for DBS: A feasibility study[END_REF] extract from a set of entry points some valid trajectories, but without ordering them. Their work, as well as those of Shamir et al. [START_REF] Shamir | A method for planning safe trajectories in image-guided keyhole neurosurgery[END_REF] consider only risk criteria, that are taken into account separately. Essert et al. [START_REF] Essert | Automatic computation of electrode trajectories for deep brain stimulation: a hybrid symbolic and numerical approach[END_REF] proposed a method based on a weighted combination of several types of criteria, applied to preoperative MRI. The method proposed by Bériault et al. [START_REF] Bériault | Automatic trajectory planning of DBS neurosurgery from multi-modal MRI datasets[END_REF] is very similar, but is applied on multi-modal MRI. This results in an accurate segmentation of anatomical structures of interest, however only a risk criteria is considered.

Yet, none of the previous planning solutions (manual or automatic) take into account brain deformation or the associated motion of its substructures: the computations are based on static pre-operative images. This is potentially risky as there is no guarantee that vital structures will not have moved before the insertion of the electrode, because of the opening of the skull, nor that the planned trajectory will actually be aiming at the subthalamic area, as it may also have shifted. The objective of our method is to consider brain shift in the preoperative trajectory planning. To this end, we propose the computation of an advanced brain shift aware risk map based not only on static pre-operative images, but on an estimation of the possible location of the vessels due to brain shift, taking into account dierent possible scenarios of deformations. This information can then be used in combination with any planning strategy, whether it is manual or automated.

To achieve this goal, we propose a framework which takes into account intracranial uid loss, subdural air invasion, brain shift and vascular motion, and can rely on patient-specic data. This work includes the following features: physicsbased models of the brain; mechanical interactions between the brain and the skull; inuence of the surrounding cerebro-spinal uid (CSF) and air invasion in the skull on the brain shift (including the asymmetry of the brain shift). We also consider the orientation of the patient on the operating table as it inuences CSF loss and brain shift. Then we couple the results of the simulation with an automatic optimal trajectory planning system. Results from our simulations and computations (performed on generic dataset), show that 1) important brain shifts result in a signicant vascular motion, 2) we can take this motion into account in the planning process with reduced computational load and no signicant changes in the software and 3) we can compute automatically optimal trajectories based on this information, and the result is not too restrictive and still allows for solutions to be found. This paper is organized as follows: Section 2 describes the simulation, and planning processes and how they are linked; Section 3 presents our preliminary results, and nally Section 4 concludes and addresses future steps for the method.

Materials & Methods

Our approach relies on a physics-based model of the brain tissue deformations, the contact response with the skull and the falx cerebri, and the interaction with the cerebro-spinal uid (CSF). We perform numerical simulations to produce a brain shift aware risk map embedding the motion of the vascular structures as well as a measure of uncertainty about their exact location. The map is then used in a planning system to determine automatically optimal trajectories for the electrode. The approaches we used are presented in the following section.

Brain Shift Simulation

Previous works have studied the brain shift with intraoperative images. [START_REF] Ferrant | Registration of 3-d intraoperative mr images of the brain using a nite-element biomechanical model[END_REF] use intraoperative MRI to register them but it requires a very costly equipment, whereas [START_REF] Letteboer | Brain shift estimation in image-guided neurosurgery using 3-d ultrasound[END_REF] estimate the brain shift with 3D-ultrasound. However, these intraoperative techniques assume the brain shift already occurred, therefore the possible trajectories are limited by the hole in the skull. Other works try to preoperatively predict the brain deformation in order to increase the electrode implantation accuracy, but these methods are based on parameters that are unknown at the planning step. We propose a method that helps surgeons to visualize a risk of brain shift using a patient specic simulation before the operation, during the planning. More specically, our method indicates which zones should be avoided because of a high risk of the presence of vessels if a brain shift occurs, with dierent possible levels of intensity.

Numerical Simulation In brain deformation simulation research, many use a nite element method (FEM) associated with linear or non-linear material model, or eventually a uid model. In our simulation we use a non-linear geometric FEM model with a linear constitutive law [START_REF] Felippa | A unied formulation of small-strain corotational nite elements: I. theory[END_REF]: our model is adapted for large displacements and small deformations, which is appropriate in our case. The volume of each brain hemisphere is meshed as a set of tetrahedral elements. We also accurately compute the contacts between the brain, the inner part of the skull and the falx cerebri using static and dynamic friction law models [START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF] which reproduce, at a macroscopic level, the eect of the connective tissue in the subarachnoid space. We also take CSF inuence into account, which is acting on the brain and balancing the brain weight. The consequence of a CSF loss is a brain shift as the brain is more inuenced by gravity than by uid forces. The CSF action is modelled as external forces acting on the surface of the brain:

f CSF = S ρgh(P )dS
where ρ is the density of CSF, g is the norm of the gravity and h is the distance between a point P on the surface and the uid level. This force is computed on each triangle S of the brain mesh that corresponds to the immersed surface.

Independently of the choice of the deformation model, we end up with the following dierential system of non-linear equations

Ma = f(x, v) + p + f CSF + H T λ
where M is the mass matrix, f gathers the internal forces. a, v and x are respectively the acceleration, the velocity and the position of the nodes from the mesh. The forces p are exerted by the gravity and f CSF by the CSF. Finally, H T λ gathers constraints response resulting from unilateral contacts and bilateral constraints.

The vascular motion is associated to the brain motion, that is why we connect the vessels to the brain by keeping constant the barycentric coordinates of the vessels positions in the elements of the brain mesh.

To ensure stability and accuracy, we use an implicit integration scheme, and a GPU implementation guarantees a fast computation for a clinical use.

Construction of the Brain-Shift Aware Risk Map Two phenomena have to be taken into account: 1) the rst one, illustrated on Fig. 1, is the unknown exact intensity of the brain shift that will occur to the patient for a particular location of the burr hole. This is due to the complexity of the model of brain shift and the large number of parameters involved. For this reason, we chose to take into account all possible intensities of brain shift from no brain shift to a maximum brain shift for a particular insertion point; 2) the second one, illustrated on Fig. 2, is the dependency between the maximum possible brain shift, the location chosen for the burr hole, and the orientation of the head compared to gravity. For this reason, we chose to browse the space of possible entry points and simulate the maximum brain shift for the dierent locations.

To address the rst phenomenon, as there is no reliable patient specic brain shift model, we can not consider only the nal position of the vessels in the worst case of brain shift. We don't know what will be the intensity of the brain shift, so we need to consider every possible position of the vessels for every possible intensity of brain shift. That is why we build, from the shapes and positions of the 3D meshes of the vascular system simulated for dierent intensities of brain shift, a larger 3D surface embedding all the possible positions of the vessels. Figs. 3a and 3b illustrate a representative example of how is built such a convoluted surface, whereas the simulated vascular network is depicted in Fig. 3c.

Addressing the second phenomenon implies to simulate the brain shift for several possible locations of the burr hole. The process presented in the previous paragraph is applied to a set of representative locations of the burr hole. This allows us to build an advanced risk map, where the risk associated to an insertion point is the distance between a trajectory and the vascular convoluted surface corresponding to this particular insertion point. By doing this, we forbid a trajectory to cross all intermediate positions of the vessels between its original and a maximal positions, to take into account every possibility of brain shift from the best to the worst.

With this method, a surgeon is able to avoid possibly moving vessels even without an accurate prediction about the exact brain shift or the CSF loss. In the following we show how we integrate this risk map in a planning software, in order to anticipate the movement of these structures. Once the brain shift aware risk map is produced, it can be easily used in any standard DBS planning software. The software we are using for this study has an additional feature: it includes a trajectory optimization process. This process is performed in two phases, thanks to a surgical constraints solver in the spirit of [START_REF] Essert | Automatic computation of electrode trajectories for deep brain stimulation: a hybrid symbolic and numerical approach[END_REF] or [START_REF] Bériault | Automatic trajectory planning of DBS neurosurgery from multi-modal MRI datasets[END_REF]. The rst phase solves hard constraints, such as critical anatomical structures avoidance, and eliminates from an initial outer surface all insertion points impossible to consider in practice without injuring the patient. The second phase solves soft constraints, that correspond to preferences. They are described under the form of cost functions to minimize, aggregated into one linear combination of weighted constraints. The result of these soft constraints are displayed as color maps representing for each candidate insertion point the degree of satisfaction.

We dened a set of soft constraints, concerning the minimization of the risks for the patient and the optimization of the orientation in relation to the shape of the target. Among those soft constraints is one favoring a maximization of the distance between the candidate electrode and risky structures (which are in the case of DBS mostly vessels and ventricles). In this study we had to split this soft constraint in two separate constraints. The rst one simply computes the distance between the electrode and the ventricles. The second one is adapted to use the advanced risk map described in section 2.1, so that the distance with the vessels takes into account the estimated brain shift.

As a consequence, instead of maximizing the distance from the candidate electrode to the vessels extracted from the pre-operative MRI, i.e. in a static position, we maximize the distance from the candidate electrode to the set of positions where the vessels could lie if a brain shift occurs. Thanks to the advanced risk map, the estimated possible positions of the vessels depend on the location of the burr-hole involved by the position of the considered candidate electrode.

Results

This section presents the preliminary results of our method applied on a template of high-delity anatomical three-dimensional models of the brain. This contains brain tissues, skull, skin, ventricles, vessels and the target (STN). In the following, we will show how we use the simulation to produce a brain shift aware risk map, then we compute an optimal trajectory. In our tests, we consider the patient lying in the supine position. The CSF level is expressed in percentage from no CSF (0%) to full CSF volume (100%).

Independence of Mechanical Parameters

First, we show that our simulation reproduces the correlation between the volume of CSF lost and the anterior commissure (AC) shift, and the posterior commissure (PC) shift, found in [START_REF] Elias | Cortical and subcortical brain shift during stereotactic procedures[END_REF] on a group of 66 patients. As there is a Obviously, the Young's modulus has an inuence on the brain deformation. However, Fig. 4 shows that a specic brain deformation based on AC/PC displacement can be obtained whatever the Young's modulus. To reach the same deformation, the larger the Young's modulus is, the more the CSF has to leak out. That is why our approach is relatively independent of the mechanical parameters of the physical model. The model enables to get a valid geometrical deformation. As we could not anticipate the CSF loss, we can now parameterize our simulation with experimental measure like the AC/PC displacement.

Integration in Planning Software

We compared the brain shift aware risk map with the regular distance map on the template. The result can be seen on Fig. 5, which contains snapshots from the planning software. Two color maps corresponding to soft constraints are shown: (a) with regular distance computation (static vessels), and (b) the brain-shift aware risk map computed by the simulation. We can see that the green areas, corresponding to safe insertion points, are narrower and less numerous using our risk map. Only the areas that would be safe even with a brain shift are kept. However, sucient number and surface of green areas are still present. On Fig. 6, the nal result of the automatic optimal trajectory planning, including several constraints other than distance to the vessels, is shown for the template. We can notice that even with the restriction of the safe areas due to our risk map, an optimal trajectory can still be found within the range of the average angles relatively to AC/PC that are used as a basis in clinical routine (around 60 • in antero-posterior and 30 • latero-median axis).

Conclusion

In this paper we described a physics-based method for simulating deformations of cerebral structures and landmarks caused by the brain shift during deep brain stimulation surgery, according to a position of the burr hole and an estimation of the CSF loss. More particularly, we estimated the movement of the vessels which need to be avoided during the insertion of the electrode. We emphasize that our method does not pretend to provide an exact estimation of the brain deformation or vessel motion at the pre-operative planning step, as it is not possible to anticipate the exact CSF loss before the surgery, or other physiological parameters. To compensate for some uncertainties in the simulation (such as exact CSF loss), and to account for other possible errors (such as segmentation of the vessels and other structures on the patient data), we proposed to embed both vessel motion and estimated error into an advanced distance map. This map was used in a path planning software to produce optimal electrode placements. The results of our experiments show the benets of such a simulation, as it does not restrict too much the possible insertion areas but provides safer trajectories regarding the possibilities of brain shift. In the future, we will continue improving the simulation and planning methods. In particular we plan to better determine the inuence of the CSF loss (as done in [START_REF] Elias | Cortical and subcortical brain shift during stereotactic procedures[END_REF]) through additional studies. We also plan to rely on intra-operative images to assess our prediction for the brain shift and vascular motion, and use it to adjust the planned trajectory just after skull opening.

Fig. 1 :

 1 Fig.1: Sagittal projection of the brain for dierent possibilities of intensities of brain shift given an insertion point, from no brain shift (a) to a maximal brain shift (c). The trajectory is in green, from an insertion point on the scalp to the target in red. The area in black is the pneumocephalus. The white circle represents a vessel. In (a) the vessel lies in its initial position. In (b) and (c), we see its new position according to the brain shift, and the whole shape of its movement during the progression of the brain shift (light dotted shape).

Fig. 2 :

 2 Fig.2: Sagittal projection of the brain for 3 dierent trajectories. We assume that the brain shift depends on the location of the insertion point and on gravity. We depict here a patient in supine position. An anterior trajectory (a) causes a smaller brain shift than more posterior trajectories (b) and (c). Trajectories are in green and the target is in red. The area in black illustrates the maximum pneumocephalus that could occur in the worst case.

  Fig. 3

3

 3 Fig. 3

2 . 2 Automatic

 22 Computation of Electrode Tra jectories

Fig. 4 :

 4 Fig. 4: Shifts of AC/PC vs CSF level with three Young's modulus (E). Units of shifts are in mm. AC displacement in blue. PC displacement in orange.

Fig. 5 :

 5 Fig. 5: Color maps representing the risk in relation to the proximity of the vessels, computed without/with taking into account a possible brain shift. In green the safest zones, in red the zones to avoid.

Fig. 6 :

 6 Fig. 6: Optimal trajectory (red cylinder): trajectory satisfying at best the combination of soft constraints. In green the best zones, in red the zones to avoid.

Acknowledgment The authors would like to thank the French Research Agency (ANR) for funding this study through the ACouStiC project.


Fabian Pedregosa
email: fabian.pedregosa@inria.fr

Elodie Cauvet

Gaël Varoquaux

Christophe Pallier

Bertrand Thirion

Alexandre Gramfort

Learning to rank from medical imaging data

Keywords: fMRI, supervised learning, decoding, ranking

Medical images can be used to predict a clinical score coding for the severity of a disease, a pain level or the complexity of a cognitive task. In all these cases, the predicted variable has a natural order. While a standard classifier discards this information, we would like to take it into account in order to improve prediction performance. A standard linear regression does model such information, however the linearity assumption is likely not be satisfied when predicting from pixel intensities in an image. In this paper we address these modeling challenges with a supervised learning procedure where the model aims to order or rank images. We use a linear model for its robustness in high dimension and its possible interpretation. We show on simulations and two fMRI datasets that this approach is able to predict the correct ordering on pairs of images, yielding higher prediction accuracy than standard regression and multiclass classification techniques.

Introduction

Statistical machine learning has recently gained interest in the field of medical image analysis. It is particularly useful for instance to learn imaging biomarkers for computer aided diagnosis or as a way to obtain additional therapeutic indications. These techniques are particularly relevant in brain imaging [START_REF] Cuingnet | Spatial regularization of SVM for the detection of diffusion alterations associated with stroke outcome[END_REF], where the complexity of the data renders visual inspection or univariate statistics unreliable. A spectacular application of learning from medical images is the prediction of behavior from functional MRI activations [START_REF] Haynes | Decoding mental states from brain activity in humans[END_REF]. Such a supervised learning problem can be based on regression [START_REF] Kay | Identifying natural images from human brain activity[END_REF] or classification tasks [START_REF] Haxby | Distributed and Overlapping Representations of Faces and Objects in Ventral Temporal Cortex[END_REF].

In a classification setting, e.g. using support vector machines (SVM) [START_REF] Laconte | Support vector machines for temporal classification of block design fMRI data[END_REF], class labels are treated as an unordered set. However, it is often the case that labels corresponding to physical quantities can be naturally ordered: clinical scores, pain levels, the intensity of a stimulus or the complexity of a cognitive task are examples of such naturally ordered quantities. Because classification models treat these as a set of classes, the intrinsic order is ignored leading to suboptimal results. On the other hand, in the case of linear regression models such as Lasso [START_REF] Liu | Blockwise coordinate descent procedures for the multi-task lasso, with applications to neural semantic basis discovery[END_REF] or Elastic Net [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF], the explained variable is obtained by a linear combination of the variables, the pixel intensities in the present case, which benefits to the interpretability of the model [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF][START_REF] Michel | Total variation regularization for fMRI-based prediction of behaviour[END_REF]. The limitation of regression models is that they assume linear relationship between the data and the predicted variable. This assumption is a strong limitation in practical cases. For instance, in stroke studies, light disabilities are not well captured by the standard NIHSS score. For this reason most studies use a classification approach, forgoing the quantitative assessment of stroke severity and splitting patients in a small number of different classes. A challenge is therefore to work with a model that is able to learn a non-linear relationship between the data and the target variable.

In this paper, we propose to use a ranking strategy to learn from medical imaging data. Ranking is a type of supervised machine learning problem that has been widely used in web search and information retrieval [START_REF] Burges | From RankNet to LambdaRank to LambdaMART : An overview[END_REF][START_REF] Richardson | Beyond PageRank: machine learning for static ranking[END_REF] and whose goal is to automatically construct an order from the training data. We first detail how the ranking problem can be solved using binary classifiers applied to pairs of images and then provide empirical evidence on simulated data that our approach outperforms standard regression techniques. Finally, we provide results on two fMRI datasets.

Notations We write vectors in bold, a ∈ R n , matrices with capital bold letters, A ∈ R n×n . The dot product between two vectors is denoted a, b . We denote by a = a, a the 2 norm of a vector.

2 Method: Learning to rank with a linear model A dataset consists of n images (resp. volumes) containing p pixels (resp. voxels).

The matrix formed by all images is denoted X ∈ R n×p .

In the supervised learning setting we want to estimate a function f that predicts a target variable from an image, f : R p → Y. For a classification task, Y = {1, 2, 3, ..., k} is a discrete unordered set of labels. The classification error is then given by the number of misclassified images (0-1 loss). On the other hand, for a regression task, Y is a metric space, typically R, and the loss function can take into account the full metric structure, e.g. using the mean squared error.

Here, we consider a problem which shares properties of both cases. As in the classification setting, the class labels form a finite set and as in the regression setting there exists a natural ordering among its elements. One option is to ignore this order and classify each data point into one of the classes. However, this approach ignores valuable structure in the data, which together with the high number of classes and the limited number of images, leads in practice to poor performance. In order to exploit the order of labels, we will use an approach known as ranking or ordinal regression.

Ranking with binary classifiers Suppose now that our output space Y = {r 1 , r 2 , ..., r k } verifies the ordering r 1 ≤ r 2 ≤ .. ≤ r k and that f : R → Y is our prediction function. As in [START_REF] Herbrich | Large margin rank boundaries for ordinal regression[END_REF], we introduce an increasing function θ : R → R and a linear function g(x) = x, w which is related to f by

f (x) = r i ⇐⇒ g(x) ∈ [θ(r i-1 ), θ(r i )[ (1)
Given two images (x i , x j ) and their associated labels (y i , y j ) (y i = y j ) we form a new image x i -x j with label sign(y i -y j ). Because of the linearity of g, predicting the correct ordering of these two images, is equivalent to predicting the sign of

g(x i ) -g(x j ) = x i -x j , w [8].
The learning problem is now cast into a binary classification task that can be solved using standard supervised classification techniques. If the classifier used in this task is a Support Vector Machine Classifier, the model is also known as RankSVM. One of the possible drawbacks of this method is that it requires to consider all possible pairs of images. This scales quadratically with the number of training samples, and the problem soon becomes intractable as the number of samples increases. However, specialized algorithms exist with better asymptotic properties [START_REF] Joachims | Training linear SVMs in linear time[END_REF]. For our study, we used the Support Vector Machine algorithms proposed by scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF].

The main benefit of this approach is that it outputs a linear model even when the function θ is non-linear, and thus ranking approaches are applicable to a wider set of problems than linear regression. Compared to multi-label classification task, where the number of coefficient vectors increase with the number of labels (ranks), the number of coefficients to learn in the pairwise ranking is constant, yielding better-conditioned problems as the number of unique labels increases.

Performance evaluation Using the linear model previously introduced, we denote the estimated coefficients as ŵ ∈ R p . In this case, the prediction function corresponds to the sign of x i -x j , ŵ . This means that the larger x i , ŵ , the more likely the label associated to x i is to be high. Because the function θ is non-decreasing, one can project along the vector ŵ to order a sequence of images. The function θ is generally unknown, so this procedure does not directly give the class labels. However, under special circumstances (as is the case in our empirical study), for example when there is a fixed number of samples per class this can be used to recover the target values.

Since our ranking model operates naturally on pairs of images, we will define an evaluation measure as the mean number of label inversions. Formally, let (x i , y i ) i=1,...,n denote the validation dataset and P = {(i, j) s.t. y i = y j } the set of pairs with different labels. The prediction accuracy is defined as the percentage of incorrect orderings for pairs of images. When working with such a performance metric, the chance level is at 50%:

Error = # {(i, j) ∈ P s.t. (y i -y j )(f (x j ) -f (x i )) < 0} /#P .
(

) 2
Model selection The Ranking SVM model has one regularization parameter denoted C, which we set by nested cross-validation on a grid of 50 geometricallyspaced values between 10 -3 and 10 3 . We use 5-folds splitting of the data: 60% of the data is used for training, 20% for parameter selection and 20% for validation.

To establish significant differences between methods we perform 20 such random data splits.

Results

We present results on simulated fMRI data and two fMRI datasets.

Simulation study

Data generation

The simulated data X contains n = 300 volumes (size 7×7×7), each one consisting of Gaussian white noise smoothed by a Gaussian kernel with standard deviation of 2 voxels. This mimics the spatial correlation structure observed in real fMRI data. The simulated vector of coefficients w has a support restricted to four cubic Regions of Interest (ROIs) of size (2 × 2 × 2). The values of w restricted to these ROIs are {5, 5, -5, -5}.

We define the target value y ∈ R n as a logistic function of Xw:

y = 1 1 + exp (-Xw) + (3)
where ∈ R n is a Gaussian noise with standard deviation γ > 0 chosen such that the signal-to-noise ratio verifies / √ Xw = 10%. Finally, we split the 300 generated images into a training set of 240 images and a validation set of other 60 images.

Results

We compare the ranking framework presented previously with standard approaches. Ridge regression was chosen for its widespread use as a regression technique applied to fMRI data. Due to the non-linear relationship between the data and the target values, we also selected a non-linear regression model: support vector regression (SVR) with a Gaussian kernel [START_REF] Drucker | Support vector regression machines[END_REF]. Finally, we also considered classification models such as multi-class support vector machines. However, due to the large number of classes and the limited number of training samples, these methods were not competitive against its regression counterpart and were not included in the final comparison.

One issue when comparing different models is the qualitatively different variables they estimate: in the regression case it is a continuous variable whereas in the ranking settings it is a discrete set of class labels. To make both comparable, a score function that is applicable to both models must be used. In this case, we used as performance measure the percentage of incorrect orderings for pairs as defined in [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF]. non zero due to noise. All of the model parameters were set by cross-validation on the training set. The figure not only shows that Ranking SVM converges to an optimal model (statistical consistency for prediction), but also that it converges faster, i.e. has lower sample complexity than alternative approaches, thus it should be able to detect statistical effects with less data. Gaussian kernel SVR performs better than ridge regression and also reaches the optimal error. However, the underlying model of SVR is not linear and is therefore not wellsuited for interpretation [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF]; moreover, it is less stable to high-dimensional data.

As stated previously, the function θ can be estimated from the data. In Fig. 1-b we use the knowledge of target values from the validation dataset to estimate this function. Specifically, we display class labels as a function of X ŵ and regularize the result using a local regression (LOWESS). Both estimated function and ground truth overlap for most part of the domain.

Results on two functional MRI datasets

To assess the performance of ranking strategy on real data, we investigate two fMRI datasets. The first dataset, described in [START_REF] Cauvet | Traitement des Structures Syntaxiques dans le langage et dans la musique[END_REF], consists of 34 healthy volunteers scanned while listening to 16 words sentences with five different levels of complexity. These were 1 word constituent phrases (the simplest), 2 words, 4 words, 8 words and 16 words respectively, corresponding to 5 levels of complexity which was used as class label in our experiments. To clarify, a sentence with 16 words using 2 words constituents is formed by a series of 8 pairs of words. Words in each pair have a common meaning but there is meaning between each pair. A sentence has therefore the highest complexity when all the 16 words form a meaningful sentence.

The second dataset is described in [START_REF] Tom | The neural basis of loss aversion in decision-making under risk[END_REF] and further studied in [START_REF] Jimura | Analyses of regional-average activation and multivoxel pattern information tell complementary stories[END_REF] is a gambling task where each of the 17 subjects was asked to accept or reject gambles that offered a 50/50 chance of gaining or losing money. The magnitude of the potential gain and loss was independently varied across 16 levels between trials. No outcomes of these gambles were presented during scanning, but after the scan three gambles were selected at random and played for real money. Each gamble has an amount that can be used as class label. In this experiment, we only considered gain levels, yielding 8 different class labels. This dataset is publicly available from http://openfmri.org as the mixed-gambles task dataset. The features used in both experiments are SPM β-maps, a.k.a. GLM regression coefficients.

Both of these datasets can be investigated with ranking, multi-label classification or regression. We compared the approaches of regression and ranking to test if the added flexibility of the ranking model translates into greater statistical power. Due to the high number of classes and limited number of samples, we found out that multi-label classification did not perform significantly better than chance and thus was not further considered.

To minimize the effects that are non-specific to the task we only consider pairs of images from the same subject.

The first dataset contains four manually labeled regions of interest: Anterior Superior Temporal Sulcus (aSTS), Temporal Pole (TP), Inferior Frontal Gyrus Orbitalis (IFGorb) and Inferior Frontal Gyrus triangularis (IFG tri). We then compare ranking, ridge regression and Gaussian kernel SVR models on each ROI separately. Those results appear in the first three rows of Table 1 and are denoted as language complexity with its corresponding ROI in parenthesis. We observe that the ranking model obtains a significant advantage on 3 of the 4 ROIs. This could be explained by a relatively linear effect in IFGorb or a higher noise level.The last row concerns the second dataset, denoted gambling, where we selected the gain experiment (8 class labels). In this case, since we were not presented manually labeled regions of interest, we performed univariate dimensionality reduction to 500 voxels using ANOVA before fitting the learning models. It can be seen that the prediction accuracy are lower compared to the first dataset. Ranking SVM however still outperforms alternative methods. Locations in the brain of the ROIs for the first dataset, with a color coding for their predictive power are presented in Fig. 2-a. As shown previously in the simulated dataset, Fig. 2-b shows the validation data projected along the coefficients of the linear model and regularized using LOWESS local regression for each one of the four highest ranked regions of interest. Results show that the link function between X ŵ and the target variable y (denoted θ in the methods section) varies in shape across ROIs, suggesting that the BOLD response is not unique over the brain. 1. Prediction accuracy of the ranking strategy on two real fMRI datasets: language complexity (lang. comp.) in 3 ROIs and gambles. As a comparison, scores obtained with alternative regression techniques (ridge regression and SVR using a nonlinear Gaussian kernel) are presented. As confirmed by a Wilcoxon paired test between errors obtained for each fold using Ranking SVM and ridge regression, Ranking SVM leads to significantly better scores than other approaches on 4 of the 5 experiments. We observe that the shape of the curves varies across brain regions.

Discussion and conclusion

In this paper, we describe a ranking strategy that addresses a common use case in the statistical analysis of medical images, which is the prediction of an ordered target variable. Our contribution is to formulate the variable quantification problem as a ranking problem. We present a formulation of the ranking problem that transforms the task into a binary classification problem over pairs of images. This approach makes it possible to use efficient linear classifiers while coping with non-linearities in the data. By doing so we retain the interpretability and favorable behavior in high-dimension of linear models.

From a statistical standpoint, mining medical images is challenging due to the high dimensionality of the data, often thousands of variables, while the number of images available for training is small, typically a few hundreds. In this regard, the benefit of our approach is to retain a linear model with few parameters. It is thus better suited to medical images than multi-class classification.

On simulations we have shown that our problem formulation leads to a better prediction accuracy and lower sample complexity than alternative approaches such as regression or classification techniques. We apply this method to two fMRI datasets and discuss practical considerations when dealing with fMRI data. We confirm the superior prediction accuracy compared to standard regression techniques.

Figure 1 -Fig. 1 .

 11 Fig.1. a) Prediction accuracy as the training size increases. Ranking SVM performs consistently better than the alternative methods, converging faster to the empirical optimal error. b) Estimation of the θ function using non-parametric local regression and ranking SVM. As expected, we recover the logistic function introduced in the data generation section.

  RankSVM Ridge SVR P-val lang. comp (aSTS) 0.706 0.661 0.625 2e-3*** lang. comp. (TP) 0.687 0.645 0.618 7e-4*** lang. comp. (IFGorb) 0.619 0.609 0.539 0.3 lang. comp. (IFG tri) 0.585 0.566 0.533 5e-2* gambling 0.58 0.56 0.53 1e-2** Table

Fig. 2 .

 2 Fig. 2. a) Scores obtained with the Ranking SVM on the 4 different ROIs. The regions with the best predictive power are the temporal pole the anterior superior temporal sulcus. b) The target variable y as a function of X ŵ for the four regions of interest.We observe that the shape of the curves varies across brain regions.