Florent Perez

Joel Cibert

Maria Vladimirova

Denis Scalbert

Spin waves in magnetic quantum wells with Coulomb interaction and sd exchange coupling

Keywords: numbers: 75.30.Ds, 73.21.-b, 85.75.-d, 76.50.+g, 76.30.-v

come

Spin waves in magnetic quantum wells with Coulomb interaction and sd exchange coupling

I. INTRODUCTION

Collective spin dynamics in dilute magnetic semiconductors (DMS) has recently drawn lots of attention. [1][2][3][4][5] This field provides an insight into the origins of carrier-induced ferromagnetism in semiconductors [6][7][8] and particular features in the spin excitation spectrum 8,9 due to the presence of two spin sub-systems that are dynamically coupled by Coulomb-exchange interaction: that of the itinerant carrier and that of the localized magnetic impurities. The transverse spin excitation spectrum has been theoretically found to be composed of three types of excitations. These are: two collective spin waves corresponding to itinerant and localized spins precessing in phase or out of phase to each other, and single-particle (or Stoner-like) excitations of the itinerant carriers. 4,[8][9][10] . If the DMS is in the ferromagnetic state, the in-phase spin wave (IPW) becomes the Goldstone-like mode with an acoustic type dispersion responsible for long-range spin order in the ground state. The out of phase spin wave (OPW) develops an optical branch, with a zone-center energy determined by the strength of Coulomb-exchange interaction between carriers and the spins of magnetic impurities.

Experimental evidence of the entire spectrum in a ferromagnetic DMS like GaMnAs is not available. What has been reported so far are features related to the zone-center IPW, dominated by the Mn spin precession, its dynamics. 2,11,12 and its ferromagnetic resonance. 13 We find no experimental data available for the out of phase mode. Indeed, ferromagnetism in GaMnAs systems requires a high Mn concentration, which destroys the periodicity of the crystal potential and smooths out all optical resonances. More insight into the DMS spin excitation spectrum has been gained in CdMnTe doped quantum wells (QW), which constitute a clean test-bed system, appropriate to capture general properties of the collective spin dynamics in DMS materials. Evidence for carrierinduced ferromagnetism has been found in CdMnTe quantum wells doped with holes 14 .

When doped with electrons, due to the very low Curie temperature, only the paramagnetic phase is available to most experiments. The OPW mode dispersion and single-particle excitations have been probed by Raman measurements in the paramagnetic state 15,16 . The mixed nature of the IPW and OPW waves has been evidenced in the frequency 17 and time domain 1,4 . Neithertheless, there is a lack of a full theoretical description of the spin excitation spectrum in CdMnTe QW. Indeed, so far two approaches were followed to describe the spin excitations : in the first one 10 , the sd-exchange dynamical coupling between Mn and electron spins was considered, but the Coulomb interaction between electrons was dropped out. In the second 18 , the reverse point of view was adopted : spin resolved Coulomb interaction between electrons was taken into account, but only the static mean-field sd contribution from the Mn was kept to form a highly spin-polarized two dimensional electron gas (SP2DEG). This work fills the gap between the two theoretical approaches, by solving the spin dynamical equations in presence of both the sd-exchange dynamical coupling between Mn and electron spins and the Coulomb interaction between electrons. Starting from the full DMS Hamiltonian, the approach combines exact commutation rules and standard generalized Random Phase Approximation (RPA). We also include the intrinsic damping of the pure electron spin waves due to the delocalized character of the electrons 19,20 . We show that the introduction of Coulomb electron-electron interaction induces strong qualitative changes in the spectrum compared to the approach of Ref. 10 and that inclusion of the intrinsic damping diminishes the strength of the coupling between the two spin subsystems for non-zero wavevectors. Generalization of this model to hole systems might be considered : then one should take into account the fact that the hole spin states are not isotropic.

The paper is divided as follows : in Sec. II, we detail the Hamiltonian of the system and rewrite it in terms of collective variables, in Sec. III, we use transverse spin dynamics equations to derive spin response functions, and in the last section, we study the spectrum of spin mixed electron-Mn modes.

II. THE 2DEG DMS HAMILTONIAN UNDER STATIC FIELD

We consider a QW of width w containing x ef f N 0 unpaired 21 Mn spins per unit volume.

The first subband is populated by n 2D electrons per unit surface. The 2DEG-DMS Hamiltonian under the influence of a static magnetic field B = Be z applied in the plane of the QW writes :

Ĥ = ĤKin + ĤCoulomb + Ĥs-d + ĤZeeman (1) Ĥs-d = -α Ŝ (r) • M (r) d 3 r ĤZeeman = g e µ B Ŝ (r) • Bd 3 r + g Mn µ B M (r) • Bd 3 r
where α is the exchange coupling between conduction electrons and Mn spins (α > 0) ,

and g e and g Mn are normal g-factor of, respectively, conduction electrons and Mn electrons.

In the convention where µ B > 0, we have g e -1.44 and g Mn 2.00. We have introduced two vector operators : Ŝ (r) = χ 2 (y) i ŝi δ r // -r i// is the 3D electron spin density in a splitted coordinates frame r = r // , y with r // , the in-plane position and y the out of plane coordinate. χ (y) is the electron envelope-function of the first subband of the QW. The i index accounts for the i-th electron of the 2DEG, its spin 1 2 is described by the operator ŝi and its position is r i// . M (r) = j Îj δ (r -R j ) is the Mn 3D spin density. The j-th 5 2 -spin Îj of a single Mn impurities is localized on the cation site R j . In the equilibrium state at temperature T , each Mn spin has the average value Îz (B, T ) , which is the thermodynamic average over the five occupied states of the Mn atom d-shell, given by the modified Brillouin function 21 . The 2DEG has the equilibrium spin polarization ζ = (n ↑ -n ↓ ) /n 2D .

We, now, rewrite the s -d Hamiltonian using the electron (and Mn) spin fluctuations operators at in-plane wave-vector q, respectively, Ŝq = Ŝ (r) e -iq•r // d 2 r // dy = i ŝi e -iq•r i//

and Mq = j Îj e -iq•R j// . Due to the 2D and 3D characters of, respectively, the conduction electron and Mn spins subsystems, the electron spin-degrees of freedom naturally couples to Mn spin profile weighted by the squared electron wave-function. For later convenience, we introduce the following n-profile Mn spin fluctuations operators :

M(n) q = w n χ 2n (y) M (r) e -iq•r // d 2 r // dy = w n j χ 2n (y j ) Îj e -iq•R j//
Hence, Mq = M(0) q . M(n) q are vector operators verifying the following commuting relations :

M (n) α,q , M (p) β,q = i α,β,γ M (n+p) γ,q+q (2)
where α, β, γ = x, y, z and α,β,γ is the Levi-Cevita tensor. It follows :

Ĥs-d = - α wL 2 q Ŝq • M(1) -q (3a) = -α          δ Ŝz,q=0 • M (1) z,q=0 0 + Ŝz,q=0 0 • δ M (1) z,q=0 + δ Ŝz,q=0 • δ M (1)
z,q=0

+ 1 2 Ŝ+,q=0 • M (1) -,q=0 + 1 2 Ŝ-,q=0 • M (1)
+,q=0

+ q =0 Ŝq • M(1) -q          (3b)
In Eq. (3a), we have defined the exchange coupling constant α = α/wL 2 , equilibrium averaging 0 and fluctuation operators :

δ Â = Â -Â 0 .
Finally, we get :

Ĥs-d + H Zeeman = Z e (B) • δ Ŝz,q=0 + g Mn µ B B • δ M (0) z,q=0 + K • δ M (1) z,q=0 (4a)
- α 1 2 Ŝ+,q=0 • M (1) -,q=0 -Ŝ-,q=0 • M (1) +,q=0 (4b)
-αδ Ŝz,q=0 • δ M (1) z,q=0 (4c)
-

α q =0 Ŝq • M (1) -q (4d)
with the total bare Zeeman energy of conduction electrons given by :

Z e (B) = αγ 1 N Mn Îz (B, T ) ∆ -|g e | µ B B (5)
where γ 1 = w 0 χ 2 (y) dy is the probability to find the electron in the QW, N Mn = x ef f N 0 wL 2 is the number of Mn spins available in the QW. In Eq. (5), we evidence the "Overhauser shift" ∆ and the normal Zeeman contribution, opposite to ∆ (the sd contribution). Indeed, as g e < 0, g Mn > 0 and α > 0, the Mn spins are anti-parallel to the field, thus, the sd coupling aligns the electron spins anti-parallel to the field, while the normal Zeeman aligns the electron spin parallel to the field. When the sd coupling dominates over the normal Zeeman contribution and B > 0, both Mn and electron spins are ↓ : Îz < 0 and the 2DEG spin polarization degree ζ is also negative.

The first line, Eq. (4a) gives the mean-field "effective Zeeman" Hamiltonian, where nodynamical coupling between electrons and Mn spins appears. The mean-field Hamiltonian naturally introduces the "Knight shift", due to the equilibrium electron spin polarization S z,q=0 0 = n 2D L 2 ζ/2, which shifts the Mn spin precession energy:

K = -α S z,q=0 0 = 1 2 α w n 2D |ζ| (6)
Eq. (4b) and Eq. (4d) give rise to a first order sd-dynamical coupling between transverse spin degrees of freedom and induce spin-mixed electron-Mn modes of precession. They also contain higher orders correlation terms which have been indentified to be responsible for a damping 5 . The effects of the correlations contained in Eqs. (4b)-(4d) are out of the scope of the present work and will be neglected when they appear.

III. TRANSVERSE SPIN DYNAMICS

A. SP2DEG dynamics without the s-d dynamical coupling

In this paragraph, we take into account only the first line of Eqs.(4a)-(4d). Hence conduction electron and Mn dynamics are independent. The 2DEG is polarized by the static exchange field of Mn spins, this forms a spin-polarized 2DEG (SP2DEG) as described in Ref. 18 . The Hamiltonian which rules the electron dynamics in the SP2DEG reduces to :

ĤSP2DEG = ĤKin + ĤCoulomb + Z e (B) Ŝz,q=0 (7)
Introducing the the electron creation-anihilation operators, a spin-flip single particle excitation (SF-SPE) is described by a single electron-hole pair operator c + k-q,↑ c k,↓ , and electrons spin-wave operators, introduced above, are given by Ŝ+,q = Ŝx,q + i Ŝy,q = k c + k-q,↑ c k,↓ . Let's notice that Ŝ+,q , c + k-q ,↑ c k,↓ = 0, such that collective and single particle modes are not intrinsically coupled. In the following we will use exact commutation rules to write equation of motions for these normal modes of the SP2DEG. We will note d Â/dt 2DEG = Â, ĤSP2DEG /i the time derivative of  related to ĤSP2DEG only. Further, we will make use of linear response theory to derive quantities like :

Â; B ω = - i lim ε→0 + +∞ 0 Â (t) , B 0 e -iωt-εt dt (8)
which gives the linear response of an observable  to a perturbation coupled linearly to B in the considered Hamiltonian 22 . Â; B ω has the following equations of motion :

Â; B ω = - i ω Ȧ; B ω - 1 ω Â, B 0 (9a) = i ω Â; Ḃ ω - 1 ω Â, B 0 (9b)
where Ȧ = Â, Ĥ /i is the time derivative of  related to the considered Hamiltonian.

We will note Â; B

2DEG ω when Ȧ is replaced by d Â/dt 2DEG .

Single particle modes dynamics

The kinetic Hamiltonian ĤKin = k,σ E k c + k,σ c k,σ and the mean-field Zeeman Hamiltonian conserve the single-particle modes :

c + k-q,↑ c k,↓ , ĤKin + Z e (B) Ŝz,q=0 = (E k -E k-q -Z e ) c + k-q,↑ c k,↓ But the Coulomb Hamiltonian ĤCoulomb = 1 2 q =0,k ,σ ,k,σ V q c + k+q,σ c + k -q,σ c k ,σ c k,σ
couples a single particule mode to multi-pair modes having a spin +1 :

c + k-q,↑ c k,↓ , ĤCoulomb = (10)
+ q =0,k ,σ V q c + k-q+q ,↑ c k ,σ c + k -q ,σ c k,↓ - q =0,k ,σ V q c + k-q,↑ c k ,σ c + k -q ,σ c k-q ,↓
where V q is the space Fourier transform of the bare Coulomb interaction. It follows from Eq. ( 10) that ĤCoulomb does not conserve the SPE motion, but introduces an infinite hierarchy where a single electron-hole pair of the Fermi sea (a SPE) couples to multiple pairs having the same global spin. Approximations can be made : in the rhs terms of Eq.

(10), some conserve the SPE motion and renormalize it, others introduce a scattering effect, the so-called spin-Coulomb drag 23 , and can be described by an electron-electron scattering time 24 τ e-e . The former consists in making the random phase approximation (RPA) on single mode dynamics 25 , i.e., keeping in Eq. ( 10) only terms which can be written as a product of a SF-SPE with an occupation number nk,σ , and replacing nk,σ by its average value nk,σ 0 .

Then the Coulomb factor V q has to be replaced by a local field factor G xc which accounts for the effective dynamical exchange-field produced by other electrons 26 (a part of what has been suppressed in making the RPA). Adding a damping rate η = /τ e-e due to the scattering leads to the SF-SPE equation of motion that we will use in the following :

i d dt c + k-q,↑ c k,↓ 2DEG = E k -E k-q -Z e + 2G xc Ŝz,q=0 0 c + k-q,↑ c k,↓ (11)
+G xc nk,↓ 0 -nk-q,↑ 0 Ŝ+,q -iηc + k-q,↑ c k,↓
Eq. ( 11) evidences the renormalized Zeeman energy, i.e., the spin-flip energy of single electrons :

Z * = Z e -2G xc Ŝz,q=0 0 (12)
Compared to the bare Zeeman energy Z e , Z * is enhanced by Coulomb-exchange between spin-polarized electrons, a phenomenon linked to the spin-susceptibility enhancement 27 .

Each SF-SPE is characterized by two wavevectors k and q. At q = 0, SF-SPE are degenerate to Z * . When q = 0 the degeneracy is lifted by the kinetic spread of velocities which depend on the initial momentum k.

Collective modes dynamics

Along the SF-SPE, the above spin polarized SP2DEG develops collective modes, the socalled spin-flip waves (SFW). SFW dynamics are described by the Ŝ+,q operators dynamics.

As ĤCoulomb conserves the macroscopic spin, it follows :

Ŝ+,q , ĤCoulomb + Z e (B) Ŝz,q=0 = -Z e (B) Ŝ+,q (13)
But, the kinetic Hamiltonian couples collective states to the spin current Ĵ+,q = m * k k-q 2 c + k-q,↑ c k,↓ carried by single particle states :

Ŝ+,q , ĤKin = q • Ĵ+,q (14)
Finally the equation of motion of collective SFW modes writes :

d dt Ŝ+,q 2DEG = iω e Ŝ+,q -iq • Ĵ+,q (15)
where ω e = Z e (B) / is the frequency of the Larmor's electron mode.

One is left with evaluating the spin-current dynamics to find the 2DEG electron spin waves. The spin current evolution is dominated by single particle states dynamics as ĤCoulomb does not conserve Ĵ+,q and destroys the coherence between the single particle objects composing Ĵ+,q . As seen from Eq. ( 11), the exchange field produced by the spin fluctuation Ŝ+,q drives the spin current. The interplay between the spin-current dynamics and the Ŝ+,q dynamics then determines both the SFW dispersion and its damping. The relevant spincurrent response is the transverse spin-conductivity 20 σ+ , which links the spin current to the gradient of the exciting exchange field :

Ĵ+,q ω = qσ + (q, ω) G xc Ŝ+,q ω ( 16
)
where ω is the expectation value at frequency ω. The spin-conductivity has an imaginary part originating from the damping of SF-SPE, intrinsically due to ĤCoulomb or any source of disorder acting on transverse spin degrees of freedom. Consequently, the real part of the spin conductivity determines the SFW dispersion, while the imaginary part determines its damping. It is worth noting, that the spin wave damping originates from the kinetic motion of the conduction electrons and from the topoly of the conduction band, a 2D parabolla.

In a Luttinger liquid 28 , the conduction band is linear in k and 1D, thus Eq. ( 14) conserves the macroscopic spin. It breaks the coupling between Ŝ+,q and SF-SPEs, which are coupled to charge degrees of freedom by ĤCoulomb . This property is at the origin of the well known spin-charge separation 29 occuring in Luttinger liquids. Injecting Eq. ( 16) into Eq. ( 15) and solving it in the frequency domain for long wavelength (q k F ) leads to :

d dt Ŝ+,q = iω q Ŝ+,q (17)
with ωq a complex pulsation:

Re ωq = ω e -q 2 G xc lim q→0,ω→0

Re σ+ (18a)

Im ωq = q 2 G xc lim q→0,ω→0
Im σ+ (18b)

In the following we note σ + = lim q→0,ω→0 Im σ+ , the imaginary spin-conductivity. It was calculated in Ref. 20 and some corrections were added in Ref. 19 which gave also an experimental evidence of the kinetic damping law found in Eq. (18b). We highlight that these q 2 laws are valid in the longwavelength limit when the SFW propagates far from the SF-SPE continuum (see Ref. 18 ). When close to this continuum, the stronger coupling with SF-SPEs introduces corrections to the above laws and one should better replace the dispersion law with the pole appearing in the transverse spin susceptibility (see Ref. 18 ) which will be derived in the next paragraph.

Electron spin-susceptibility

The transverse spin-susceptibility, defined by the ratio of the expectation value Ŝ+,q ω to the perturbing potentiel g e µ B b +,qω , where b +,qω is the amplitude at the same pulsation ω and wavevector q of the exciting magnetic field, is given by :

χ + (q, ω) = Ŝ+,q ; Ŝ-,-q 2DEG ω (19)
Straightforward calculations using the equation of motion Eq. ( 11) lead to :

χ + (q, ω) = Ŝ+,q ; Ŝ-,-q 2DEG ω = Π ↓↑ (q, ω) 1 + G xc Π ↓↑ (q, ω) (20)
where we have introduced the transverse Lindhardt-type response 18 :

Π ↓↑ (q, ω) = k nk-q,↑ 0 -nk,↓ 0 Z * + E k-q -E k -ω -η (21)
A comparison between the above spin-susceptibility expression and the one given by local spin-density approximation 18 , gives the expression of the local field factor G xc :

G xc = - 2 n 2 2D L 2 1 ζ ∂E xc ∂ζ (22)
where E xc is the exchange-correlation part of the ground state energy 30 .

SFW appear as poles of χ + (q, ω) , one finds in the long wavelength limit, another expression for Re ωq :

Re ωq = ω e - 1 |ζ| Z e Z * -Z e 2m * q 2 (23)
Alternatively, if one uses the approximated equation of motion Eq. ( 15), one finds the spin susceptibility in the long wavelength limit :

χ + (q, ω) = - 2 Ŝz,q=0 0 ω -ωq (24)

B. Transverse spin dynamics equations with s-d dynamical coupling

Now, we keep lines (4b) to (4d) in the s-d Hamiltonian, and we reconsider collective transverse spin dynamics. In the following, the derivative d Â/dt = Â, Ĥ /i takes into account the coherent coupled dynamics due to lines (4b) and (4d) in the s-d Hamiltonian, but reduced to first order terms : higher order correlation terms like q δ Ŝz,q+q • M

+,-q have been dropped.

Electron dynamics

We find :

d dt Ŝ+,q = Ṡ+,q - i K M (1) +,q (25)
where Ṡ+,q = iω q Ŝ+,q . Compared to the SP2DEG dynamics, the sd-dynamical coupling adds the second term of Eq. ( 25) which is a coherent coupling with Mn transverse degrees of freedom. One key feature is that the collective electron motion naturally couples with M (1) +,q Mn-modes, a Mn precession having a profile, out of the QW plane, following the electron probability distribution. We are left with deriving the equation of motion for these Mn-modes.

Manganese dynamics

We obtain the first order equation of motion for Mn spins :

d dt M (n) +,q = i g M n µ B B M (n) +,q + i K M (n+1) +,q - i ∆ n+1 Ŝ+,q (26)
where we have introduced n-profile Overhauser shifts :

∆ n = α M (n) z,q=0 0 = γ n /γ 1 ∆
with γ n = w n-1 w 0 χ 2n (y) dy. The important features are the second and third terms in Eq. ( 26). The later couples the Mn-precession with collective electron modes. The former couples a n-profile Mn mode to a (n + 1)-profile mode, because this coupling is mediated by the 2DEG. Thus, the Mn-dynamics is given by an infinite serie of equations. This is a consequence of the 3D nature of the Mn dynamics. A variable like M (n) +,q describes an oscillation propagating in the plane with a rigid profile in the normal direction, but the out of plane degree of freedom is restored by the possibility for Mn spins to build modes which are combinations of M (n) +,q resulting in different out of plane profiles 10 . Obviously, the M (n) +,q

are not independent variables because they don't correspond to orthogonal out of plane profiles. Solving the serie of infinite equations requires a projection of M (n) +,q over a set of modes with orthogonal profiles as it was carried out in Ref. 10 . Along with modes having a strong mixed nature (electron-Mn modes), we then expect to find a high number of modes having essentially a Mn character, but with orthogonal profiles (Mn modes). The number of Mn modes has to be consistent with the initial number of degrees of freedom present in the system. Ref. 10 found a high number of Mn modes branches which were separated by energies of the order of 0.1∆. However, in the experimental data of Ref. 4 , only one branch of these Mn modes was apparent. It appears then, that the set of modes chosen in Ref. 10 is not the most appropriate to describe properly all the modes contained in Eqs.( 25)-( 26), at least in the vicinity of the anticrossing gap (see below). Anyway, this point requires further developments out of the scope of the present study. Indeed, we are particularly interested in discussing mixed electron-Mn modes which are strongly coupled to electrons rather than modes specific to the 3D nature of the Mn dynamics. We can remark that the coupling between M (n) +,q and M (n+1)

+,q has a strength given by K, which is very small compared to

∆ n due to the ratio M (n) z,q=0 0 / Ŝz,q=0 0
1. Hence, considering only modes strongly coupled with electron modes is reasonable. As electron modes are naturally coupled to M

+,q modes, we will consider the dynamics for these ones only by cutting the infinite serie with Ŝ+,q ; Ŝ-,-q

ω = ( ω -ω M n ) χ + (q, ω) ω -ω M n -α∆ 2 2 χ + (q, ω) (30)
and

M (1) +,q ; M (1) -,-q ω = -2 M (1) z,q=0 ω -ω M n -α∆ 2 2 χ + (q, ω) (31)
Consequently, e-Mn mixed spin excitations appear as poles of the above responses, i.e., are zeros of the propagator :

ω -ω M n - α∆ 2 2 χ + (q, ω) (32)
with χ + (q,ω) being the spin-susceptibility of the uncoupled SP2DEG described in Section III B.

We can understand the above equation as follows. Consider the Mn point of view ;

in the presence of the SP2DEG, the precession frequency of Mn spins is shifted from the normal precession (g M n µ B B/ ) by two quantities : a blue shift due to static exchange field with spin polarized electrons (K/ ) and an additional shift due to the dynamic change of the electron spin-polarization. The later is induced by the Mn precession itself. Finally Eq. ( 32) describes a recursive closed loop where : Mn transverse precession induces electron transverse precession proportional to ∆ 2 χ + (q,ω) , this dynamically changes the electron spin polarization which in turn shifts the Mn precession frequency by an amount α∆ 2 χ + (q,ω) .

Finally, in dropping correlation terms given by Eqs. (4c) and (4d), one finds a collective behavior where electrons and Mn respond adiabatically to the dynamical perturbation from the opposite spin-subsystem.

Similar expressions for the coupled modes propagator have been obtained in previous works. To our knowledge it was first derived in Ref. 9 for bulk DMS, and more recently using a spin-path integral approach in DMS quantum wells with electrons 6,10 or bulk DMS with holes 8 . However, none of these works did include the influence of the Coulomb interaction between carriers. Instead of Eq. (32), they resulted in the following propagator :

ω -ω M n - α∆ 2 Π ↓↑ (q,ω) (33)
where the electron spin-susceptibility was replaced by the non-interacting single-particle response introduced in Eq. ( 21). Introduction of Coulomb interaction results in strong qualitative changes in the spectrum which have been partially adressed in Ref. 15 and will be detailed below. Eq. (33) was used to successfully fit the experiment of Teran et al. (Ref. 17 ) where homogenous modes (q=0) were probed and shown to experience an anticrossing at a magnetic field B R such that ω M n = ω e . It is a consequence of the Larmor's theorem that the homogeneous electron mode behaves as if electron were not interacting (Eq. ( 15) for q = 0). Indeed, setting q = 0 in Eq. (32), leads to the homogenous precession modes equation :

B. Homogeneous modes

( ω -ω M n ) ( ω -ω e ) -K∆ 2 = 0 (34)
which solutions are real :

ω ± q=0 = ω e + ω M n 2 ± 1 2 (ω e -ω M n ) 2 + 4K∆ 2 / 2 (35)
Figure 1 shows the magnetic field dependence of these modes and the gap opening at the resonant field B R . The upper branch ω + q=0 (resp ω - q=0 ) has an electron character (resp. Mn character) when B < B R and vice-versa for B > B R . The amplitude of the homogeneous anticrossing gap ∆ω q=0 = √ 4K∆ 2 denotes the strength of the dynamical coupling between the two spin subsystems. Detailed discussions on this gap have been given in Ref. 4,6 and Ref. 1 . In particular, Ref. 1 identifies the anticrossing Mn mode as M (1) +,q consistent with the homothetic approximation of Eq. ( 27) used here. The anticrossing gap was found to be

∆ω q=0 = 4K∆ 2 -T 2e 2
where T 2e is the damping rate of the homogeneous uncoupled electron mode, a quantity that we have neglected here in dropping electron-Mn correlation terms contained in Eqs.(4b)-(4d). In Ref. 1 , T 2e was estimated from measurements of the electronic spin wave damping at q = 0. A rigorous simultaneous determination of ∆ 2 and T 2e , lead to the extraction of K from the anticrossing gap and furthermore to the spin-polarization degree ζ of the 2DEG. Data showed that the so-extracted ζ was slightly exceeding the prediction 18 made for ζ in contradiction with other determinations of ζ performed in the same type of samples 31 , which showed that the model used to predict ζwas reliable. A more accurate description of the anticrossing gap taking into account the infinite set of coupled n-profile mode M (n) +,q equation of motion might overcome this discrepancy. One should also mention, that the infinite serie of equations must be cut in order to conserve the initial number of degrees of freedom (number of available spins in the system). But finding the right number of Eq. ( 26)-like to be taken depends on how the total number of degrees of freedom separates into a number of (quasi-) individual modes 4 and a number of collective electron-Mn modes. Determining this separation is also an important and interesting issue.

C. Spin waves

For q > 0, Eqs(32) and (33) give very different qualitative results as illustrated on Fig. 2.

Without Coulomb interaction between electrons, uncoupled modes of the electrons are the SF-SPE which are degenerate to Z e at q = 0. The sd dynamical coupling introduces two additional collective mode : the OPW propagating above the SF-SPE domain with a positive dispersion and the IPW propagating below with a negative dispersion. Introducing the In ususal conditions (see Ref. 15 ) both the OPW and IPW propagate below the SF-SPE continuum with negative dispersions. Dipersions were calculated after setting to zero the kinetic damping rate.

Coulomb interaction between electrons shifts the SF-SPE to higher energies (Z * ) and givs rise to the collective wave SFW propagating below the SF-SPE continuum. The SFW is further coupled to Mn modes through the sd interaction. An evaluation of the coupling between SF-SPE and Mn modes was given in Ref. 1 and found to be negligible. Thus the Coulomb interaction introduces a shift between the SF-SPE and the SFW energies, and the later is further shifted by the sd dynamical coupling. In realistic conditions, it was shown in Ref. 15 that the Coulomb shifts dominates over the sd dynamical shift. Hence, when Coulomb interaction is taken into account, Eq. (32), except under unrealistic conditions, gives rise to two spin wave modes propagating below the SF-SPE continuum, the IPW and OPW. An anticrossing gap opens at a specific wavevector q R (B) given by :

q R (B) = |ζ| (Z * /Z e -1) 2m * (ω e -ω M n ) (36) Note that q R (B R ) = 0 and that ωq R = ω M n + iq 2 R G xc σ + . If q R (
B) > 0, compared to the homogenous gap, the anticrossing gap at q R is dramatically reduced by the kinetic damping of the electron wave and is given by :

ω + q R -ω - q R = 4K∆ 2 / 2 -η 2 q R (37)
where η q R = q 2 R G xc σ + . We note that the kinetic damping is the only one considered here. Other sources of damping, as e.g., the ones dropped in Eq. (4d), will of course further reduce the amplitude of the gap. Fig. 3 illustrates the variation of q R (B) with the magnetic field. It is always smaller than q m = k F ↓ -k F ↑ , the wavevector delimiting the window where the SFW propagates 18 .

Overlaid in Fig. 3, are the anticrossing gap at q R and the corresponding damping rate η q R . In the absence of the kinetic damping, the gap would be given by √ 4K∆ 2 . One sees the dramatic effect of this intrinsic kinetic damping, which kills the gap outside a very narrow range of magnetic fields. As coupling between spin waves of the electron and the Mn spin systems is responsible for the appearance of the carrier induced ferromagnetism 6 , we might conclude that the above disappearance of the gap diminishes the possibilities for ferromagnetic transitions with complex order (out of q = 0).

The disappearance of the gap is illustrated in Fig. 4 by comparing the dispersions obtained from the zeros of the propagator in Eq. (32) in the presence or absence of the kinetic damping. In the presence of the damping, the solutions have a non-zero imaginary part for q > 0. The corresponding damping rate is plotted in the lower insets of Fig. 4. It is well known that when the frequencies of two coupled oscillators anticross each other, their corresponding damping rates cross themselves. Clearly for B < 5.7 T, the mixed modes do instead of τ e-e (∼ 150ps) to match the experimental conditions of Ref. 19 .

not anticross at any q and each branch conserves its former character, Mn-like or electronlike. On the contrary, for 5.7 T< B < B R , the modes anticross at q R , and the OPW transfers the kinetic damping (q 2 law) of the SFW to the IPW when q > q R . It is worth to note that this q 2 law for the IPW damping rate was also found in GaMnAs compounds in the ferromagnetic state 5 .

V. CONCLUSION

In conclusion, we have introduced equations for the spin dynamics in a test-bed diluted magnetic system that allow to take into account the interplay between the Coulomb interaction dynamics and the sd dynamical coupling between the electrons and the localized spins.

We have shown how the Coulomb interaction introduces strong qualitative changes : the mixed electron-Mn modes propagate below the SF-SPE continuum and an anticrossing gap is open for a given range of magnetic field. Because of Coulomb interaction, the intrinsic kinetic damping due to the electron motion is always present (the SF-SPE scattering time can not be longer than τ e-e ), this damping kills the anticrossing gap outside a very narrow range

FIG. 1 :

 1 FIG. 1: (Color online). Zone center electron-Mn modes. Dotted lines are the uncoupled electron (upper curve) and Mn (lower curve) modes. Full lines are the solutions ω ± q=0 of Eq.(34). Out of the resonant field B R where the modes anticross, branches have electron or Mn character. Sample parameters are, x ef f = 0.23%, T = 2 K, w = 150 Å and n 2D = 3.1 × 10 11 cm -2 .

FIG. 2 :

 2 FIG. 2: (Color online). Illustration of the changes introduced by the Coulomb interaction on the mixed modes dispersions for the same sample parameters as in Fig.1. (a) Solutions without Coulomb given by Eq.(33) for B = 5.8 T. At q = 0 SF-SPE are degenerate to Z e . The sd dynamical shift introduces two propagating waves : above (resp. below) the SF-SPE domain, the OPW (resp.IPW) propagates with a positive (resp. negative) dispersion. The dashed line is the uncoupled Mn mode (degenerate to ω M n ). (b) and (c) The Coulomb interaction between electrons is included as in Eq.(32) of this work. Dashed lines are the uncoupled electron and Mn modes. At q = 0 the SF-SPE energy is shifted to Z * . An anticrossing gap opens at the wavevector q R (B) if B ≤ B R .

FIG. 3 :

 3 FIG. 3: (Color online). Left axis : variation of the wavevector q R (B) and the wavevector q m = k F ↓ -k F ↑ . Right axis : variation of the anticrossing gap and the kinetic damping rate η q R = Im ωq R .The anticrossing gap is killed by the damping rate when the later is of the same order of magnitude as √ 4K∆ 2 . A break in the horizontal axis scale has been introduced to zoom the region close to B B R . To calculate the damping rate, we have used a typical SF-SPE scattering time τ = 2ps

FIG. 4 :

 4 FIG. 4: (Color online). Illustration of the progressive disappearance of the anticrossing gap when B goes away from B R . Continuous lines (open symbols) are dispersions calculated in absence (presence) of the kinetic damping. Dashed lines are the uncoupled modes. Upper insets : zoom on the anticrossing gap region. Lower insets : variation with q of the damping rate of both the OPW (dashed line) and IPW(straight line).

Acknowledgments

The authors would like to thank I. d'Amico, E. Hankiewicz and the GOSPININFO consortium for fruitfull discussions as well as the grant ANR GOSPININFO for financial support.

an homothetic approximation :

+,q = (γ 2 /γ 1 ) M (1)

Consequently the set of coupled electron-Mn equations reduces to :

d dt

where :

is the natural precession pulsation of the free M (1) +,q mode.

IV. MIXED MN-ELECTRON SPIN WAVES A. Spin susceptibilities

To find the dynamically coupled modes, we will derive the electron spin susceptibility with help of equations of motion (9a) and (28a)-(28b).

From Eqs.(9a) and (28a), we first get :

Ŝz,q=0 0 then from Eq. (28b), we get :

+,q ; Ŝ-,-q

which finally leads to : of magnetic fields. Our calculations illustrate also how this kinetic damping is transfered to the IPW, a phenomenon found in GaMnAs compounds.


Mouloud Ourak

Brahim Tamadazte

Nicolas Andreff Femto-St

Registration using wavelet coefficients in spectral domain

This paper focuses on the development of an efficient registration method of Optical Coherence Tomography (OCT) images. The proposed technique is based on the use of a wavelet representation of the OCT images in the registration process. The validation tests show that beside non-stationary and highly inhomogeneous OCT images, the wavelet-based registration demonstrates very interesting results.

Introduction

A biopsy is a medical intervention usually performed by a clinician who takes, from the body, a sample of a suspect tissue for an ex-situ examination [START_REF] Zerbino | Biopsy: its history, current and future outlook[END_REF]. The latter is performed using usually a microscope to detect signs of damage or disease. Biopsies can be performed on all parts of the human body even in the brain. Thus, there are various types of biopsy processes: needle biopsy, endoscopic biopsy, surgical biopsy, etc. Most of these methods are invasive procedures even minimally and generally follow four essential steps: sampling, transfer, histological analysis, and then diagnostic. This takes time, requires a short hospitalization of the patient, a general or local anesthesia, and especially can cause infections in the tissue collection [START_REF] Greco | Breast cancer patients treated without axillary surgery: clinical implications and biologic analysis[END_REF].

Furthermore, a biopsy may be performed in another way using specific medical imaging systems. This is usually known as optical biopsy and consists of an in-situ and real-time optical analysis of the suspect tissue without any removal. To perform an optical biopsy, it requires miniature and high resolution imaging systems as an optical endomicroscopic device. Doubtless, the two most popular techniques used in optical biopsy are confocal microscopy [START_REF] Pawley | Handbook of biological confocal microscopy[END_REF] and OCT [START_REF] Huang | Optical coherence tomography[END_REF]. Over the past two decades, OCT has emerged as the most suitable technique for carrying optical biopsies. This is mainly due to a better compromise between resolution and depth of penetration of tissues compared to the confocal microscopy technique.

This paper deals with the use of OCT images as a means for an in-situ chronological monitoring diagnosis of suspect tissues by tacking several optical biopsies of the same tissues at different time intervals. To be able to compare the different optical biopsies and detect the evolution of the suspect tissue, it is necessary to represent all the biopsies in the same reference. This amounts to perform a registration between successive OCT images [START_REF] Zitová | Image registration methods: a survey[END_REF]. To ensure an accurate registration of the high resolution OCT images, we provide a novel method based on wavelet representation in the spectral domain of the OCT signals.

Wavelet Transform

The wavelet transform allows to represent a 1D or 2D indeed 3D signal by orthonormal series [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The wavelet is widely used for the study of non-stationary and highly inhomogeneous signals in both time and frequency domains. It decomposes the signal into a set of elementary waveforms called wavelets defined by the wavelet coefficients [START_REF] Zerbino | Biopsy: its history, current and future outlook[END_REF].

WI(b, a, θ) +∞ -∞ I(x)ψ r -θ (x -b) a d 2 x, (1)
where WI are the wavelet coefficients, ψ is the mother wavelet, I(x) ∈ R 2 is an image, a ∈ R + * is the scale parameter of the wavelet, b ∈ R 2 translation parameters along x and y, respectively, and r is the 3×3 rotation matrix carried by the angle θ. Thereby, r is given by

r =   cos(θ) -sin(θ) 0 sin(θ) cos(θ) 0 0 0 1   (2)
In order to maintain the continuity of the signal I, we choose ψ as the mother complex oriented Morlet wavelet [START_REF] Weiss | Wavelets and wideband correlation processing[END_REF] defined in the spatial domain by considering the correction exponential term e -1 2 k 2 ψ as follows:

ψ(x, y) = e ik ψ x -e -1 2 k 2 ψ e -1 2 (x 2 +y 2 ) (3)
where k ψ is the wavenumber associate to the Morlet wavelet.

Wavelet-based Registration

The registration defined as the process to align current image I c with reference image I r . This returns to find the minimum of the cost-function defined by:

T = arg min[WI r , WI c ] (4)
where T ∈ (SE(2)×SO( 1)) is the estimated rigid transformation and WI r and WI c are the wavelet coefficients of I r and I c respectively. Therefore, using a log-polar phase correlation method [START_REF] Reddy | An fft-based technique for translation, rotation, and scale-invariant image registration[END_REF] it is possible to represent the cost-function between two images to be registered by the convex shape curve shown in Figure 1. Note the minimum of this cost-function is well defined and can be easily computed by a simple gradient-based optimization algorithm.

Experimental Validation

Having not yet experimental set-up for a physical registration, we opted for a 4 dof (x, y, z, and θ) numerical registration. To ensure this, we consider a first OCT image I r (reference image shown in Figure 3(A)) taken under nominal conditions and a second image I c (image to register shown in Figure 3(B)) taken using another point of view and under extreme conditions of saturation. Figure 3(C) and (D) represent the corresponding wavelet coefficients of I r and I c , respectively. The registration result is illustrated in Figure 3.(E) which is the combination "⊕" of both images I r and I c after transformation. By analyzing the image result, it can be highlighted the high accuracy of the proposed registration method despite the high difference between I r and I c . The registration process takes only 130 milliseconds for 256 × 256 pixels images, this using a 3.5 GHz PC.

Conclusion

In this paper, it was shown using a wavelet-based registration method in the case of OCT images. Despite of the large difference between the OCT images to be registered (i.e., nominal acquisition conditions vs. saturation acquisition conditions), the proposed method was illustrated that the wavelet-based registration algorithm remains efficient with high level of accuracy of convergence.

The next stages of this work will involve to implementation of the developed method in an experimental set-up (including a multiple degrees of freedom) in order to achieve a physical registration during repetitive optical biopsies.

Figure 1 :Figure 2 :

 12 Figure 1: Illustration of the cost-function used in the developed registration method.

Figure 3 :

 3 Figure 3: Registration results: (A) OCT reference image, (B) OCT image to be registered, (C) and (D) the corresponding wavelet coefficients, respectively, and (E) the result of the registration process.

Acknowledgment

This work was supported by the French ANR Labex-ACTION no ANR-11-LABX-01-01.


Matthieu Lê

Hervé Delingette

Jayashree Kalpathy-Cramer

Elizabeth R Gerstner

Tracy Batchelor

Jan Unkelbach

Nicholas Ayache

MRI Based Bayesian Personalization of a Tumor Growth Model

Keywords: tumor growth, glioblastoma, reaction-diffusion, personalization, bayesian, LBM, Monte Carlo

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION A. Motivations

T HE objectives of the mathematical modeling of brain tumor growth are fourfold. First, it could help to better understand the mechanisms behind the disease progression [START_REF] Raman | Computational trials: Unraveling motility phenotypes, progression patterns, and treatment options for glioblastoma multiforme[END_REF]. Second, the personalization of such models to specific patients could allow one to quantify the aggressiveness of the tumor, which has been shown to be correlated with clinically relevant information [START_REF] Neal | Response classification based on a minimal model of glioblastoma growth is prognostic for clinical outcomes and distinguishes progression from pseudoprogression[END_REF], [START_REF] Wang | Prognostic significance of growth kinetics in newly diagnosed GBM revealed by combining serial imaging with a novel biomathematical model[END_REF]. Third, personalized models could help predicting the future evolution of a given tumor [START_REF] Stretton | Predicting the location of glioma recurrence after a resection surgery[END_REF]. Fourth, personalized models could lead the way toward objective and more efficient personalized therapy. For instance, it has already been used in order to personalized radiotherapy plans [START_REF] Corwin | Toward patient-specific, biologically optimized radiation therapy plans for the treatment of glioblastoma[END_REF], [START_REF] Rockne | A mathematical model for brain tumor response to radiation therapy[END_REF]. In this article, we present a method to study the uncertainty in the patient specific parameters of a tumor growth model. To do so, we estimate the posterior probability of the parameters, which allows us to analyze the correlation between the personalized parameters of the model. Moreover, it gives valuable information on the confidence one has in the estimation and the ability of the model to explain the data. The method is based on the Bayesian personalization of a tumor growth model. We specifically apply it to the personalization of glioblastoma growth using a reaction-diffusion model. We detail the results on one synthetic -but realistic -case, and 7 patients.

B. Clinical Background

Gliomas account for 30% of primary brain and central nervous system tumors. They are characterized by their infiltrative nature: malignant cells invade the tissue by progressing along the white matter fiber tracts. They can either be low grade or high grade, and very well localized or extremely diffused. We focus on grade IV gliomas called glioblastomas multiformes (GBM). They are the most common glioma -at least 50% of the recorded cases, and an incidence of 3.19 per 100,000 [START_REF] Dolecek | CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United States in 2005-2009[END_REF]. They grow extremely rapidly, and are the most malignant primary brain tumors.

C. Modeling: Previous Work

The modeling of the tumor evolution can be done at different scales [START_REF] Angelini | Glioma dynamics and computational models: a review of segmentation, registration, and in silico growth algorithms and their clinical applications[END_REF]. On the microscopic scale, the model describes the evolution of individual cells based on division and invasion rules. On the macroscopic scale, the model describes the spatio-temporal evolution of the tumor cell density in the considered domain, based on partial differential equations. Over the last 20 years, particular attention has been given to the reaction-diffusion equation, (1)

D∇u. - → n ∂Ω = 0 (2) θ = (D, ρ) (3)
Equation ( 1) describes the spatio-temporal evolution of the tumor cell density u, which infiltrates neighboring tissues with a diffusion tensor D, and proliferates with a net proliferation rate ρ. Equation ( 2) enforces Neumann boundary conditions on the brain domain Ω, and θ denotes the parameters of interest of the model. In 1989, Murray et al. [9] pioneered on the use of mathematics applied to brain tumor growth, using the reaction-diffusion equation with an exponential growth. Assuming an isotropic diffusion in an infinite domain, he relates the velocity v of the tumor growth to the proliferation rate ρ and the diffusion coefficient D, v = 2 √ ρD. This work has lead the way to numerous other studies [START_REF] Tracqui | A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth[END_REF], [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF], [START_REF] Chaplain | Avascular growth, angiogenesis and vascular growth in solid tumours: The mathematical modelling of the stages of tumour development[END_REF].

Improving the computational efficiency of tumor growth models has been the topic of several studies. In 2010, Konukoglu et al. [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] introduced an approximation of the asymptotic solution of the reaction-diffusion equation based on an Eikonal equation, which can be solved with a Fast Marching algorithm, resulting in drastically reduced computation times. In 2012, Mosayebi et al. [START_REF] Mosayebi | Tumor invasion margin on the riemannian space of brain fibers[END_REF] introduced a similar method to compute the tumor invasion in the brain parenchyma based on a geodesic distance computed from DTI information.

The mass effect applied by the tumor to the surrounding tissues has been modeled in 2005 by Clatz et al. [START_REF] Clatz | Realistic simulation of the 3-D growth of brain tumors in MR images coupling diffusion with biomechanical deformation[END_REF]. They coupled the reaction-diffusion equation with a mechanical model, to take into account the mass effect, i.e. the displacement of adjacent brain tissue due to the growing tumor. They relate the diffusion tensor D to the Diffusion Tensor Images (DTI), thus forcing the tumor cells to preferably follow the white matter fiber tracts. In 2008, Hogea et al. [START_REF] Hogea | An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects[END_REF] proposed another approach to include mass effect, using a reactiondiffusion-advection model coupled with a Eulerian mechanical framework.

Finally, research efforts recently targeted the development of more complex models. In 2011, Swanson et al. [START_REF] Swanson | Quantifying the role of angiogenesis in malignant progression of gliomas: In silico modeling integrates imaging and histology[END_REF] developed a model dividing the tumor cells into sub-categories: the proliferative, the invasive and the hypoxic cells. This type of complex models has been further developed by Saut et al. [START_REF] Saut | A multilayer grow-or-go model for gbm: effects of invasive cells and anti-angiogenesis on growth[END_REF] in 2014, in which they developed a coupled PDE-based model to simulate the growth of high grade gliomas seeded in real images. Finally, the models were extended to take into account the apparition of edema during the growth [START_REF] Hawkins-Daarud | Modeling tumor-associated edema in gliomas during anti-angiogenic therapy and its impact on imageable tumor[END_REF], [START_REF] Badoual | Oedema-based model for diffuse low-grade gliomas: application to clinical cases under radiotherapy[END_REF].

D. Personalization: Previous Work

It is of particular interest to solve the inverse problem of estimating the model's parameter θ = (D, ρ), which can describe a specific dynamic observed for a patient. As such, the modeling effort goes hand in hand with the personalization work. For personalization, most methods rely on comparing the model estimation with the abnormalities visible on Magnetic Resonance Images (MRIs) at two different time points.

A common approach to personalize reaction-diffusion equations is based on the observation that they admit wave-like solutions whose parameters can be related to the volume of the abnormalities observed on MRIs. In 2007, Harpold et al. [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF] described the personalization of a glioma growth model using the reaction-diffusion equation with an exponential growth.

Another popular approach relies on solving the inverse problem using an optimization strategy. Hogea et al. [START_REF] Hogea | An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects[END_REF] formulated a PDE-constrained optimization problem to estimate the parameters of the model based on image registration. Konukoglu et al. [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] used the derivative-free optimization algorithm BOBYQA (for bounded optimization by quadratic approximation) [START_REF] Powell | The BOBYQA algorithm for bound constrained optimization without derivatives[END_REF] to minimize the distance between the tumor segmentations observed on the MRI, and the output of the model. In 2012, Gooya et al. [START_REF] Gooya | GLISTR: glioma image segmentation and registration[END_REF] proposed a method for the segmentation and registration of MRIs presenting glioblastoma based on the personalization of a reaction-diffusion-advection model using a single time point acquisition.

Finally, probabilistic approaches recently gained interest to solve complex inverse problems. Menze et al. [START_REF] Menze | A generative approach for image-based modeling of tumor growth[END_REF] was among the first to propose a Bayesian formulation for brain tumor growth personalization, based on the approximation of the posterior using sparse grids. It is interesting to note that Bayesian formulations have been used for a few years in the field of cardiac modeling [START_REF] Konukoglu | Efficient probabilistic model personalization integrating uncertainty on data and parameters: Application to eikonal-diffusion models in cardiac electrophysiology[END_REF], [START_REF] Neumann | Robust image-based estimation of cardiac tissue parameters and their uncertainty from noisy data[END_REF]. However, the methods usually rely on approximations of the forward model using reduced order model, such as the polynomial chaos, to make the estimation of the posterior using Markov Chain Monte Carlo (MCMC) computationally tractable.

E. Contributions of the Paper

In this paper, we propose a Bayesian method for the personalization of a tumor growth model based on the reactiondiffusion equation. We propose the use of the Lattice Boltzmann Method (LBM) to implement the tumor growth model which results in reduced computation times. This is combined with a high acceptance rate Monte Carlo technique called the Gaussian Process Hamiltonian Monte Carlo (GPHMC). Contrary to previous approaches, our method does not rely on approximations of the forward model (resp. posterior probability) using reduced order models [START_REF] Konukoglu | Efficient probabilistic model personalization integrating uncertainty on data and parameters: Application to eikonal-diffusion models in cardiac electrophysiology[END_REF], [START_REF] Neumann | Robust image-based estimation of cardiac tissue parameters and their uncertainty from noisy data[END_REF] (resp. sparse grid methods [START_REF] Menze | A generative approach for image-based modeling of tumor growth[END_REF]). We compare our approach to two methods adapted from the literature. The former is based on the spherical asymptotic analysis of the forward model, inspired by the work of Swanson et al. [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF], [START_REF] Corwin | Toward patient-specific, biologically optimized radiation therapy plans for the treatment of glioblastoma[END_REF]. The latter is based on the gradient-free optimization method BOBYQA, and is used in the work of Konukoglu et al. [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF].

This paper extends [START_REF] Lê | Bayesian personalization of brain tumor growth model[END_REF] with a comparison with a spherical asymptotic analysis of the personalization, and more com-prehensive analysis on 3 additional patients. Also, a new likelihood model based on the 95th percentile of the Hausdorff distance is used, as well as a new log-uniform prior for the parameters of interest. Finally, the parameters of the GPHMC have been updated to increase the robustness of the personalization: i) the initialization of the Gaussian process is now done with a coarse grid to ensure that the whole space of parameters is covered, ii) the parameters of the Gaussian process are set by maximizing the likelihood, iii) the noise level of the likelihood model has been reduced from σ = 10 mm in [START_REF] Lê | Bayesian personalization of brain tumor growth model[END_REF] to σ = 5 mm in order to increase the focus on the region of interest (i.e. the posterior is more peaked when the noise level is lower).

II. MODEL

A. The Reaction-Diffusion Model

Because of the logistic reaction term, equation (1) does not have a closed form solution. However, in an infinite domain, with constant proliferation rate and diffusion coefficient, equation (1) admits solutions which asymptotically behave like traveling waves with speed v = 2 √ ρD [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF]. This solution has the asymptotic form u(x, t) = u(x -vt) = u(ζ). Plugging it into equation (1), we get the ordinary differential equation

n Dn d 2 u dζ 2 + v du dζ + ρu(1 -u) = 0 (4)
where n is the direction of motion in the infinite domain. Moreover, linearizing equation ( 4) for small u << 1, meaning that u(1 -u) ∼ u, we get the second order linear equation

n Dn d 2 u dζ 2 + v du dζ + ρu = 0 (5)
which admits solutions of the form u(ζ) = (Aζ + B) exp(-x/λ) with A, B two constants, and λ = D/ρ. This equation shows that the parameter λ is of particular importance: it is related to the spatial decay of the tumor cell density, and it is referred to as the invisibility index [START_REF] Corwin | Toward patient-specific, biologically optimized radiation therapy plans for the treatment of glioblastoma[END_REF].

In the case of GBM growth modeling, it is clinically admitted that tumor cells have higher motility in white matter, compared to gray matter. Some work has been conducted to relate the diffusion tensor D to DTI -see [START_REF] Stretton | Importance of patient DTI's to accurately model glioma growth using the reaction diffusion equation[END_REF], [START_REF] Dittmann | Use of diffusion tensor images in glioma growth modeling for radiotherapy target delineation[END_REF] for a detailed discussion. In this work, for simplification, we follow [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF], and define the diffusion tensor as D = d w I in the white matter, and D = d w /10 I in the gray matter, where I is the 3x3 identity matrix, and d w a scalar parametrizing the diffusion tensor. We further identify the parameter D with d w . As such, the diffusion is heterogeneous and locally isotropic. This model reproduces the infiltrative nature of the GBM, takes into account anatomical barriers (ventricles, sulci, falx cerebri), and the tumor's preferential progression along white matter tracts such as the corpus callosum.

B. Model and MRIs

Data driven GBM growth modeling is based on the acquisition of different MRI modalities routinely acquired: T1 with Gadolinium contrast agent (T1Gd) and T2 Fluid-Attenuated Inversion Recovery (T2-FLAIR). They reveal different parts of the tumor: the active part, or proliferative rim can be observed on the T1Gd, while the T2-FLAIR reveals the disruption of the extacellular matrix, or edema. The T1Gd abnormality is usually included inside the larger T2-FLAIR abnormality (Figure 1). The personalization of the growth of the tumor is based on two consecutive time points, resulting in a total of four abnormalities segmented by the clinician.

We can immediately see one of the pitfalls of most models: the tumor cell density is not directly visible on T1Gd and T2-FLAIR MRIs. They rather reveal the impact of the tumor growth on the brain. In order to relate the tumor cell density u -solutions of (1) -to the MRIs, the frontier of the visible abnormalities is assumed to correspond to a threshold value of the tumor cell density u. We note τ 1 the value of the tumor cell density u corresponding to the frontier of the T1Gd abnormality, and τ 2 the value corresponding to the frontier of the T2-FLAIR abnormality. The invisibility index defined in the previous section is related to the distance between the boundaries of the T1Gd abnormality and the T2-FLAIR abnormality: the larger the distance, the larger the invisibility index is.

III. SIMULATION

A. Preprocessing of the data

The pre-processing of the MRIs is of particular importance (Figure 2). The T1Gd and T2-FLAIR abnormalities were segmented by a clinician, and the fractional anisotropy (FA) MRI was extracted from the DTI. First, the segmentations of the second time point t 2 are mapped on the space of the first time point t 1 as follows. i) For each visit, the T2-FLAIR and FA MRIs are rigidly registered to the T1Gd MRI [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF]. ii) The T1Gd MRI at t 2 is non-linearly registered to the T1Gd MRI at t 1 using the FNIRT function of the FSL software [START_REF] Andersson | Non-linear registration, aka spatial normalisation fmrib technical report tr07ja2[END_REF], [START_REF] Jenkinson | FSL[END_REF]. The segmentation of the tumor at t 2 is used to exclude the tumor from the similarity criterion (sum-of-squared difference) during the registration. As such, the intensity values in the tumor are ignored, and the tumor is warped in accordance with the surrounding (non-masked) tissues. iii) The resulting transformations are applied to transport the T2-FLAIR abnormalities at t 1 and t 2 , and the T1Gd abnormality at t 2 , on the T1Gd MRI at t 1 . Second, white and gray matter, and cerebrospinal fluid (CSF) are extracted from MRIs at t 1 . The pipeline is the following i) Extract the brain from the skull [START_REF] Iglesias | Robust brain extraction across datasets and comparison with publicly available methods[END_REF]. ii) Correct for the bias and segment the white matter, gray matter, and CSF [START_REF] Zhang | Segmentation of brain mr images through a hidden markov random field model and the expectationmaximization algorithm[END_REF]. iii) Separate the left and right hemispheres [START_REF] Zhao | Automatic cerebral and cerebellar hemisphere segmentation in 3D MRI: adaptive disconnection algorithm[END_REF], iv) The voxels at the boundary between the two hemispheres are tagged as CSF, in order to prevent the tumor from invading the contra-lateral hemisphere through the falx cerebri. The voxels with high FA value (> 0.45) and which were tagged as white matter are not affected by this process in order to ensure that the corpus callosum stays segmented as white matter. The red circle on Figure 2 emphasizes the importance of the hemisphere separation to label the falx cerebri as CSF, and the corpus callosum as white matter.

B. Lattice Boltzmann Method

A typical approach to implement the reaction-diffusion equation is to discretize the equation using the Crank-Nicolson scheme [START_REF] Özugurlu | A note on the numerical approach for the reactiondiffusion problem to model the density of the tumor growth dynamics[END_REF]. This requires the inversion of a large sparse matrix n × n where n is the number of voxels in the image, using a preconditioned gradient method like the biconjugate gradient stabilized method. For 3D MRIs with n ∼ 10 6 , this approach is computationally expensive. In this paper, we use the more recent explicit method called the Lattice Boltzmann Method (LBM). LBM has been successfully applied to implement the reaction-diffusion equation in the fields of cardiac electrophysiology modeling [START_REF] Rapaka | LBM-EP: lattice-Boltzmann method for fast cardiac electrophysiology simulation from 3D images[END_REF], [START_REF] Neumann | Robust image-based estimation of cardiac tissue parameters and their uncertainty from noisy data[END_REF], and liver tumor radiofrequency-ablation [START_REF] Audigier | Efficient lattice Boltzmann solver for patient-specific radiofrequency ablation of hepatic tumors[END_REF]. The idea is to model the reaction-diffusion equation as a set of fictitious particles which collide and stream on the cartesian grid. Using the D3Q7 scheme (3 dimensions and 7 velocities), the solution of the reaction-diffusion equation ( 1) is discretized u = 6 α=0 u α , where α indexes the 7 different velocity directions {e α } 6 α=0 in the 3D cartesian grid (i.e. the considered voxel and its six immediate neighbors).

[e 0 , e 1 , e 2 , e 3 , e 4 , e 5 , e 6 , ] =   0 1 -1 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 1 -1   Furthermore, according to the connectivity of the grid, a weight vector {w α } 6 α=0 = [1/4, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8] is defined . The algorithm, as described in [START_REF] Yoshida | Multiple-Relaxation-Time LBM for the convection and anisotropic diffusion equation[END_REF], consists of 4 steps: i) Initialization, ii) Collision, iii) Streaming, and iv) Neumann boundary conditions. During the initialization, we set u α = w α u, for every α. The collisions of the fictitious particles is then modeled with

∀α u * α = u α - 6 i=0 A α,i u i + ∆t w α ρu(1 -u)
where A = M -1 SM ∈ R 7×7 is the collision matrix. M projects a vector on the moment space, and S is the relaxation time matrix,

M =           1 1 1 1 1 1 1 0 1 -1 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 1 -1 6 -1 -1 -1 -1 -1 -1 0 2 2 -1 -1 -1 -1 0 0 0 1 1 -1 -1           S -1 =           τ 0 0 0 0 0 0 0 0 τ xx τ xy τ xz 0 0 0 0 τ xy τ yy τ yz 0 0 0 0 τ xz τ yz τ zz 0 0 0 0 0 0 0 τ 4 0 0 0 0 0 0 0 τ 5 0 0 0 0 0 0 0 τ 6           with τ ij = 1/2δ ij + 4∆t/(∆x 2 )D ij , τ k = 1.
33, and D is the diffusion tensor. The value of τ k does not impact the convergence of the solution but rather the stability of the method. A value of 1.33 has been found to lead to reasonable results [START_REF] Rapaka | LBM-EP: lattice-Boltzmann method for fast cardiac electrophysiology simulation from 3D images[END_REF]. During the streaming step, for every α, we have u α (x + δx e α , t + δt) = u * α (x, t). Finally, we enforce the Neumann boundary conditions with the equations

u α(x, t + δt) = 1 2∆ u α (x, t) + 2∆ -1 2∆ u α(x, t) (6)
u α(x, t + δt) =2∆u α (x, t) + (1 -2∆)u α (x -∆x e α , t) (7)
where α refers to the velocity of opposite directions than α, and ∆ is the distance between the fictitious particle and the boundary [START_REF] Yu | Viscous flow computations with the method of lattice Boltzmann equation[END_REF] (see Figure 3). We use equation ( 6) if ∆ > 1/2 and equation ( 7) otherwise. With a diffusion coefficient constant in time, the LBM does not require costly inversion of matrices. The LBM is easily parallelized such that simulating 30 days of growth, with δt = 0.1 day, takes approximately 50 seconds on a 2.3 Ghz 50 core machine for a 1 mm isotropic 155 × 182 × 157 grid.

C. Initialization

The initialization of the tumor cell density u(t = t 1 , x) at the time of the first acquisition is of particular importance, as it impacts the rest of the simulation. In this work, the tumor tail extrapolation algorithm described in [START_REF] Konukoglu | Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins[END_REF] is used. The tumor cell density is computed outward (and inward) of the T1Gd abnormality borders as a static approximation of the wave-like solution of equation ( 1) with parameter θ. The algorithm is based on the recursive approximation of equation ( 4) for current values of the solution u. It only depends on the invisibility index λ = D/ρ. The details of the implementation can be found in [START_REF] Konukoglu | Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins[END_REF]. By construction Figure 3: Cartesian grid with a boundary. Picture adapted from [START_REF] Yu | Viscous flow computations with the method of lattice Boltzmann equation[END_REF] of the initialization, the T1Gd abnormality falls exactly on the threshold τ 1 of the tumor cell density at the first time point (Figure 1).

IV. PERSONALIZATION

We describe in this section the personalization of the parameter θ = (D, ρ) of the reaction-diffusion equation [START_REF] Raman | Computational trials: Unraveling motility phenotypes, progression patterns, and treatment options for glioblastoma multiforme[END_REF]. Three different methods with increasing levels of complexity are presented. The first two are adapted from the literature, and are compared to the last one which we propose for the Bayesian personalization.

The first method (Section IV-A) is based on Swanson's approach [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF], [START_REF] Corwin | Toward patient-specific, biologically optimized radiation therapy plans for the treatment of glioblastoma[END_REF]. The spherical asymptotic analysis of the solution in a large domain is used to relate the invisibility index to the measured radii of the T1Gd and T2-FLAIR abnormalities, using multiple runs of the forward model.

The second method (Section IV-B) is based on [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] where the derivative-free optimization method BOBYQA is used to infer the parameters D and ρ. In this article, we apply it to the case of glioblastoma where two abnormalities are observable at each time point. We do not optimize for the source point of the tumor since we initialize the model with the contour of the first T1Gd abnormality.

Finally, we propose to estimate the posterior probability of the parameters D and ρ, knowing the clinician segmentations, using a MCMC method, the Gaussian Process Hamiltonian Monte Carlo.

A. Spherical Asymptotic Analysis

The parameters D and ρ can be related to the asymptotic velocity v = 2 √ Dρ, and the invisibility index λ = D/ρ. The knowledge of the velocity and invisibility index uniquely identifies D and ρ since D = vλ/2 and ρ = v/(2λ). As such, given a patient and the segmented T1Gd and T2-FLAIR abnormalities, we want to infer the parameters (v, λ) from the measured radius of the abnormality through their measured volume. We hypothesize that there is a simple relationship between those radii, and the velocity and the invisibility index [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF]. We consider the special case of an isotropic and homogeneous domain with a Gaussian initialization. The speed of growth can be measured as the temporal variation of the radial expansion of the abnormalities over time (Figure 4 left). Since the distinct speed of growth of the T1Gd and T2-FLAIR abnormalities can be measured (and can be different), we use the geometric mean between the two, v = √ v 1 v 2 , as an estimate. We propose in this section to find a simple relationship between the invisibility index, and the measured radii of the T1Gd and T2-FLAIR abnormalities (see Figure 4). The tumor growth model was run on a large 201×201×201 1 mm isotropic grid, with different parameters D and ρ, for 200 days, and initialized with a symmetric Gaussian. Using 15 equally spaced values of D ∈ [0.02 1.5] mm 2 .day -1 , and 15 equally spaced values of ρ ∈ [0.002 0.2] day -1 , 225 simulations were performed, keeping 11 time points per simulations. For each time point, we considered the value of the radii of the T1Gd and T2-FLAIR abnormality (resp. r T2-FLAIR and r T1Gd ), using thresholds of τ 1 = 80% and τ 2 = 16%. We observe a good linear relationship between the invisibility index λ, and the radius difference (r T2-FLAIR -r T1Gd ) (Figure 5),

λ = D ρ = a (r T2-FLAIR -r T1Gd ) + b (8)
with a = 0.13 and b = 0.23 mm for D in mm 2 .day -1 and ρ in day -1 . We further checked that the measured velocity corresponded to the asymptotic one (Figure 5). In this paper, we consider the personalization using two time points. Since the invisibility index can be measured on both time points, similarly to the speed, we use the geometric mean between the two as an estimate. This personalization is based on volumetric consideration, and valid in large homogeneous domain, much like Swanson's approach [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF], [START_REF] Corwin | Toward patient-specific, biologically optimized radiation therapy plans for the treatment of glioblastoma[END_REF].

B. BOBYQA Optimization

The previous method does not take into account the inhomogeneity and anisotropy of the growth, and does not account for anatomical barriers such as CSF. In order to derive a finer estimation, we use an optimization method. We need to define the error corresponding to a simulation with associated parameter θ. The simulation is initialized at t 1 using the first T1Gd abnormality, and simulated with the LBM until the second time point t 2 is reached. We extract the simulated contours corresponding to the thresholds τ 1 and τ 2 , at t 1 and t 2 . The 95th percentile of the symmetric Hausdorff distance between the borders of the clinician segmentation and the extracted contours is computed for: i) the T2-FLAIR abnormality at time t 1 , ii) the T2-FLAIR abnormality at time t 2 , iii) the T1Gd abnormality at time t 2 . The mean of these distances H mean is used as an error measure for the simulation. We use this error measure because the Hausdorff distance is sensitive and independent to the size of the abnormality. As such, the T1Gd and T2-FLAIR will be penalized in a similar fashion. Other error measures could be considered such as the DICE coefficient, which is however sensitive to the size of the abnormalities. We minimize the error using the derivative-free optimization algorithm BOBYQA [START_REF] Powell | The BOBYQA algorithm for bound constrained optimization without derivatives[END_REF], using the implementation of the library NLopt 1 . The algorithm is run 9 times with 9 different initializations to explore various local minima. We keep the best of the 9 solutions. This personalization is in the spirit of Konukoglu's work [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF], but applied to glioblastomas.

C. Bayesian Personalization

We denote by S the set of clinician segmentations the model should fit. We are interested in the posterior probability of the model parameter θ = (D, ρ), knowing the observations S. To cast the problem in a probabilistic framework, we follow Bayes rule: P (θ|S) ∝ P (S|θ) P (θ). The likelihood is modeled as P (S|θ) ∝ exp(-H 2 mean /σ 2 ), where the distance H mean is the mean of the 95th percentile symmetric Hausdorff distance between the border of the segmentations S and the isolines of the simulated tumor cell density u using θ, and the thresholds τ 1 and τ 2 . As such, the negative log-likelihood H 2 mean /σ 2 is the error term optimized during BOBYQA, normalized with the noise level σ. P (θ) is the prior on the parameters of the model. We want to estimate the posterior distribution P (θ|S). To do so, samples are drawn from the posterior probability using the Gaussian Process Hamiltonian Monte Carlo [START_REF] Rasmussen | Gaussian processes to speed up hybrid Monte Carlo for expensive Bayesian integrals[END_REF].

Hamiltonian Monte Carlo (HMC). HMC is a Markov Chain Monte Carlo algorithm which uses a refined proposal density function based on the Hamiltonian dynamics [START_REF] Neal | MCMC using Hamiltonian dynamics[END_REF]. The idea is to have a high acceptance rate while proposing points relatively far from the current point. The problem is augmented with a momentum variable p ∼ N (0, I). By randomly sampling p, we define a current state (θ, p). The energy of the state is H(θ, p) = E pot + E kin , with potential energy E pot = -log(P (S|θ)P (θ)), and kinetic energy ]. The conservation of the energy during the Hamiltonian dynamics -up to the numerical discretization accuracy of the Leapfrog scheme -insures a high acceptance rate A, which is the ratio of proposed samples which are accepted. The boundary conditions on the bounded parameter space θ are enforced using a bounce back condition [START_REF] Betancourt | Nested sampling with constrained hamiltonian monte carlo[END_REF] during the Leapfrog scheme used for the time integration. More precisely, if during the computations, the parameter θ crosses a boundary, its moment p is reversed and projected on the normal of the boundary.

Gaussian Process Hamiltonian Monte Carlo (GPHMC). In the HMC, computing the Hamiltonian dynamics -equations ( 9) -requires a significant amount of model evalutations. To circumvent this difficulty, E pot is approximated with a Gaussian process [START_REF] Rasmussen | Gaussian processes to speed up hybrid Monte Carlo for expensive Bayesian integrals[END_REF]. During the initialization phase, the forward model is evaluated on a coarse grid to initialize the Gaussian process. During the exploration phase, the forward model is evaluated at locations of low E pot and high uncertainty on the Gaussian process interpolation (details can be found in [START_REF] Rasmussen | Gaussian processes to speed up hybrid Monte Carlo for expensive Bayesian integrals[END_REF]). HMC is then run using the Gaussian process interpolation of E pot to compute the Hamiltonian dynamics. Given that the Gaussian process well captures E pot , the GPHMC benefits from the high acceptance rate of the HMC, with far less model evaluations.

Parameters. The parameter θ is constrained such that D ∈ [10 -4 , 10] mm 2 /days, and ρ ∈ [10 -5 , 10] days -1 , which encloses expected values for glioblastomas [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF]. The prior P (θ) is assumed log-uniform within this bounded box. The Hamiltonian dynamics is run for ∆t H = 60 days. For the likelihood, the noise level σ has been empirically set to 5 mm to provide a reasonable acceptance rate. The noise level influences how peaked the estimated posterior will be. For the GPHMC, the Gaussian process is defined with a squared exponential covariance matrix

C(θ 1 , θ 2 ) = w 0 exp(-||θ 1 -θ 2 )|| 2 2 /w 2 1 ) [45]
, where ||.|| 2 is the l 2 norm. These parameters are estimated with maximum likelihood [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. During the GPHMC, the initialization is done with 49 forward model evaluations, while 50 evaluations are used for the exploratory phase, then 1000 samples are generated.

V. RESULTS

A. Synthetic Case

We perform a thorough analysis of a synthetic -but realistic -case. A simulation is run on the MNI atlas [START_REF] Mazziotta | A probabilistic atlas and reference system for the human brain: International consortium for brain mapping (icbm)[END_REF], with parameters D = 1 mm 2 .day -1 and ρ = 0.18 day -1 , for 30 days (Figure 6). The result of the 3 different personalization methods can be seen on Figure 8. The result of the BOBYQA optimization is close to the mode of the computed posterior density. However, the spherical asymptotic analysis personalization is largely under-estimating the diffusion parameter D with a solution of D = 0.41 mm 2 .day -1 , when D = 1.0 and 0.96 mm 2 .day -1 for the BOBYQA and the Maximum a Posteriori (MAP) sample of the Bayesian solution respectively. The asymptotic personalization method makes the assumption that the T1Gd and T2-FLAIR abnormalities are growing concentrically in an infinite domain. The presence of boundaries in realistic growth leads to the under-estimation of the difference of the radius of the abnormalities, leading to the under-estimation of the infiltration length (Figure 6).The acceptance rate during the GPHMC is 83%, reflecting the fact that the Gaussian process is a good interpolation of the potential energy E pot for the sampling.

The impact of the different parameters is analyzed on the synthetic case (Figure 8). We can see that using a uniform prior -like in [START_REF] Lê | Bayesian personalization of brain tumor growth model[END_REF] -results in samples which are more concentrated in the regions of high values of D and ρ. However, it does not change much the shape and the location of the mode. On the other hand, using a noise level σ = √ 5 mm instead of 5 mm (with a log-uniform prior) results in a posterior which is much more peaked around the same mode. Finally, using the median Hausdorff distance instead of the 95th percentile (with a log-uniform prior and σ = 5 mm) results in samples which are more spread because this distance is less discriminative. The effect is very similar to increasing the noise level σ. Finally, the method is compared to a simple grid evaluation of the potential energy for the synthetic case. The GPHMC automatically sample an extremely high density of points in the region of high probability, compared to the simplest grid approach (Figure 7).

B. Glioblastoma Patients

The method is applied to 7 patients. The T1Gd and T2-FLAIR abnormalities were segmented by a clinician. A summary of the results can be found on Table I. The patients come from diverse clinical studies. We selected patients who were not treated with resection, and where two time points with visible growth were available. The patients were under chemotherapy, radiotherapy, and other specific types of therapy such as anti-angiogenic drug. This is one of the reasons why we focus on personalization, and not on prediction. The complex therapy schedule makes the prediction of the future behavior of the tumor difficult.

The result of the estimation can be seen on Figure 9. The personalization based on the asymptotic analysis tends to provide low values of D and ρ. We observe the variety of behavior of the posterior for the different patients. Figure 9 shows the best BOBYQA solution. The second and third best solutions are only shown if they are distinguishable with the first one at this scale, or if they fall in regions of low potential energy. The best BOBYQA solution always falls close to a mode of the posterior, but does not reflect the correlation between the parameters D and ρ. Moreover, it is worth mentioning that for every patient, some of the differently initialized BOBYQA solutions lead to very poor results.

Patients 1, 2, and 3 present only one mode. The BOBYQA solution gives good results, while the asymptotic analysis always gives slightly off solutions. Patients 3, 6, and 7 present two modes. The two modes are very close to each other for patients 6 and 7 while they are far apart for patient 3. This reflects the fact that the model cannot explain simultaneously the four available segmentations, and different sets of parameters lead to equivalently plausible solutions. For patient 5, the BOBYQA solutions and spherical asymptotic analysis give completely false result. This seems to be due to the fact that the solutions lie near the boundary of the domain, which lead to false convergence of the BOBYQA algorithm. For this patient, the parameters tend to lie at the edge of the admissible domain defined in Section IV-C. The acceptance rate during the GPHMC is on average 72%. This high acceptance rate is an indication that the Gaussian process is a good interpolation of the potential energy. However, for some patients, the interpolation is not as good resulting in acceptance rate of 52% (Patient 3) and 56% (Patient 7). This is mainly due to narrow regions of high probability, where there are relatively few points to interpolate the potential energy. Note that the parameters of the GPHMC were the same for the seven different patients.

For a given patient, we can gather the samples in a matrix X ∈ R n×2 , where n is the number of samples. We can project the samples in the (log(λ), log(v)) space, and compute the empirical covariance C,

C = cov X 1 2 1 2 -1 2 1 2 (10)
We then compute the normalized variances σ inv and σ speed of the posterior samples along the speed and invisibility index axis respectively,

σ inv = C(1, 1) Tr [C] , σ speed = C(2, 2) Tr [C] (11)
where Tr [ . ] refers to the trace of the matrix. This way, we can quantify which of the invisibility index or speed is better captured by the model. For instance, patient 3 presents an extremely elongated posterior along the line of constant invisibility index (σ inv = 1%, σ speed = 99%), showing that the uncertainty on the fit is due to the speed rather than the invisibility index. For patient 1, the posterior is not as elongated, showing that the uncertainty due to the speed or the invisibility index are more similar (σ inv = 45%, σ speed = 55%). This is probably caused by the fact that the two time points are 105 days apart for patient 1. As such the speed of growth of the tumor is more identifiable resulting in a drop of the variance due to the speed compared to the other patients. Similarly, the synthetic case presents a larger growth which makes the identifiability of the speed easier.

We show the results of the simulation using the parameter θ corresponding to the Maximum a Posteriori (MAP) of the samples (Figure 10) for two representative patients. Furthermore, for each sample resulting from the GPHMC, one can compute the corresponding segmentation for the T1Gd at the second time point, and the T2-FLAIR at the first and second time points. We can then compute the probability for each voxel to lie within one of those segmentations. We show on Figure 11 the 10% (in blue) and 90% (in green) isolines of this probability mask which allows to visualize the uncertainty in the MAP segmentation.

VI. DISCUSSION

In this work, we presented an efficient implementation of the reaction-diffusion equation for the brain tumor growth based on the Lattice Boltzmann Method. We further presented estimation methods of the model's parameters of different levels of complexity. The simplest one is based on the asymptotic properties of the reaction-diffusion equation. It does not require complex computations, but fails when the growth is constrained by the brain boundaries. The second method is the optimization of an error term using a derivative-free algorithm. In our experiments, the method required on average 20 model evaluations per initialization, resulting in 180 model evaluation for a total of 9 different initializations. The third method is based on an efficient Monte Carlo method called the Gaussian Process Hamiltonian Monte Carlo, used to sample from a posterior, derived as the Boltzmann distribution of the previous error term. It requires a total of 1100 model evaluations. As such, the estimation of the posterior probability requires only 6 times more evaluations than the direct optimization, and provides additional valuable informations about the shape of the posterior. This reveals the possible presence of several modes, and the correlation between the parameters due to the lack of identifiability of the speed or the invisibility index.

Moreover, the samples of the posterior density could be used to estimate the density in the whole domain using for instance kernel density estimation [START_REF] Lê | Bayesian personalization of brain tumor growth model[END_REF]. This could be helpful to have access to other useful statistical indices such as the evidence of the model, which can help answer questions about model selection [START_REF] Oden | Selection and assessment of phenomenological models of tumor growth[END_REF], [START_REF] Lê | Bayesian personalization of brain tumor growth model[END_REF].

In the future, we could consider the thresholds τ 1 and τ 2 as parameters of the model to be sampled. This could also be the case for the noise level σ. Moreover, the modeling of the mass effect could be useful in order to avoid the non-linear registration between the second and the first time points.

VII. PERSPECTIVES

In the future, we intend to apply the Bayesian personalization in order to explicitly take into account the uncertainty in the expert's segmentation. More specifically, the segmentations used during each model evaluation could be sampled in the space of plausible segmentations [START_REF] Lê | GPSSI: Gaussian process for sampling segmentations of images[END_REF]. We also believe that this work could be used for automatic personalized therapy planning. Some work has already been done on relating tumor growth models to radiation response models to better define radiation therapy plans [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: improving target volume delineation[END_REF], [START_REF] Unkelbach | Radiotherapy planning for glioblastoma based on a tumor growth model: implications for spatial dose redistribution[END_REF], [START_REF] Rockne | A mathematical model for brain tumor response to radiation therapy[END_REF]. Such a method could provide personalized therapy plans taking into account the uncertainty in the model's parameters.

Figure 1 :

 1 Figure 1: (Left) The proliferative rim is outlined in orange on the T1Gd MRI at two different time points; (Middle Left) The edema is outlined in red on the T2-FLAIR MRI at two different time points. The edema encloses the proliferative rim; (Middle Right) Tumor cell density computed with the reaction-diffusion model. The black (resp. white) line is the threshold values τ1 (resp. τ2) corresponding to the T1Gd (resp. T2-FLAIR) abnormality; (Right) Comparison between the clinician segmentation and the contours from the model.

Figure 2 :

 2 Figure 2: First row: Registration pipeline at time t1 of the T1Gd, T2-FLAIR, and FA. Middle row: Registration pipeline at time t2 of the T1Gd, and T2-FLAIR. Bottom row, from left to right: segmentation of the left and right hemispheres to ensure the falx cerebri is labeled as CSF; gray matter segmentation; white matter segmentation.

Figure 4 :

 4 Figure 4: 1D graphical explanation detailing how the measured radii are related to the invisibility index and speed of growth of the tumor. (Left) Relationship between the invisibility index and the speed of the tumor, and the parameters D and ρ of the model. (Right) 1D tumor cell density at two different time points detailing how the measured radii are related to the parameters of the model.

Figure 5 :

 5 Figure 5: (Left) Example of simulations on a large cube with different parameters D and ρ used to infer the relationship between the invisibility index and the measured radii. (Middle) Invisibility index ( D/ρ) function of rT2-FLAIR -r T1Gd . The blue dots result from the LBM simulations in the square, the red line is the linear fit. (Right) Asymptotic speed (2 √ Dρ) function of the measured radial growth speed. The blue dots result from the LBM simulations, the red line is the first bisector.

Figure 6 :

 6 Figure 6: (Left) Synthetic growth of a tumor in the MNI atlas during 30 days with D = 1 mm 2 .day -1 and ρ = 0.18 day -1 . The T2-FLAIR abnormality (isodensity of τ2 = 16%) is outlined in red. The T1Gd abnormality (isodensity of τ1 = 80%) is outlined in orange. (Right) Limits of the asymptotic personalization in presence of boundaries. The infiltration length -proportional to rT2-FLAIR -r T1Gd -tends to be under-estimated by considering the equivalent spherical volumes.

Figure 7 :

 7 Figure 7: Comparison with a simple grid estimation of the posterior using 34 × 34 = 1156 model evaluations in the admissible domain for the synthetic case. The zoom to the region of low potential energy ( i.e. high probability) proves the value of the presented method. The color scale indicates the negative log-likelihood.

Figure 8 :

 8 Figure 8: In the black box, (Left) i) Initialization of the Gaussian process interpolating the potential energy Epot with a a coarse 7 × 7 grid ii) Refinement of the Gaussian process with points of low potential energy and high uncertainty (see details in [42]) (Right) Personalization of the tumor growth model for the synthetic case. The result of the asymptotic personalization is at the intersection of the dashed red lines, the result of the BOBYQA optimization is at the intersection of the full red lines, and falls on the true parameters values at this scale. The color scale indicates the negative log-likelihood. Outside of the black box, results of the GPHMC personalization for the synthetic case using (from left to right, top to bottom) a uniform prior, a noise level σ = √ 5, and the median Hausdorff distance.

Figure 9 :

 9 Figure 9: Personalization of the tumor growth model for the 7 patients. The result of the asymptotic personalization is at the intersection of the dashed red lines, the result of the BOBYQA optimization is at the intersection of the solid red lines. The presence of several BOBYQA solutions from different initializations is indicated when the best three solutions were noticeably distinct. The color scale indicates the negative log-likelihood of the posterior P (θ|S), which is equal to the potential energy Epot.

Figure 10 :

 10 Figure 10: Maximum a Posteriori fit for Patient 1 (Top) and Patient 4 (Bottom). For each patient: (Left) The proliferative rim is outlined in orange on the T1Gd MRI at two different time points; (Middle Left) The edema is outlined in red on the T2-FLAIR MRI at two different time points; (Middle Right) Tumor cell density, the black (resp. white) line is the threshold value τ1 (resp. τ2) corresponding to the T1Gd (resp. T2-FLAIR) abnormality; (Right) Comparison between the clinician segmentation and the contours from the model.

Figure 11 :

 11 Figure 11: (Top) Patient 2 (Bottom) Patient 4. From left to right: T1Gd for the second time point, T2-FLAIR for the first time point, T2-FLAIR for the second time point. The clinician segmentation are in orange for the T1Gd and red for the T2-FLAIR.The blue outline (resp. green) encloses the voxels which were present in at least 10% (resp. 90%) of the segmentations deduced from the samples.

Table I :

 I Results on 7 patients: time span between the time points; acceptance rate of the GPHMC; percentage of variance due to the speed; percentage of variance due to the invisibility index.

Asclepios Project, Inria Sophia Antipolis, France.

Martinos Center for Biomedical Imaging, Harvard-MIT Division of Health Sciences and Technology, Charlestown, MA, USA.

Department of Neurology, Massachusetts General Hospital, Boston, MA, USA.

Department of Radiation Oncology, Massachusetts General Hospital and Harvard Medical School, Boston, MA, USA.

Steven G. Johnson, The NLopt nonlinear-optimization package, http://abinitio.mit.edu/nlopt

ACKNOWLEDGMENT

Part of this work was funded by the European Research Council through the ERC Advanced Grant MedYMA 2011-291080.