Marc Bonnet
email: bonnet@lms.polytechnique.fr

Nicolas Nemitz
email: nemitz@lms.polytechnique.fr

FM-BEM and topological derivative applied to acoustic inverse scattering

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

Introduction

Defect identificaton problems are often solved by minimization of a cost function featuring the experimental data and (if available) prior information. Such cost functions are non-convex and exhibit local minima. Despite that fact, traditional iterative minimization or equation-solving methods are usually preferred to global search techniques such as evolutionary algorithms due to the prohibitive computational cost of solving large numbers of forward wave scattering problems. To perform optimally, gradient-based iterative algorithms are used in conjunction with shape sensitivity techniques see e.g. [START_REF] Kirsch | The domain derivative and two applications in inverse scattering theory[END_REF][START_REF] Nintcheu Fata | A computational basis for elastodynamic cavity identification in a semi-infinite solid[END_REF][START_REF] Potthast | On the convergence of a new Newton-type method in inverse scattering[END_REF].

Still, the stand-alone use of gradient-type minimization for such purposes is not satisfactory for its success is strongly dependent on a reliable prior information about the geometry of the hidden object. This has prompted the development of 'sampling' or 'probe' non-iterative methods, which may be used in isolation or as a preliminary step for choosing adequate initial guesses to be used in subsequent standard optimization schemes. Such methods are surveyed in a recent review article [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF] and include the linear sampling method [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF][START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF], not pursued here, and the application of topological sensitivity, which is the subject of this article. The concept of topological sensitivity consists in quantifying the sensitivity of a cost function with respect to the creation of a scatterer B ε (x s ) of small characteristic radius ε and given location x s , as a function of x s . It appeared first in [START_REF] Eschenauer | Bubble method for topology and shape optimization of structures[END_REF][START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Lochpositionierungskriterien[END_REF] in connection with topological optimization of mechanical structures, allowing to define algorithms where "excess" material is iteratively removed until a satisfactory shape and topology is reached [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF]. More recently, other investigations have studied the topological sensitivity as a sampling tool for inverse scattering problems, in the context of identification of cavities in 3D semi-infinite and infinite elastic media [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] and in elastic 3D bounded bodies [START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF], and of elastic inclusions [START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF] (see also [START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundsary integral equation[END_REF] for 2D elastostatics and [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] for 2D linear acoustics).

The distribution of topological sensitivity can be expressed in terms of a bilinear formula involving the free field and an adjoint field associated with the cost function. However, the computational cost of solving the forward and adjoint problems and evaluating the topological sensitivity distribution on a fine sampling grid increases rapidly with the non-dimensional wavenumber. The purpose of this article is to propose the topological sensitivity field computed by means of the Fast Multipole BEM (FM-BEM) [START_REF] Greengard | A fast algorithm for particle simulations[END_REF][START_REF] Gumerov | Fast multipole methods for the helmholtz equation in three dimensions[END_REF][START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF] as the basis of a computationally fast tool for probing acoustic media for hidden hard obstacles on the basis of overdetermined boundary data, within the model framework of forward scattering problems governed by the scalar Helmholtz equation. To that end, the FM-BEM is in particular applied to evaluate in a fast way the integral representation formulae expressing the free and adjoint fields at a large number of sampling points points inside the medium.

This article is organized as follows. After some preliminaries concerning the forward and inverse problems of interest (Section 2), the concept of topological sensitivity is presented in Section 3. The FM-BEM approach for the scalar Helmholtz equation is then summarized in Section 4. Finally, results of numerical experiments on qualitative scatterer identification using computed distributions of topological sensitivity are presented in Section 5, for configurations featuring O(10 5 ) nodal BE unknowns and O(10 6 ) sampling points.

Forward and inverse problems

This article is concerned with the identification of rigid obstacles embedded in acoustic media. The generic acoustic scattering problem of interest is defined as follows. Let Ω denote a three-dimensional open domain, either bounded or unbounded, with a sufficiently regular boundary S and filled with an acoustic medium characterized by wave velocity c and mass density ρ; this configuration will be referred to as the reference (i.e. obstacle-free) medium. Let B denote a rigid scatterer (or a set thereof) bounded by the closed surface Γ , so that Ω = Ω\ B is the acoustic region surrounding the scatterer. Steady-state excitations on S with angular frequency ω generate an acoustic pressure field u in the acoustic domain Ω , governed by the following set of field equations and boundary conditions (collectively denoted by P(B ) for later reference):

P(B ) : (∆ + k 2 )u = 0 (in Ω ), p[u ] = p D (on S), p[u ] = 0 (on Γ ), (1)
where k = ω/c is the wavenumber, w → p[w] ≡ w ,n = ∇w • n denotes the normal derivative operator, n is the normal on S ∪ Γ outward to Ω , and p D is the prescribed Neumann data over S (other types of boundary conditions may be considered as well). It is assumed that ω is not an eigenfrequency of any of the boundary-value problems arising in the ensuing developments.

In the inverse scattering problem of interest, an unknown obstacle B true , of boundary Γ true , is to be identified. The corresponding exact acoustic field u true is then governed by problem P(B true ). With reference to problem P(B true ), supplementary information is needed for the identification of B true . Here, measured values u obs of acoustic pressure are assumed to be available over the measurement surface S obs ⊂ S. Ideally (i.e. assuming that the physics is exactly described by the chosen linear acoustics setting and that no measurement errors are present), u obs is the trace of u true on S obs . The identification of B true may then be formulated in terms of the minimization of a cost function. Generic cost function of format

J (Ω ) = S obs ϕ u R (ξ), u I (ξ), ξ dΓ (2)
are considered, where u is the boundary trace of the solution to the forward problem P(B ) for an assumed obstacle configuration B , the subscripts 'R' and 'I' being used to indicate the real and imaginary parts of a complex quantity (i.e. w R = Re(w) and w I = Im(w)). For instance, the output leastsquares cost function associated to measurement u obs on S obs , commonly used for such purposes, corresponds to

ϕ(w R , w I , ξ) = 1 2 w(ξ)-u obs (ξ) 2 , (3)
The minimization of such cost functions can be performed using many methods, all of which are iterative and need repeated evaluations of J(Ω ). Traditional gradient-based minimization may converge within a moderate number of evaluations of J (Ω ) if the trial surface Γ can be described in terms of a few geometrical parameters, but reach a local minimum which depends on the choice of initial guess. Global search techniques, e.g. evolutionary algorithms [START_REF] Michalewicz | How to solve it: modern heuristics[END_REF] or sampling methods based on the Metropolis algorithm [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF], perform a global search (i.e. identify absolute and/or multiple minima), but at the cost of very large numbers of cost functions evaluations. In this article, the topological sensitivity is proposed as a tool for performing a qualitative global search at a computational cost which is far below that entailed by a true global optimization technique. The results may then (for example) be used as initial guesses in subsequent minimization-based inversion procedures.

3 Topological sensitivity of the cost function

Notations

Let B ε (x s ) = x s + εB, where B ⊂ R 3 is a fixed bounded open set with boundary S and volume |B| containing the origin, define the region of space occupied by a hard obstacle of (small) size ε > 0 containing a fixed sampling point x s . It is convenient to introduce the scaled position vector ξ defined by

ξ = x s + ε ξ (ξ ∈ B ε , ξ ∈ B) (4)
In particular, this mapping recasts integrals over B ε and Γ ε into integrals over B and S , respectively, and transforms the differential volume and area elements according to

dV ξ = ε 3 d Vξ (ξ ∈ B ε , ξ ∈ B), dΓ ξ = ε 2 d Γξ (ξ ∈ Γ ε , ξ ∈ S ) (5)
Without loss of generality, x s can be chosen as the center of B ε , i.e. such that

Bε (ξ -x s ) dV ξ = 0 , i.e. B ξ d Vξ = 0. (6)
Let u = u ε (ξ; x s ) denote the solution to the scattering problem P B ε (x s ) defined by [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF], where

Ω = Ω ε (x s ) = Ω \ B ε (x s ) and B ε (x s ) is the closure of B ε (x s ).
Further, let J(ε; x s ) be defined by

J(ε; x s ) = J Ω ε (x s ) = S obs ϕ u ε R (ξ), u ε I (ξ), ξ dΓ, (7)
For convenience, explicit references to x s will often be omitted in the sequel, e.g. by writing J(ε) or u ε (ξ) instead of J(ε; x s ) or u ε (ξ; x s ). The evaluation of J(ε) entails solving for u ε the forward problem P(B ε ). It is convenient, and customary, to decompose u ε according to

u ε = u + v ε , (8)
where u, the free field defined as the response of the obstacle-free (reference) medium Ω due to the given excitation p D , solves

(∆ + k 2 )u = 0 (in Ω), p[u] = p D (on S), (9)
while v ε , the scattered field, solves

(∆ + k 2 )v ε = 0 (in Ω ε ), p[v ε ] = 0 (on S), p[v ε ] = -p[u]
(on Γ ε ).

(10)

Expansion of J (ε)

To establish the topological sensitivity of J(ε), one starts with the expansion

J(ε) = J(0) + S obs Re ϕ ,u v ε dΓ + o |v ε | S obs , (11)
with

ϕ ,u ≡ ∂ϕ ∂u ε R -i ∂ϕ ∂u ε I u ε =u . ( 12
)
Let the adjoint field û be defined by

(∆ + k 2 )û = 0 (in Ω) , p[û] = ϕ N (on S obs ) , p[û] = 0 (on S \ S obs ). (13)
Then, the reciprocity identity (i.e. third Green's formula) applied to the states û and v ε over the domain Ω ε leads, by virtue of the boundary conditions in [START_REF] Greenbaum | Iterative methods for solving linear systems[END_REF] and [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], to the identity

S obs ϕ ,u v ε dΓ + Γε p[û]v ε dΓ + Γε ûp[u] dΓ = 0 (14)
As a result, the integral in the r.h.s. of ( 11) is converted into integrals over the vanishing cavity surface. Besides, since both u and û are also defined inside B ε , the last integral in ( 14) can be recast into a domain integral over B ε by means of the divergence formula. Expansion [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] then takes the form

J(ε) = J(0) + Re Bε ∇u •∇û -k 2 uû dV - Γε v ε p[û] dΓ + o(|v ε | S obs ), (15)
The first integral in [START_REF] Kirsch | The domain derivative and two applications in inverse scattering theory[END_REF] features a density whose definition does not depend on ε, and its expansion about ε = 0 can therefore be obtained by simply using the scaled coordinates (4), [START_REF] Darve | The fast multipole method: numerical implementation[END_REF] and expanding ∇u • ∇û -k 2 uû (x s + ε ξ).

In contrast, the second integral of (15) features the scattered field v ε , which depends on ε. Its asymptotic behaviour must then be obtained from that of v ε on Γ ε (taking into account the fact that Γ ε also depends on ε). This step is based on exploiting an integral equation reformulation of equations [START_REF] Greenbaum | Iterative methods for solving linear systems[END_REF].

Governing integral equation formulation for the scattered field

The governing problem [START_REF] Greenbaum | Iterative methods for solving linear systems[END_REF] for the scattered field v ε = v ε (• ; x s ) can be recast as the boundary integral equation [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF][START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF]:

1 2 v ε (x) + Γε H(x, ξ; k)v ε (ξ) dΓ ξ = - Γε G(x, ξ; k)p[u](ξ) dΓ ξ (x ∈ Γ ε ), ( 16
)
where the Green's function G(x, ξ, k) is defined by

(∆ ξ + k 2 )G(x, ξ, k) + δ(ξ -x) = 0 (ξ ∈ Ω), H(x, ξ, k) = 0 (ξ ∈ S), (17)
and

H(x, ξ; k) = ∇G(x, ξ, k) • n(ξ)
is the normal derivative of G(x, ξ, k) (the nabla symbol ∇, when used in front of such kernel, conventionally indicates a gradient with respect to the second argument ξ). Moreover, the free and adjoint fields have the explicit expressions

u(x) = S G(x, ξ; k)p D (ξ) dΓ ξ , û(x) = S G(x, ξ; k)ϕ ,u (ξ) dΓ ξ . ( 18
)
It is convenient for the present purposes to split (G, H) according to:

G(x, ξ; k) = G(x, ξ; k) + G C (x, ξ; k) H(x, ξ; k) = H(x, ξ; k) + H C (x, ξ; k), (19)
where (G, H) is the well-known singular free-space fundamental solution for the Helmholtz equation, given by

G(x, ξ; k) = 1 4πr e ikr , H(x, ξ; k) = [r • n(ξ)] ikr -1 4πr 3 e ikr , (20)
with r = ξx and r = |ξ -x| = |r|, and the complementary part (G C , H C ) is not singular at ξ = x. On using decomposition [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF], equation ( 16) becomes

1 2 f (x) + Γε H(x, ξ; k)v ε (ξ) dΓ ξ + Γε H C (x, ξ; k)v ε (ξ) dΓ ξ = - Γε G(x, ξ; k)p[u](ξ) dΓ ξ - Γε G C (x, ξ; k)p[u](ξ) dΓ ξ (x ∈ Γ ε ). ( 21
)

Leading asymptotic contribution to the scattered field

To study the asymptotic behaviour of integral equation ( 21) as ε → 0, it is useful to introduce further scaled geometric quantities:

x = εx , r = εr , r = εr (x, ξ ∈ Γ ε ; x, ξ ∈ S ) (22)
in addition to definition (4) of ξ. The leading contributions as ε → 0 to the fundamental kernels featured in equation ( 21) are

G(x, ξ; k) = ε -1 G(x, ξ) + O(1) H(x, ξ; k) = ε -2 H(x, ξ) + O(1) (x, ξ ∈ Γ ε ), (23)
for the singular kernels (G, H) defined by [START_REF] Potthast | On the convergence of a new Newton-type method in inverse scattering[END_REF], where

G(x, ξ) = 1 4πr , H(x, ξ) = - r • n( ξ) 4πr 3 (24)
are the fundamental kernels for the Laplace equation, and

G C (x, ξ; k) = G C (x s , x s ; k) + O(ε) H C (x, ξ; k) = H C (x s , x s ; k) + O(ε) (x, ξ ∈ Γ ε ). ( 25
)
for the nonsingular kernels (G C , H C ).

On performing the coordinate transformations ( 4), ( 22) and using estimates ( 23), [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF] together with (5b), one finds that

Γε H C (x, ξ; k)v ε (ξ) dΓ ξ = O(ε 2 |v ε |) Γε G C (x, ξ; k)p[u](ξ) dΓ ξ = O(ε 3 ). ( 26
)
The limiting form of integral equation ( 21) as ε → 0, retaining only the leading contributions, is hence found to be

LS v ε x s + ε ξ (x) = -ε∇u(x s ) • S 1 4πr n( ξ) d Γξ + o(ε) (27)
where LS , defined by

LS f (x) = 1 2 f (x) + S H(x, ξ)f ( ξ) d Γξ (x ∈ S ), (28)
is in fact the governing integral operator associated with exterior Neumann problems for the Laplace equation in the normalized domain R 3 \ B. Equation (27) indicates that the scattered field is of order O(ε) on Γ ε :

v ε (ξ) = v ε x s + ε ξ = ε∇u(x s ) • V( ξ) + o(ε) (ξ ∈ Γ ε , ξ ∈ S ) (29)
where the vector function V( ξ) solves the exterior Laplace problem

LS V (x) = -∇u(x s ) • S 1 4πr n( ξ) d Γξ (30) i.e. ∆ V = 0 ( ξ ∈ R 3 \ B), ∇ ξ V.n = -n ( ξ ∈ S ), |ξ| 2 V = O(1) (| ξ| → +∞) (31)
It is important to note that V( ξ) does not depend on the sampling point x s , and hence needs to be computed only once.

Topological derivative

On substituting (29) into (15) and taking (5) into account, one finally arrives at the following expansion of J(ε):

J(ε) = J(0) + ε 3 T (x s ) + o(ε 3 ) (32)
where T (x s ), the topological derivative of J(ε), is given in terms of the free and adjoint fields by

T (x s ) = Re ∇û • A(S ) •∇u -|B| k 2 ûu (x s ) (33)
and with the second-order tensor A(S ) defined by

A(S ) = |B| I - S V( ξ) ⊗ n( ξ) d Γξ (34) = |B| I - S n( ξ) ⊗ V( ξ) d Γξ
where the second equality (i.e. the fact that A(S ) is symmetric) easily stems from the third Green's identity and the definition (31) of V.

For arbitrary surfaces S which are sufficiently regular for integral equation (30) to be mathematically meaningful (this includes surfaces with edges and corners, e.g. box-shaped scatterers, but precludes infinitely-thin screens), the vector density V may be found by e.g. solving numerically three sets of BEM equations for exterior Laplace problems, a computationally modest task.

For the simplest case of a rigid spherical obstacle (where B is the unit sphere, |B| = 4π/3, and on which n( ξ) = -ξ) one easily finds by analytical means that V( ξ) = ξ/2. Then, (34) is readily found to be given by A(S ) = 2πI.

(35)

Practical computation of topological sensitivity

The developments thus far are based on the Green's function G defined by [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF], and lead to almost explicit formulae for T (x s ), their only non-explicit component being the auxiliary density V, which must be computed numerically except for simple shapes of the trial scatterer B.

In practice, this explicit character is retained only for geometries Ω such that the corresponding Green's function is known analytically. Such cases correspond with geometrically simple configurations, e.g. the acoustic halfspace. For configurations where the Green's function is not available, the free and adjoint fields and the nonsingular kernels may be sought as solutions of boundary integral equations [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF][START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF]. The free and adjoint fields, defined by ( 9) and ( 13), satisfy the well-known integral identities

c(x)u(x) + S H(x, ξ; k)u(ξ) dΓ ξ = S G(x, ξ; k)p D (ξ) dΓ ξ (36) c(x)û(x) + S H(x, ξ; k)û(ξ) dΓ ξ = - S obs G(x, ξ; k)ϕ N (ξ) dΓ ξ (37)
which provide integral equations for x ∈ S (with c(x) = 1/2, except at points x where S is only piecewise smooth, such as edges or corners, for which c(x) is also known) and integral representation formulae for x ∈ Ω (with c(x) = 1). Differentiation of the latter under the integral sign provide integral representations of ∇u(x) and ∇u(x s ) in Ω.

Hence, evaluation of the topological sensitivity field as given by (33) entails the solution of integral equations (36) and (37) for u and û on S, followed by an evaluation of u(x s ), ∇u(x s ) and û(x s ), ∇u(x s ) by means of (36) and (37) used as integral representations. The first step may involve large numbers of unknowns if the diameter of Ω spans more than a few wavelengths. Moreover, the present objective being to comprehensively explore a 3-D region of space for hidden scatterers by examining the distribution of T (x s ), formula (33) is to be evaluated at a large number of sampling points x s . Both steps can then be considerably accelerated using the FM-BEM. [START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF] The fast multipole method for Helmholtz equation

BEM discretization

Equations ( 36) and (37) are in this article solved by means of the simplest BEM discretization, which employs flat triangular boundary elements with straight edges and piecewise-linear C 0 interpolation of u and û (other choices, e.g. C 0 quadratic interpolation, would of course have been possible). All numerical results presented hereinafter have been obtained on that basis. The primary unknowns are the values of u or û at the mesh nodes, i.e. at all the element vertices. Equations (36) and (37) are collocated at all mesh nodes. All singular element integrals associated with the kernel H(x, ξ; k) are zero because of the piecewise-flat geometry representation, which of course simplifies the implementation. Denoting by N the total number of nodes (and hence of unknowns), this procedure gives rise to the linear systems of equations

[A]{u} = {b} (38) [A]{û} = { b} ( 39
)
where the N -vectors {u} and {û} collect all nodal values of u and û. The discussion to follow will focus on solving system (38), the adjoint system (39) being of course solved in exactly the same way with {b} replaced by { b}.

As the problem size N grows, direct solvers become impractical or infeasible with respect to both computing time and storage, mainly due to the fully-populated nature of the BEM matrix [A], and iterative solvers are used instead. Since [A] is a non-symmetric, invertible matrix (except when k is a eigenvalue for Ω and homogeneous Neumann BCs), the iterative solution technique most often used for systems such as (38) is the generalized minimal residual (GMRES) algorithm (see e.g. [START_REF] Greenbaum | Iterative methods for solving linear systems[END_REF]), which is applicable to general invertible square matrices. Such algorithms are based on matrix-vector evaluations, and therefore do not require actual assembly and storage of [A]. The GMRES algorithm requires repeated evaluations of the residual

{b} -[A]{u} ( 40
)
where {u} is a given solution candidate, which is updated at each GMRES iteration. Hence, one needs to compute (discretized versions of) the doublelayer and single-layer potentials featured in the left-hand and right-hand sides, respectively, of (36) and (37) for known densities.

Traditional BE methods lead to a O(N 2 ) computational cost for each residual evaluation, because element integrals computed for a collocation point cannot be reused for another collocation point. By adopting the Fast Multipole boundary element method (FM-BEM), each residual can be computed within a O(N LogN ) time. The implementation used here, concisely described in the remainder of this section, follows Darve [START_REF] Darve | The fast multipole method: numerical implementation[END_REF] and Sylvand [START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme : performances, parallélisation, applications[END_REF][START_REF] Sylvand | Complex industrial computations in electromagnetism using the fast multipole method[END_REF].

Expansion of the fundamental solution

The starting point for the FM-BEM is the following representation of the full-space fundamental solution [START_REF] Potthast | On the convergence of a new Newton-type method in inverse scattering[END_REF]:

G(x, ξ; k) = lim p→∞ ik 4π Ŝ e -ik(ŝ•x) T p (ŝ, r 0 )e ik(ŝ• ξ) dΓ ŝ (41)
where Ŝ = {ŝ, |ŝ| = 1} is the unit sphere, the position vector r = ξx has been decomposed as

r = (ξ 0 -x 0 ) + (ξ -ξ 0 ) -(x -x 0 ) = r 0 + ξ -x ( 42
)
in terms of two poles x 0 and ξ 0 , and with the transfer function T p (ŝ, r 0 ) defined by

T p (ŝ, r 0 ) = (-i) n 4π p n=0 (-1) n (2n + 1)h (1) n (kr 0 )P n (ŝ • r0 ) (43)
In (43), the h

n and P n are respectively the spherical Hankel functions of first kind and the Legendre polynomials, and r 0 ≡ |ξ 0 -x 0 |. Moreover, for any vector z ∈ R 3 , a hat symbol indicates the corresponding unit vector, i.e. ẑ = z/|z|. Representation (41) holds under the condition

|r 0 | > | ξ -x| (44)
The poles x 0 and ξ 0 are actually meant to be chosen close to ξ and x, respectively, so as to satisfy the stronger condition

|ξ -x 0 | > |x -x 0 | and |x -ξ 0 | > |ξ -ξ 0 | (45)
The representation of H(x, ξ; k) obtained by differentiating (41) is then

H(x, ξ; k) = lim p→∞ -k 2 4π Ŝ e -ik(ŝ•x) T p (ŝ, r 0 ) ŝ • n(ξ) e ik(ŝ• ξ) dΓ ŝ (46)
In practice, representations (41) and (46) are approximated by (i) using the transfer function T p (ŝ, r 0 ) with a finite index p and (ii) replacing the integral over the unit sphere with a quadrature rule with Q points ŝq and weights w q . So, one replaces (41) and ( 46) with the approximations

G(x, ξ; k) ≈ ik 4π Q q=1
w q e -ik(ŝq•x) T p (ŝ q , r 0 )e ik(ŝq• ξ) (47)

H(x, ξ; k) ≈ ik 4π Q q=1 w q e -ik(ŝq•x) T p (ŝ q , r 0 ) ŝq • n(ξ) e ik(ŝq• ξ) (48)
The choice of points ŝq and weights w q , and their number Q, depends on the truncation order p used in the transfer function T p . On parameterizing unit vectors ŝ ∈ Ŝ using spherical angular coordinates (θ, φ), a commonly used choice [START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme : performances, parallélisation, applications[END_REF] Now, let S x and S ξ denote two disjoint portions of S, and let the poles x 0 and ξ 0 be chosen close to S x and S ξ , respectively, in such a way that (45) holds for any x ∈ S x and ξ ∈ S ξ . Consider the computation of

I(x) = S ξ G(x, ξ; k)v(ξ) dΓ ξ (x ∈ S x )
for a given density v(ξ), which is a typical contribution to the evaluation of the residual (40). On substituting (47), one obtains

I(x) ≈ ik 4π Q(p) q=1
w q e -ik(ŝq•x) T p (ŝ q , r 0 )

S ξ e ik(ŝq• ξ) v(ξ) dΓ ξ (x ∈ S x ) (49)
So, the same integral over S ξ can be re-used for all collocation points x ∈ S x . Moreover, for a chosen portion S ξ , this is true for any portion S x and associated pole x 0 such that condition (45) holds. Computations of the form (49) can be decomposed into three stages: (i) compute for each quadrature point of Ŝ the multipole moment R(ŝ q ; ξ 0 ):

R(ŝ q ; ξ 0 ) = S ξ e ik(ŝq• ξ) v(ξ) dΓ ξ (50)
(ii) multiply R(ŝ q ; ξ 0 ) by the transfer function T p (ŝ q , r 0 ), to obtain local expansion coefficients L(ŝ q ; x 0 ) at x 0 :

L(ŝ q ; x 0 ) = T p (ŝ q , ξ 0 -x 0 )R(ŝ q ; ξ 0 ) (51)
(iii) for all x ∈ S x , transfer L(ŝ q ; x 0 ) locally from x 0 to x and perform the numerical quadrature over Ŝ, to obtain (an approximation of) I(x):

I(x) ≈ ik 4π Q(p) q=1
w q e -ik(ŝq•x) L(ŝ q ; x 0 ) (52)

The one-level fast multipole method consists in partitioning the spatial region containing S into cubic cells of identical sizes whose vertices lie on a regular cubic grid. Each pair (S x , S ξ ) is such that S x = S∩C x and S ξ = S∩C ξ , where (C x , C ξ ) is any pair of disjoint cubic cells. The poles x 0 and ξ 0 are the respective cell centroids. The one-level FM-BEM has a complexity of O(N 3/2 ) per iteration for equations of type (38), which is of course better than the O(N 2 ) complexity of traditional BEM, but not optimal. Further acceleration is provided by using the multi-level fast multipole method, where the size of clusters S x , S ξ depends of their distance.

Multilevel FM-BEM algorithm

To exploit optimally the acceleration afforded by (49), a hierarchical oct-tree structure of elements is introduced. For that purpose, a cube containing the boundary S, called 'level-0 cell', is divided into eight cubes (level-1 cells), each of which is divided in the same fashion, and so on. A level-cell is divided into level-( + 1) cells unless it contains less than a preset (relatively small) number E of boundary elements (such cells are termed leaves). A noteworthy feature of the FM-BEM applied to Helmholtz-type equations is that to achieve the same accuracy in approximations (47), (48) at all levels, the truncation parameter p is level-dependent. A often-used formula [START_REF] Darve | The fast multipole method: numerical implementation[END_REF] for the adjustment of p is of the form

p( ) = √ 3kd( ) + CLog 10 ( √ 3kd( ) + π) (53)
where d( ) is the size of a cubic -level cell (so, d( + 1) = 2d( ) and p( + 1) is roughly 2p( )) and C is a constant. As a result, the set of quadrature points ŝq on Ŝ is also level-dependent: each level necessitates a distinct set ŝ q of quadrature points and associated weights. For the present implementation, values of C such that 2 ≤ C ≤ 8 were found to provide a acceptable compromise between accuracy and cost, and C = 4 was actually used. The FM-BEM algorithm then consists of:

• An upward pass where multipole moments (50) are first computed for the lowest-level cells and then recursively aggregated by moving upward in the tree until level 2 (for which there are 4 × 4 × 4 cells) is reached. Letting C (ξ 0 ) denote the set of children of a given -level cell C(ξ 0 ), i.e. of ( +1)level cells C(ξ 0 ) contained in C(ξ 0 ), this operation relies on the identity

R(ŝ +1 q ; ξ 0 ) = C(y 0 )∈C (x0) e ik(ŝ +1 q •(ξ 0 -ξ 0 )) R(ŝ +1 q ; ξ 0 ) (54)
for shifting the origin from the center ξ 0 of a level-( +1) cell to the center ξ 0 of a level-cell in the the contribution of a given cell. Then, it is necessary to interpolate, i.e. compute R(ŝ q ; ξ 0 ) at the quadrature points ŝ q from the previously determined values R(ŝ +1 q ; ξ 0 ). The procedure used follows [START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Sylvand | La méthode multipôle rapide en électromagnétisme : performances, parallélisation, applications[END_REF] and is not detailed here.

• A downward pass where the coefficients of local expansions are first computed at level = 2 and then evaluated at lower-level cells by tracing the tree structure downwards. This operation relies on the identity

L(ŝ q ; x 0 ) = e -ik(ŝ +1 q •(x0-x 0 )) L(ŝ q ; x 0 ) (55)
for shifting the origin from the center x 0 of a level-( ) cell to the center x 0 of a level-( +1) cell. This operation is not performed at the root level, i.e. when + 1 = 2.

Then, the contributions of all level-( +1) cells belonging to the interaction list I(x 0 ) of the level-+ 1 cell C(x 0 ) (i.e. all such cells which are nonadjacent to C(x 0 ) but whose father is adjacent to the father of C(x 0 )) are aggregated:

C(y 0 )∈I(x0)
T p (ŝ q , y 0 -x 0 )R(ŝ q ; y 0 ) and the result is added to L(ŝ q ; x 0 ) given by (55) Then, the values L(ŝ q ; x 0 ) are converted to values at the quadrature points ŝ +1 q by a 'reverse interpolation', or 'anterpolation', procedure.

• When the lowest level is reached, all quadratures of the form (52) are finally performed, where x 0 is the centers of a leaf cell, thus evaluating all far-field contributions to the residual at all collocation points. Moreover, for all leaf cells C(x 0 ) and all collocation points x ∈ C(x 0 ), the near-field contributions are computed by evaluating using traditional integration methods the element contributions for all elements located in C(x 0 ) and all cells of same level adjacent to C(x 0 ).

The computation of integral representation formulae for u(x s ), ∇u(x s ) and û(x s ), ∇u(x s ) at all chosen sampling points follows the same approach, with collocation points x ∈ S replaced with sampling points x s ∈ Ω. For sampling points lying in leaf cells not adjacent to any same-level cell intersecting S, the integral representations result from far-field interactions (i.e. fast-multipole contributions) only. Besides, all multipole moments used in this step are those corresponding to the solution of (38) or (39), i.e. those evaluated at the last iteration of the GMRES solution algorithm

Numerical verification of theoretical complexity

To check that the theoretical complexity of the FM-BEM is achieved and illustrate the computational advantage brought by the FM-BEM over the conventional BEM, the simple situation of a spherical acoustic domain subjected on its surface to a uniform normal velocity, is considered. BEM solutions for this problem have been computed for a sequence of meshes with decreasing element size. For each mesh, the prescribed frequency is selected so that the number of nodes per wavelength is (approximately) the same for all meshes. Figure 1 shows the CPU time used for one evaluation of residual (40), using either classical BEM integration techniques (with theoretical complexity O(N 2 )) or the FM-BEM (with theoretical complexity N LogN ), as a function of the number N of nodal unknowns. Functions of the form C 1 N 2 and C 2 N LogN closest to the actual recorded values of CPU(N ) are also plotted on the same graph. The theoretical complexity for both the classical BEM and the FM-BEM are very well verified by the actual CPU times. The FM-BEM, as expected, performs much better for large problems.

Preliminary identification via topological sensitivity: numerical examples

To illustrate the approach described in Sections 3 and 4, the following configuration has been considered: the bounded acoustic domain is the cube defined by each of which consists in applying a uniform excitation over a small region S q (1 ≤ q ≤ 30) of S and recording the acoustic pressure over the whole boundary S (i.e. at all BE mesh nodes), so that S obs = S. The acoustic excitation is such that the wavelength is λ = 3a. Each of the six faces of the cubical domain Ω supports five excitation surfaces S q , each a disk of radius a, arranged as depicted on Fig. 2. The cost function for the inverse problem is defined by

Ω(L) = { |ξ i | ≤ L (i = 1, 2, 3) }, with L = 8a or L = 16a
J (Ω ) = 1 2 30 q=1 S |u -u obs q | 2 dΓ ξ (56)
where u obs q is the data obtained for the q-th applied excitation, with u obs q = u true in the absence of data noise. The centroid x true of true scatterer B true to be identified is located at x true = (2a, 3a, 2a).

To facilitate the graphical interpretation, a thresholded variant T (x s ) of T (x s ) is introduced according to

T (x s ) = T (x s ) , T ≤ C T min 0 , T > C T min with C = 0.25. (57)
The BE meshes used for computing the free field u, the adjoint field û and the simulated data u true q are made of three-noded flat triangular elements, arranged in a regular mesh with approximately 15 nodes per wavelength. For the purposes of computing the simulated error-free data u true for each synthetic experiment, BE meshes of Γ true have been set up as well. Table 1 indicates the numbers of nodes and elements supported by the BE meshes.

First, the identification of one spherical scatterer is considered, for three cases with increasing scatterer radii 0.2a, 0.4a and 0.8a. The field T has been computed for each case on the basis of error-free synthetic data, over a sampling grid S made of 100 × 100 × 100 sampling points located on the vertices of a regular cubic grid, centered at the origin and with grid spacing ∆x s = 16a/101, uniformly filling the whole acoustic domain bounded by S. Figures 3, 4 and 5 show, for each scatterer configuration considered in turn, the distribution of the thresholded topological sensitivity T defined by (57) in the three coordinate planes containing the true scatterer centroid x true . In all cases, T (and hence T ) is seen to attain its lowest values in zones corresponding to, or close to, the actual true scatterer location. It should however be mentioned that T has been observed to also achieve low values in zones located near edges or corners of S, where there is no scatterer. To emphasize and illustrate this remark, Fig. 6 shows the iso-surfaces of the field T corresponding to T = 0.55T min , computed on the entire search grid S (Fig. 6a) and on truncated search grids such that { |ξ i | ≤ 6.5a (i = 1, 2, 3) } (Fig. 6b) and { |ξ i | ≤ 5a (i = 1, 2, 3) } (Fig. 6c), i.e. in which sampling points situated at distances less than 1.5a and 3a, respectively, from S have been taken out of S. Figure 6c shows that, at least in the central region { |ξ i | ≤ 5a (i = 1, 2, 3) } of the acoustic domain, low negative values of T occur only in a small region which is consistent with the actual location of the scatterer.

Then, the effect of data errors is considered for the spherical scatterer of radius 0.8a, by using synthetic data u obs in the form u obs = u true (1 + η), where η are random numbers with zero mean and uniform distribution over the interval [-0.1, 0.1]. Figure 7 shows the distribution of T in the three coordinate planes containing x true . A comparison between Figs. 5 and 7 reveals that the distribution of T is only marginally affected by the data noise.

The examples shown thus far illustrate the capability of T , here defined on the basis of an asymptotic analysis involving vanishing spherical obstacles, to identify the location of spherical obstacles of finite extent. Now, the identification of a non-spherical, box-shaped obstacle whose sides are aligned along the coordinale axes and of finite size 0.8a × 0.8a × 1.6a and whose centroid is still x true = (2a, 3a, 2a), is considered. Figure 8 shows the distribution of T in the three coordinate planes containing x true . The true obstacle is again satisfactorily located.

A last example illustrates the case of a larger acoustic domain Ω(16a), instead of Ω(8a) considered up to this point, with the same wavelength λ = 3a as before. The 'true' scatterer (again a sphere of radius 0.8a) is still located at x true = (2a, 3a, 2a), and hence is located at a larger distance (expressed in wavelengths) from the measurement surfaces. The sampling grid S is now made of 150×150×150 regularly-spaced sampling points, with a grid spacing now of ∆x s = 32a/151. Figure 9 shows the distribution of T in the three coordinate planes containing x true . Presumably as a result of greater remoteness (and hence lower identifiability) of the scatterer, these distributions show, in addition to the correct one, secondary spatial zones where the presence of a small scatterer is also consistent with the data. Still, the lowest values of T furnish a reasonable indication of the true obstacle location, as seen on the 3-D plots of iso-surfaces T = 0.55T min and T = 0.7T min of Fig. 10 for the truncated grid defined by { |ξ i | ≤ 14a (i = 1, 2, 3) }. Again, one notes that values of T close to its minimum T min occur only in the vicinity of the correct obstacle location (excluding, as before, regions close to the external surface S).

Typical CPU times and GMRES iteration counts are provided in Table 2. It is interesting to observe that the overall CPU times for Ω(16a), which involves roughly 4 times as many nodal unknowns as Ω(8a), are about 5 times higher than those for Ω(8a) (while the expected ratio for traditional BEM would be 4 3 = 64), and that the GMRES iteration counts are only fractionally higher for Ω(16a). All computations have been performed on a Linux PC computer with one 3 GHz processor.

Conclusion

In this article, a computationally fast qualitative technique for probing acoustic media for hidden hard obstacles on the basis of overdetermined boundary data, based on the computation via the FM-BEM of the distribution of topological sensitivity of the cost functional associated with the inverse problem, has been presented. Its usefulness has been demonstrated through results of numerical experiments on the qualitative identification of a hard obstacle in a bounded acoustic domain, for configurations featuring O(10 5 ) nodal unknowns and O(10 6 ) sampling points, resulting in overall computing times of a few hours on a 3 GHz PC computer. There is ample scope for increasing these computational sizes. Besides, the proposed approach can be developed for many other physical models, e.g. elastodynamics or electromagnetic waves.

  consists in using the Q(p) = (p+1)(2p+1) quadrature points of the form ŝq = (θ a , φ b ), where cos θ a (1 ≤ a ≤ p+1) are the abscissae for the (p+1)point Gauss-Legendre 1-D quadrature rule over [-1, 1], φ b = 2πb/(2p + 1) (1 ≤ b ≤ 2p + 1) are uniformly-spaced abscissae on [0, 2π]. The associated weights are w q = 2πw θ a /(2p+1), where w θ a are the Gauss-Legendre weights for the (p + 1)-point rule.

2 Fig. 1 .

 21 Fig. 1. CPU time for one evaluation of residual (40) as a function of the number N of nodal unknowns, with matrix-vector products computed using classical BEM techniques (lozenges, theoretical complexity O(N 2 ) or the FM-BEM (circles, theoretical complexity O(N LogN )).

Fig. 2 .

 2 Fig. 2. Pattern of excitation surfaces Sq on each face of external boundary S.

Fig. 3 .x true 1 (a), x2 = x true 2 (b) and x3 = x true 3 (

 3123 Fig. 3. Identification of spherical hard scatterer of radius 0.2a: distribution of thresholded topological sensitivity T (xs) for sampling points xs in coordinate planes ξ1 = x true 1 (a), x2 = x true 2 (b) and x3 = x true 3 (c).

Fig. 4 .x true 1 (a), x2 = x true 2 (b) and x3 = x true 3 (

 4123 Fig. 4. Identification of spherical hard scatterer of radius 0.4a: distribution of thresholded topological sensitivity T (xs) for sampling points xs in coordinate planes ξ1 = x true 1 (a), x2 = x true 2 (b) and x3 = x true 3 (c).

Fig. 5 .x true 1 (a), x2 = x true 2 (b) and x3 = x true 3 (Fig. 6 .

 51236 Fig. 5. Identification of spherical hard scatterer of radius 0.8a: distribution of thresholded topological sensitivity T (xs) for sampling points xs in coordinate planes ξ1 = x true 1 (a), x2 = x true 2 (b) and x3 = x true 3 (c).

Fig. 7 . 1 (a), x2 = x true 2 (b) and x3 = x true 3 (

 7123 Fig. 7. Identification of spherical hard scatterer of radius 0.4a, synthetic data with 10% noise: distribution of thresholded topological sensitivity T (xs) for sampling points xs in coordinate planes ξ1 = x true 1 (a), x2 = x true 2 (b) and x3 = x true 3 (c).

Fig. 8 . 1 (a), x2 = x true 2 (b) and x3 = x true 3 (

 8123 Fig. 8. Identification of box-shaped scatterer: distribution of thresholded topological sensitivity T (xs) for sampling points xs in coordinate planes ξ1 = x true 1 (a), x2 = x true 2

Fig. 9 . 1 (a), x2 = x true 2 (b) and x3 = x true 3 (Fig. 10 .

 912310 Fig. 9. Identification of spherical hard scatterer of radius 0.4a in large domain: distribution of thresholded topological sensitivity T (xs) for sampling points xs in coordinate planes ξ1 = x true 1 (a), x2 = x true 2 (b) and x3 = x true 3 (c).

Table 1 .

 1 Number of element and DOFs supported by the BE meshes.

	Cube size	Cube		Obstacle		Total
		Elements	nodes Elements nodes Elements	nodes
	2L = 16a	76800	38402	336	170	77136	38572
	2L = 32a	307200 153602	336	170	307536 153772

Table 2 .

 2 CPU times and (in parentheses) GMRES iteration count for computing the true, free and adjoint solutions on the boundary, and CPU times for computing the topological sensitivity over the whole sampling grid S.

		u true on S ∪ Γ true	u on S	û on S	T on S
	2L = 16a	1444s (435)	969s (282) 1163s (342)	852s
	2L = 32a	6461s (439) 5615s (388) 6818s (476)	1860s


Florian Lesaint
email: lesaint@isir.upmc.fr

Olivier Sigaud

Jeremy J Clark

Shelly B Flagel

Mehdi Khamassi

Experimental predictions drawn from a computational model of sign-trackers and goal-trackers

Keywords: Reinforcement Learning, Dopamine, Pavlovian conditioning, Autoshaping, Model-based, Model-free, Factored representation, Sign-tracker, Goal-tracker, Conditioned approach 1. Highlights

Gaining a better understanding of the biological mechanisms underlying the individual variation observed in response to rewards and reward cues could help to identify and treat individuals more prone to disorders of impulsive control, such as addiction. Variation in response to reward cues is captured in rats undergoing autoshaping experiments where the appearance of a lever precedes food delivery. Although no response is required for food to be delivered, some rats (goal-trackers) learn to approach and avidly engage the magazine until food delivery, whereas other rats (sign-trackers) come to approach and engage avidly the lever. The impulsive and often maladaptive characteristics of the latter response are reminiscent of addictive behaviour in humans, In a previous article, we developed a computational model accounting for a set of experimental data regarding sign-trackers and goal-trackers. Here we show new simulations of the model to draw experimental predictions that could help further validate or refute the model. In particular, we apply the model to new experimental protocols such as injecting flupentixol locally into the core of the nucleus accumbens rather than systemically, and lesioning of the core of the nucleus accumbens before or after conditioning. In addition, we discuss the possibility of removing the food magazine during the inter-trial interval. The predictions from this revised model will help us better understand the role of different brain regions in the behaviours expressed by sign-trackers and goal-trackers.

• Phasic Dopaminergic activity can be explained by previously engaged features

Introduction

A significant number of models have been developed since the 1970s to describe Pavlovian and instrumental phenomena. Early models were mostly focusing on repro-ducing the averaged behaviour expressed within a population, neglecting inter-individual variations and possibly smoothing the true behaviour of individuals [START_REF] Gallistel | The learning curve: Implications of a quantitative analysis[END_REF], or even masking the variation in behaviour.

However, this variation is of particular interest when trying to identify those individuals within population prone to impulsive behaviours or having a higher risk of addiction (Flagel et al., 2011b;[START_REF] Saunders | Individual variation in resisting temptation: implications for addiction[END_REF]Huys et al., in press).

Recent studies have investigated such intervariability among rats undergoing an autoshaping experiment [START_REF] Flagel | Individual differences in the propensity to approach signals vs goals promote different adaptations in the dopamine system of rats[END_REF][START_REF] Flagel | Individual differences in the attribution of incentive salience to reward-related cues: Implications for addiction[END_REF](Flagel et al., , 2011b,a;,a;[START_REF] Difeliceantonio | Which cue to 'want' ? Opioid stimulation of central amygdala makes goal-trackers show stronger goal-tracking, just as sign-trackers show stronger signtracking[END_REF][START_REF] Mahler | Which cue to "want?" Central amygdala opioid activation enhances and focuses incentive salience on a prepotent reward cue[END_REF][START_REF] Robinson | Dissociating the predictive and incentive motivational properties of reward-related cues through the study of individual differences[END_REF][START_REF] Meyer | Quantifying individual vari-ation in the propensity to attribute incentive salience to reward cues[END_REF][START_REF] Fitzpatrick | Variation in the form of pavlovian conditioned approach behavior among outbred male sprague-dawley rats from different vendors and colonies: Sign-tracking vs. goaltracking[END_REF], where a lever (conditioned stimulus, CS) was presented for 8 seconds, followed immediately by delivery of a food pellet (unconditioned stimulus, US) into an adjacent food magazine. Although no response was required to receive the reward, with training, some rats (sign-trackers; STs) learned to rapidly approach and engage the lever-CS. However, others (goal-trackers; GTs) learned to approach the food magazine upon CS presentation, and made anticipatory head entries into it. Some rats (intermediate group; IG) presented a mixed behaviour, switching between lever and magazine during presentation of the CS, and sometimes engaging both during one trial. Furthermore, in STs, phasic dopamine release in the core of the nucleus accumbens, measured with Fast Scan Cyclic Voltammetry (FSCV), matched the pattern that would be predicted by reward prediction error (RPE) signalling, and dopamine was necessary for the acquisition of a sign-tracking conditioned response (CR). In contrast, despite the fact that GTs acquired a Pavlovian conditioned approach response, this was not accompanied with the expected RPE-like dopamine signal, nor was the acquisition of a goal-tracking CR blocked by administration of a dopamine antagonist (see also [START_REF] Danna | Disruption of conditioned reward association by typical and atypical antipsychotics[END_REF]). While the proportion of STs and GTs in the population varies [START_REF] Fitzpatrick | Variation in the form of pavlovian conditioned approach behavior among outbred male sprague-dawley rats from different vendors and colonies: Sign-tracking vs. goaltracking[END_REF], both phenotypes are typically represented in an outbred population.

To our knowledge, only one model [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF] accounts for these experimental results and has been validated with existing data. This model is built on a combination of model-free and model-based systems [START_REF] Daw | Uncertainty-based competition between prefrontal and dorsolateral striatal systems for behavioral control[END_REF][START_REF] Clark | Pavlovian valuation systems in learning and decision making[END_REF]Huys et al., in press) and extended with state factored representations. Combining multiple systems enables the model to express a large repertoire of behaviours and considering features within states enables the model to learn Pavlovian impetuses [START_REF] Dayan | The misbehavior of value and the discipline of the will[END_REF] specific to the Pavlovian features within the task.

In this paper, we review the model described by [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF], extending it with a new tool to improve its reliability. We suggest new experimental protocols and some new analyses of the data that would further validate the model and strengthen its explanatory power, refine our understanding of the role of the nucleus accumbens in the described behaviours, and help clarify the impact of some choices made in the original protocol.

Material and methods

The model from which the present results are extracted is described in depth in a previous article [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF]. It is composed of two distinct reinforcement learning systems that collaborate to define the action to be selected at each step of the experiment (see Figure 1 A; [START_REF] Clark | Pavlovian valuation systems in learning and decision making[END_REF]).

The first system, a model-based system (MB), incrementally learns a model of the world (a transition function T and a reward function R) from which it infers values (A) for each action in each situation, given the classical following formulas:

Q(s, a) ← R(s, a) + γ s T (s |s, a) max a Q(s , a ) (1) A(s, a) ← Q(s, a) -max a Q(s, a ) (2)
where the discount rate 0 ≤ γ ≤ 1 classically represents the preference for immediate versus distant rewards. At each step, the most valued action is the most rewarding on the long run (e.g. approaching the magazine to be ready to consume the food as soon as its delivery). It favours goal-tracking because this is the shortest path towards the rewarding state (see Figure 1 B).

The second system, a revised model-free system, learns values (V) over features (e.g. food, lever or magazine).

Contrary to the first system, which uses a classical abstract state representation, it relies on the features that compose these abstract states. In traditional reinforcement learning, each situation that can be encountered by the agent is defined as an abstract state (e.g. arbitrarily defined as s 1 , s 2 . . . s x ), such that similarities between situations (e.g. presence of a magazine) are lost. By using features, we reintroduce the capacity to use and benefit from these similarities. The second system is further defined as the feature model-free system (FMF). It relies on a RPE signal δ, computed as follows:

V(c(s, a)) ← V(c(s, a)) + αδ (3) δ ← r + γ max a V(c(s , a )) -V(c(s, a))
where c : S × A → {lever, magazine, f ood, ∅} is a featurefunction that returns the feature c(s, a) the action a was focusing on in state s (e.g. it returns the lever when the action was to engage with the lever). We hypothesized that, similarly to classical model-free systems, δ parallels phasic dopaminergic activity [START_REF] Schultz | Predictive reward signal of dopamine neurons[END_REF]. This signal enables to revise and attribute values, seen as motivational, to features without the need of the internal model of the world used by the MB system. When an event is fully expected, there should be no RPE as its value is fully anticipated.

When an event is positively surprising, there should be a positive RPE. Actions are then valued by the motivational value of the feature they are focusing on (e.g. engaging with the lever would be valued given the general motiva-tional value of the lever). Hence, it favours actions that engage with the most motivational features. This might lead to favour suboptimal actions with regard to maximizing rewards (e.g. engaging with the lever keeps the rat away from the soon to be rewarded magazine). It favours sign-tracking (a suboptimal path, see Figure 1 B) as the lever, being a full predictor of reward, earns a strong motivational value relative to the magazine.

The model does not base its decision on a single system at a time, rather the values of the MB system (A M B ) and the FMF system (V F M F ) are integrated such that a single decision is made at each time step: producing a sort of cooperation between the two systems. The values computed by these two systems are then integrated through a weighted sum and passed to a softmax action selection mechanism that converts them into probabilities of selecting the action given a situation (see Figure 1 A).

The integration is done as follows:

P(s, a) = (1 -ω)A M B (s, a) + ωV F M F (c(s, a)) (4)
where 0 ≤ ω ≤ 1 is a combination parameter which defines the importance of each system in the overall model. Varying ω (while leaving the other parameters of the model unchanged) is sufficient to reproduce the characteristics of the different subgroups of rats [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF]. .

V s ( a 1

) .

Ps(a1)

.

A s ( a 2 )
.

V s ( a 2

) .

Ps(a2)

.

A s ( a 3 )
.

V s ( a 3

) .

Ps

Simulations of experimental protocols

The experiment is described through an episodic Markov Decision Process (MDP) that represents one trial of the session (see Figure 1 B,C). The inter-trial interval (ITI), not being part of the MDP, is simulated between each run by revising downward the magazine value (V F M F (M ) ←

(1u IT I ) × V F M F (M ), u IT I being a parameter of the model). This simulates the hypothesis that the presence of the magazine in the absence of food delivery reduces its value. If the magazine were removed during ITI, we would expect no revision of its value.

The model is used to simulate experiments that involved injections of flupentixol, an antagonist of dopamine, either systemically or within the core of the nucleus accumbens. In the case of local injections, assuming that the FMF system relies on the core of the nucleus accumbens, we simulate the impact of flupentixol on phasic dopamine by degrading the reward predictions errors as follows:

δ ←    δ -f if δ-f δ ≥ 0 0 otherwise (5)
where 0 ≤ f < 1 represents the impact of flupentixol. Its effect is defined such that flupentixol injections cannot lead to negative learning when RPE is positive, but at most blocks it. In the case of systemic injections, we also assume an additional impact on tonic dopamine [START_REF] Humphries | Dopaminergic control of the exploration-exploitation trade-off via the basal ganglia[END_REF], which affects the action selection process. We simulate this impact by revising the temperature parameter (β ← β 1-f ). Hence, flupentixol favours random exploration instead of using learned values to take a decision. Some predictions presented here suggest to lesion the core of the nucleus accumbens. Such a lesion is simulated by removing the FMF system from the model , i.e. all values that would have come from the system are replaced by 0. The rest of the model is left intact. Equation 4 can be replaced by:

P(s, a) = (1 -ω)A M B (s, a) (6)

Index Score

Introduced by [START_REF] Meyer | Quantifying individual vari-ation in the propensity to attribute incentive salience to reward cues[END_REF], the Pavlovian Conditioned Approach (PCA) Index Score provides a metric to categorize rats as STs, GTs or IGs independent of the rest of the population. That is, instead of ordering rats based on their engagement with the lever and splitting the population in 3 groups of approximately equal size, as done in previous studies [START_REF] Flagel | Individual differences in the propensity to approach signals vs goals promote different adaptations in the dopamine system of rats[END_REF][START_REF] Robinson | Dissociating the predictive and incentive motivational properties of reward-related cues through the study of individual differences[END_REF], classifying rats based on PCA Index minimizes the chances of misclassification and allows one to compare across studies or populations of rats. The PCA Index relies on the number of contacts with the lever and the magazine, the probability to engage with one versus the other and the latencies to act towards each (Table 1 in [START_REF] Meyer | Quantifying individual vari-ation in the propensity to attribute incentive salience to reward cues[END_REF]).

We developed a similar Index Score as it provides a good metric for some of the predictions described here.

Simulated rats whose score is > 0.5 are defined as STs.

Simulated rats that have a score < -0.5 are defined as GTs. Remaining rats are defined as IGs. Table 1 explains how it is computed based on the last two sessions of simulations. Contrary to the PCA Index Score, it cannot use latencies as they are not accounted for by the model.

Estimation of model parameters

The model relies on a set of 8 parameters (a shared learning rate, a shared discount rate, a selection temperature, an integration parameter, a ITI impact parameter and 3 initial conditions) that need to be tuned for simulations to fit experimental data. We use the multi-objective algorithm NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: Nsga-ii[END_REF][START_REF] Mouret | SFERESv2: Evolvin' in the Multi-Core World[END_REF] to find the best values (solutions) for the parameters. This method is an efficient tool to fully explore the high dimensional parameter space and avoid local minima.

As in [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF], we search a set of parameter values per group. The two first objectives of the fitness function are to fit the averaged behaviours of the simulated group to the averaged behaviours of the experimental group. More formally, for each group, we try to minimize the least square error between the probabilities of rats and simulated rats to engage with the magazine and the lever over time (see Table 2). This results in multiple solutions that are compromises between these two objectives. We subsequently select one of the solutions that is visually acceptable (no misclassification, and a good compromise between the two other criteria).

We noticed however, that without further constraints, as we are fitting averaged data, some of the resulting solutions could induce great variability of behaviour within a group, leading to misclassification. For example, a simulated rat classified as a GT by its parameters could have behaved as a ST and went undetected as its behaviour would have been diluted in the averaged behaviours of the simulated GT group.

The fitness function was extended with a new criterion based on the Index Score (see Table 2), to favour sets of parameter values that lead to groups of rats that did not introduce such errors, hence without strong inter variability. This is consistent with experimental data [START_REF] Meyer | Quantifying individual vari-ation in the propensity to attribute incentive salience to reward cues[END_REF]. The resulting new sets of parameter values (see Table 3) did not affect the explanatory power of the model.

This metric ensures, for example, that using a set of parameter values for sign-tracking will produce a sign-tracker when applying the model in a simulation reproducing the original experiment. Interestingly, it allows us to predict qualitatively what the behaviour of such a rat (ST in normal conditions) would be in new experimental conditions:

for example, whether the acquisition or the expression of the behaviour would be blocked or shifted to intermediate or even a goal-tracking behaviour, according to the Index Score defined above.

Response Bias(n)

= (LeverP resses-M agazineEntries)/(LeverP resses + M agazineEntries) Probability Difference(n) = p(LeverP ress)-p(M agazineEntry) Type

Score(n) = [ResponseBias(n) + P robabilityDif f erence(n)]/2 Index Score = [Score(6) + Score(7)]/2
ω β α γ u IT I Q i (s 1 , L) Q i (s 1 , ∅) Q i (s 1 , M ) ST 0.

Results

The model has already been validated on a set of behavioural, physiological and pharmacological data [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF]. Interestingly, while the model was only tuned to fit the behavioural data for each group, simulations of additional experiments without changing the parameters were consistent with the remaining experimental data.

The model accounts for the respective engagements of STs and GTs towards distinct specific features [START_REF] Flagel | Individual differences in the propensity to approach signals vs goals promote different adaptations in the dopamine system of rats[END_REF][START_REF] Flagel | Individual differences in the attribution of incentive salience to reward-related cues: Implications for addiction[END_REF](Flagel et al., , 2011b)). It reproduces the difference in patterns of dopaminergic activity for GTs and STs (Flagel et al., 2011b). It also reproduces behaviours indicative of incentive salience attribution, including the conditioned reinforcement effect of the lever shown to a greater extent in STs than GTs [START_REF] Robinson | Dissociating the predictive and incentive motivational properties of reward-related cues through the study of individual differences[END_REF], and the consumption-like engagement of the lever or magazine [START_REF] Mahler | Which cue to "want?" Central amygdala opioid activation enhances and focuses incentive salience on a prepotent reward cue[END_REF][START_REF] Difeliceantonio | Which cue to 'want' ? Opioid stimulation of central amygdala makes goal-trackers show stronger goal-tracking, just as sign-trackers show stronger signtracking[END_REF]. Finally, it also reproduces the impact of flupentixol injected either systemically prior to training (Flagel et al., 2011b), i.e. during acquisition, or locally after the rats have acquired their respective conditioned responses [START_REF] Saunders | The role of dopamine in the accumbens core in the expression of pavlovian-conditioned responses[END_REF], i.e. expression.

In the following sections, taking inspiration from the set of studies used to validate the model, we generate predictions that new experiments or extended analyses of the data could confirm.

Dopaminergic patterns of activity

The model parallels the dopaminergic activity recorded in the core of the nucleus accumbens by Fast Scan Cyclic

Voltammetry with the RPE signal used in the FMF system. At US time, the RPE signal within the FMF system comes from the difference between the value of the previously engaged cue and the value of the delivered food. At CS time, it mainly reflects the value of the most rewarding cue between the lever and the magazine.

STs and GTs dopaminergic patterns at CS and US time are very distinct (Flagel et al., 2011b). While we observe a clear propagation of the signal from US to CS in STs (as expected from the classical RPE theory [START_REF] Schultz | Predictive reward signal of dopamine neurons[END_REF]), this is not the case for GTs for which the CS and US signals are similar to one another and remain relatively constant across sessions (hence, in discrepancy with the classical theory).

In the model, the RPE signal is dependent of the feature previously focused on by the simulated rat. Thus, RPE patterns, averaged over sessions, strongly depend on the dominant path taken by the simulated rats before food delivery. Simulated STs, that mainly engage with the lever before food delivery, have an averaged signal that propagates from US to CS. This reflects that any rat that engages with the lever, eventually learns that it is a full predictor of food delivery. Simulated GTs, that mainly engage with the magazine before food delivery, have an averaged signal that do not show such a propagation. Indeed, the magazine is not fully informative of food delivery for any rat, hence a persistent reward prediction error remains at food delivery when engaging with the magazine during CS.

In Flagel et al. (2011b) The initial analysis (Flagel et al., 2011b) and its reproduction [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF] was done without taking into account the features engaged by animals prior to food delivery, possibly averaging very distinct patterns. The model predicts that if we were to organize the data per groups and actions rather than only per groups, we would observe patterns as shown in Figure 3. At the time the CS is presented, there should be no differences as all rats are exploring the world and not expecting the lever appearance, hence the positive RPE common to all rats. The difference would be at US time.

STs previously engaged with the lever (Figure 3

Removal of magazine during the ITI

In the present model, the simulation of the ITI has a significant impact on the data. We hypothesize that the permanent presence of the magazine during the whole experiment lead animals to revise its associated motivational value, upward at lever retraction (i.e. food delivery) and downward during the ITI as there is no reward to be found then. Hence, on average, its presence does not guarantee access to food. In contrast, the time-locked presence of the lever before food delivery would lead to learn and maintain the motivational value of the lever to a certain level, as its presence guarantees food to be delivered.

First, by keeping the motivational value of the lever higher than that of the magazine in the FMF system, it makes simulated rats favouring this system (STs) to follow a sign-tracking policy. The small contribution of the MB system, which would attract rats towards the magazine does not compensate. Thus, the presence of the magazine in ITI is central for the emergence of STs in the model.

Second, by revising downward the magazine value between episodes, it maintains a discrepancy between the expectation (value) and the observation (reward) at food delivery in simulated rats being engaged with the magazine. This leads to the persistent positive RPE at US time and prevents a full propagation of the signal to CS time.

Thus, the presence of the magazine in ITI is also central for the model, to explain the distinct dopaminergic patterns of activity in STs versus GTs that have been observed in Flagel et al. (2011b).

Third, we also hypothesize that values of the FMF system account for the motivational engagement, i.e. incentive salience, observed in rats towards either the lever or the magazine. The higher motivational value of the lever relative to that of the magazine implies that simulated rats . . chew/bite more the lever than the magazine. While not central to the model, it is consistent with experimental observations [START_REF] Mahler | Which cue to "want?" Central amygdala opioid activation enhances and focuses incentive salience on a prepotent reward cue[END_REF][START_REF] Difeliceantonio | Which cue to 'want' ? Opioid stimulation of central amygdala makes goal-trackers show stronger goal-tracking, just as sign-trackers show stronger signtracking[END_REF].

If no magazine were available during the ITI then, according to the model, the magazine would not loose its motivational value, as it would become a full predictor of food delivery and be highly valued. Hence, we would expect [START_REF] Daw | Model-based influences on humans' choices and striatal prediction errors[END_REF] an increased motivational engagement (chew/bite) towards the magazine, (2) a decreased tendency in signtracking within the population and (3) a different pattern of dopamine activity when goal-tracking for all rats.

As the motivational value of a feature accounts for the level of motivational engagement towards it, a higher motivational value of the magazine, relative to a control group, would necessarily lead to a relatively stronger motivational engagement towards it.

As the motivational value of the magazine would be as high as that of the lever, there should be no reason for rats relying mainly on the FMF system (STs) to favour one over the other, hence shifting to behaviours similar to those of IGs and GTs (see Figure 4). GTs, relying mainly on the MB system would not be deeply affected (see Figure 4 B).

Finally, as the presence of the magazine would be timelocked to the moments before the delivery of food, we would expect a propagation of the dopamine signal from US time to CS time (see Figure 5). At some point (after the value of the food has been fully learned) the signal at US time should start decreasing. Note that if we would have used the same parameters (except for the weighting parameter) to simulate STs and GTs, we would have ex-. . pected an identical RPE signal for STs and GTs, and we know this is not the case based on existing data (Flagel et al., 2011b).

Sign-Trackers

Intermediate

The expected decreased tendency in sign-tracking within the simulated population does not mean that simulated rats would not be attracted any more by the lever. Simulated rats would indeed be attracted by both the lever and magazine because their FMF system attributes a high motivational value to all signs preceding reward delivery.

Combined with the contribution of the MB system which attracts rats towards to magazine, it could make the simulated animal engage more with the magazine than with the lever. Thus if the computational model is valid, this would mean that the tendency to sign-track in real animals can be gradually changed by affecting some of the signs or features present in the context of the task (here the magazine during the ITI).

Injections of flupentixol in the core of the nucleus accumbens

In the model, flupentixol, an antagonist of dopamine, is hypothesized to impact the RPE (hypothesized to paral-lel phasic dopamine) used in the FMF system, putatively based within the core of the nucleus accumbens. Flupentixol is also assumed to affect any action selection process, relying on tonic dopamine [START_REF] Humphries | Dopaminergic control of the exploration-exploitation trade-off via the basal ganglia[END_REF].

Hence, under systemic injections of flupentixol, the learning process of the FMF system is disrupted and actions are almost randomly picked barely using learned values.

With systemic injections of flupentixol (Flagel et al., 2011b), no goal-tracking nor sign-tracking is expressed in the population. However, when afterwards released from flupentixol, GTs fully express goal-tracking, whereas STs behave as untrained rats.

The model accounts for the absence of behaviours under flupentixol by the hypothesized impact of flupentixol on the action selection process, blocking the expression of any acquired behaviour [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF]. The subsequent absence of sign-tracking on a last session free of flupentixol is explained by the disruption of the FMF system during the 7 first sessions, blocking behaviour acquisition.

The full expression of goal-tracking as soon as flupentixol is removed, relies on the unaffected learning process in . . signal at CS and US seems to follow similar trends and there is no propagation of signal from US to CS. (B) When magazine is time-locked to CS presentation, the value of US is propagated to the CS. Thus, the signal at US time, after sufficient learning (2 first sessions) start decreasing in favour of the CS time.

the MB system, assumed to be dopamine-independent and hence keeps learning under flupentixol, but which values are simply not used by the softmax function.

The model predicts that if flupentixol were injected locally in the core of the nucleus accumbens rather than systemically prior to acquisition, GTs would normally express their behaviour, as the action selection mechanism would not be disrupted and make use of the values learned in the MB system; whereas STs' behaviour would remain blocked because of the disruption of the FMF system (see (Flagel et al., 2011b) on the probability to approach the lever for STs (A) and the magazine for GTs (B) during lever presentation. Last session is without flupentixol. Under flupentixol (first 7 sessions), both sign-tracking and goal-tracking are blocked. On the last flupentixol-free session (8 session), STs are unable to express sign-tracking, its learning having been blocked, whereas GTs fully express goal-tracking, which learning was only covert. (C,D) Simulation of the impact of local injections of flupentixol in the core of the nucleus accumbens, hypothesised to impact only the FMF-system. Contrary to the initial experiment, the injections being localized to the FMF-system, the action selection mechanism is not impacted. Hence, GTs fully express goal-tracking during the first 7 sessions (C). STs are still unable to express sign-tracking (D).

mechanism remains functional and makes use of the MB system values, such that the behaviour of GTs is preserved while the one of STs is disturbed and leads to a decrease in sign-tracking in the overall population.

We expect that lesions of the core of the nucleus accumbens would lead to similar effects as the above experiments.

Lesions of the core of the nucleus accumbens prior to the experiment would (1) block the expression of signtracking responses and (2) stop the motivational engagement towards the magazine or the lever during approaches.

By disabling the FMF system (setting and keeping all values to 0), it cannot favour the lever over the magazine any more. STs would therefore act randomly , approaching lever and magazine indifferently, as observed in IGs. We would expect a shift towards goal-tracking similar to the one expected for removing the magazine during the ITI (as in Figure 4). However, while a magazine removal would lead to an increase in motivational engagement, we expect such a lesion to block any consumption-like behaviour. Especially, we would expect GTs' approach behaviour to remain similar to control group, but without subsequent chewing and biting of the magazine.

We would expect that lesions of the core of the nucleus accumbens after the experiment would disrupt the ten-dency to sign-track in the overall population, while leaving the tendency to goal-track intact (see Figure 7). However, contrary to flupentixol injections, that needed 35 min of infusion for a visible effect, we would expect the effect to be immediate with a lesion. Such a lesion would disrupt the FMF system, hence [START_REF] Daw | Model-based influences on humans' choices and striatal prediction errors[END_REF] suppressing any consumptionlike engagement towards the features (motivational values being kept to 0), and (2) stop favouring engagements towards the lever. The lesion would leave the MB system unaffected and have no impact on the general tendency to goal-track.

Discussion

Relying on a model that was previously validated using experimental data to account for variability in rats undergoing an autoshaping paradigm [START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF], we generate an additional set of behavioural, physiological and pharmacological predictions.

We predict that dopaminergic patterns for IGs should be a mixed signal between those observed for STs and GTs.

We predict that looking separately at the DA patterns given the prior engagement towards either the lever or the magazine should lead to clearly distinct patterns. We predict that the removal of the magazine during the ITI should lead to an increased motivational engagement towards the magazine, a decreased tendency in sign-tracking within the population and a different pattern of dopaminergic activity when goal-tracking. Finally, we predict that local injections of flupentixol to the core of the nucleus accumbens would preserve goal-tracking and prevent the learning of a sign-tracking response, a result that should also be observed following lesions of the core of the nucleus accumbens prior to conditioning. Lesions after conditioning, would only block the expression of the learned sign-tracking behaviour.

An important limitation of the present predictions is that most of them are based on the behaviour that is expected to emerge from naive rats trained in a revised proto-col, assuming that they would have behaved in a specific manner in the standard protocol (e.g. expecting a supposed ST to goal-track). To overcome this difficulty, one must look at the population level rather than the individual level [START_REF] Saunders | The role of dopamine in the accumbens core in the expression of pavlovian-conditioned responses[END_REF], which might be problematic as the proportion of GTs, STs and IGs is highly variable in a population [START_REF] Meyer | Quantifying individual vari-ation in the propensity to attribute incentive salience to reward cues[END_REF][START_REF] Fitzpatrick | Variation in the form of pavlovian conditioned approach behavior among outbred male sprague-dawley rats from different vendors and colonies: Sign-tracking vs. goaltracking[END_REF]. An alternative would be to use selectively bred rats that can more or less be ensured to behave as STs or GTs in experimental conditions (Flagel et al., 2011b).

Another limit of the present predictions are the hypotheses on which they are based. It cannot be excluded that the core of the nucleus accumbens also contributes to the MB system, but not by its dopaminergic activity [START_REF] Khamassi | Integrating cortico-limbicbasal ganglia architectures for learning model-based and modelfree navigation strategies[END_REF][START_REF] Van Der Meer | Ventral striatum: a critical look at models of learning and evaluation[END_REF][START_REF] Mcdannald | Ventral striatum and orbitofrontal cortex are both required for model-based, but not model-free, reinforcement learning[END_REF] (but see van der One could argue that, to some extent, describing STs with a MF system and GTs with a MB system could be sufficient to explain dominant behaviours [START_REF] Clark | Pavlovian valuation systems in learning and decision making[END_REF]. However, it would fail to explain the full and continuous spectrum of observed behaviours [START_REF] Meyer | Quantifying individual vari-ation in the propensity to attribute incentive salience to reward cues[END_REF]. If the predictions that we make about IGs (which have an intermediate behaviour between STs and GTs) . . are correct, this would argue in favour of a continuum in the weighting between MB and MF systems rather than a pure dichotomy.

Control Lesion

An alternative to the collaboration of both systems (through a weighted sum) would be a reciprocal inhibition, such that only one system would be working at a time. This would be sufficient to account for the previous point and may even be able to account for the absence of RPE pattern in the dopamine signal measures in GTs (Flagel et al., 2011b) without requiring a revision of the magazine value during ITI. The inhibition of the MF system in GTs would indeed prevent any RPE signal from being observed. However, it would be unable to properly account for the consumption-like engagement observed in both STs and GTs without some kind of extension (see [START_REF] Zhang | A neural computational model of incentive salience[END_REF] for a computational model of incentive salience). It would also fail to explain why the pharmacological disruptions of one system does not seem to let the other take control (Flagel et al., 2011b).

Another possibility would be that the two systems run in parallel but that only one is used to make the decision during a trial. Assuming that one system leads to the lever and the other to the magazine, we would ex-pect IGs to behave as STs when engaging with the lever and GTs when engaging with the magazine. Experimental data goes against such interpretation. Meyer et al.

(2012) observed that contrary to STs or GTs, IGs tend to approach both the magazine and the lever during single trials. Some rats even hold on to the lever while putting their head into the magazine (which no model that selects a single action at a time can reproduce). While the task representation does not allow multiple engagements in a trial, this suggests that both systems are active and contribute actively to their behaviour at all time. We would also expect rats to behave differently when using one system over the other, such that, for example, rats would actively engage with the lever but quietly wait in front of the magazine, which is not the case. Finally, the recent literature seems consistent with multiple systems working in parallel and partially contributing to a global decision (e.g. [START_REF] Daw | Model-based influences on humans' choices and striatal prediction errors[END_REF]). Hence, this does not suggest take-over competition between the systems. Trial-by-trial analyses [START_REF] Daw | Trial-by-trial data analysis using computational models[END_REF] would allow us to definitely rule out such alternatives. Finally, if only the output of the MF system was inhibited, given that the lever appearance is fully predictive of food delivery, no classical MF system This suboptimal behaviour is, however, in accordance with behavioural observations in rats.

Subsequent studies could benefit from a different approach to estimate parameters. We are currently fitting the model on the behavioural data per sessions and groups, using trial-by-trial analyses could prove a better tool to fit the parameters at the individual level [START_REF] Daw | Model-based influences on humans' choices and striatal prediction errors[END_REF] and comfort some choices in the architecture of the model.

It has been suggested that individuals for whom cues become powerful incentives (i.e. STs) are more prone to develop addiction [START_REF] Saunders | The role of dopamine in the accumbens core in the expression of pavlovian-conditioned responses[END_REF]. Thus, the current model and its predictions will allow us to further investigate and possibly identify the neural mechanisms that underlie addiction and related disorders. For example, the current model predicts that some manipulations could alter the behaviour of STs towards that of GTs, and the neurobiological targets of these manipula-

  , recordings of dopaminergic activity in outbred rats were made to parallel those of the selectively bred STs and GTs but no recordings were made in outbred IGs. We would expect that IGs, whose behaviour fluctuate between sign-tracking and goal-tracking, would have a kind of mixed signal, averaging between those following from sign-tracking and goal-tracking. The current parameters values used in the model suggest that we would expect a high signal at CS time that would converge to a certain point, while at the meantime, the signal at US time would keep fluctuating without fully disappearing (see Figure 2).Note that the visual results of this prediction are not identical with those in[START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF]. Contrary to ST and GT behaviours that deeply rely on the mechanisms, IG results strongly depend on the parameter values, which are significantly different with the introduction of the new score. Experimental recordings could help us refine the appropriate set of values for further predictions.

Figure 2 :

 2 Figure 2: Prediction of the model about expected patterns of dopaminergic activity in intermediate rats. Data are expressed as mean ± SEM. Average RPE computed by the FMF system in response to CS and US presentation for each session of conditioning in the intermediate group.

Figure 3 :

 3 Figure 3: Predictions about patterns of dopaminergic activity per groups and per actions. Average RPE computed by the FMF system in response to CS and US presentation for each session of conditioning for STs (A,C) and GTs (D,B) when engaged with the lever (A,B) or with the magazine (C,D). The model predicts very distinct patterns of activity depending on the feature engaged with prior to food delivery.

Figure 4 :

 4 Figure4: Distribution of rats given the removal of the magazine during the ITI. Simulated rats ordered by their Index Score. In blue simulated rats using parameter values tuned for GTs, in red for STs and in white for IGs in the classical condition (control). Rats with a score < -0.5 are GTs, with a score > 0.5 are STs and remaining rats are IGs. (A) As expected, rats using parameters' value for GTs are classified as GTs. Same for STs and IGs. (B) Without magazine during the ITI, simulated rats that would have been classified as GTs in normal conditions are still classified GTs. However, rats that would have been classified as STs (red) have a score that classify them as GTs or IGs. One IG (white) is now classified as GT.

Figure 5 :

 5 Figure 5: Patterns of dopaminergic activity for GTs given the removal of the magazine during ITI. Average RPE computed by the FMF system in response to CS (black) and US (blue) presentation for GTs for each session of conditioning. It is hypothesized to parallel the patterns of dopaminergic activity observed by FSCV in the core of the nucleus accumbens. (A) With the classical protocol (control),

Figure 6 )Figure 6 :

 66 Figure 6), and this is indeed what happend when Saunders and Robinson (2012) locally injected flupetixol after the behaviours were already acquired.

  [START_REF] Van Der Meer | Triple dissociation of information processing in dorsal striatum, ventral striatum, and hippocampus on a learned spatial decision task[END_REF];[START_REF] Bornstein | Dissociating hippocampal and striatal contributions to sequential prediction learning[END_REF];[START_REF] Penner | Neural systems analysis of decision making during goal-directed navigation[END_REF]). Hence, completely disrupting it might unexpectedly affect goal-tracking. Validating these predictions would help to confirm this hypothesis. In the initial model[START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF] we interpreted the parameter which simulates the ITI as accounting for the engagement of the rats towards the magazine during the ITI. Preliminary analyses of experimental data (not shown), while still inconclusive, tend to mitigate such a strong hypothesis. Hence, in the current article, we only assume that the presence of the magazine during the ITI impacts its general motivational value within the experiment. Validating such predictions would definitely help to clarify the impact of the ITI context on the expressed behaviours.

Figure 7 :

 7 Figure7: Predictions of the impact of lesions of the core of the nucleus accumbens after conditioning. General tendencies to sign-track (A) and goal-track (B) in a population of rats after training. Lesion of the core of the nucleus accumbens is simulated by a blockade of the FMF system in the model. We expect a decrease in the tendency to sign-track (A) with a lesion (purple) relative to a control group (black). General goal-tracking tendencies should remain unchanged (B).

(

  relying on classical state representation) would reproduce the differences observed in phasic dopaminergic patterns between STs and GTs nor explain the differences of focused features. Hence, the model suggests to take features into consideration. The interest of the current computational model lies in its combination of simple concepts actively used and accepted in the current field (Dual reinforcement learning and factored representations) but rarely used together, to account for a variability of experimental data, without resorting to arbitrary additions. As a result the current model does not behave as state of the art algorithms would on the same task and produces a suboptimal behaviour.

  tions may be used to alter drug-cue dependency and prevent relapse (For further discussion regarding the role of learning-related dopamine signals in addiction vulnerability, seeHuys et al. (in press)).To conclude, the current article refines the model previously described by[START_REF] Lesaint | Modelling individual differences in the form of pavlovian conditioned approach responses: A dual learning systems approach with factored representations[END_REF] with an additional metric that strengthens its explanatory power. It mainly suggests a set of predictions with which to further confront the model. The new proposed experiments would help to better localize the anatomical counterparts of the mechanisms involved and disentangle their contributions to the observed behaviours. It would also help in refining the hypotheses and simplifications of the model and we hope would confirm the interest and necessity of considering the features rather than the general situations encountered by rats when modelling this kind of phenomena.

  Free system (FMF, in red) which provide respectively an advantage function A and a value function V values for actions a i given a state s. These values are integrated in P, prior to be used into an action selection mechanism. The various elements may rely on parameters (in purple). The impact of flupentixol on dopamine is represented by a parameter f that influences the action selection mechanism and/or any reward prediction error that might be computed in the model. (B) MDP accounting for the experiments described in[START_REF] Flagel | Individual differences in the attribution of incentive salience to reward-related cues: Implications for addiction[END_REF]Flagel et al. ( , 2011b));[START_REF] Robinson | Dissociating the predictive and incentive motivational properties of reward-related cues through the study of individual differences[END_REF];[START_REF] Meyer | Quantifying individual vari-ation in the propensity to attribute incentive salience to reward cues[END_REF]. States are described by a set of variables: L/F -Lever/Food is available, cM/cL -close to the Magazine/Lever, La -Lever appearance. The initial state is double circled, the dashed state is terminal and ends the current episode. Actions are engage with the proximal stimuli, explore, or go to the M agazine/Lever and eat. The path that STs should favour is in red. The path that GTs should favour is in dashed blue. (C) Time line corresponding to the unfolding of the MDP.

	A	.	B	close to Lever .	.	g o L	s 2	. . . . cL	eng	s 5	. . cL . .	.
				s 0	. . . . .	exp	s 1	. . . . . La	exp	s 3	. . . . . . . . . . . . L . . . . . . . . . . . . . . . . . . . . .	∅	s 6	. . . . F .	goM s 0	. . . . .	ITI
	. .	(a3)	C	close to Magazine . .	g o M	s 4	. . . . L . cM	eng	. . . . cM F . goM s 7	e a t ( R )	.
				Lever (CS)
		a i		Food (US)
	Figure 1: Model and Markov Decision Process used for simulations. (A) The model is composed of a model-based system (MB, in
	blue) and a Feature-Model-

L . F .

Table 1 :

 1 Formulas for deriving the Index Score. The Index Score provides a way to classify rats as STs, GTs or IGs, independently of the rest of the population. It relies on averaging scores computed for the last two sessions of the simulations. The Score for session n is derived by averaging its Response Bias and its Probability Difference. Responses Bias is a ratio between the difference in lever presses versus magazine entries and the total number of entries. Probability Difference is the difference between the probability to engage with the lever and the probability to engage with the magazine.

	Objective

i ∈animals |ref P CA(Group) -IndexScore ai (Group)|)

Table 2 :

 2 Revised fitness function. Lists of the multiple objective/criterion of the fitness function applied to each simulated group. refPCA is 1, 0 and -1 for STs, IGs, GTs respectively. This function is combined with NSGA-II to retrieve parameter values that best reproduce the experimental results. It results in a Pareto front of parameters from which we select by hand the solution that is consistent (no agent being misclassified) and that best visually fits the observed data (between engaging with the lever versus engaging with the magazine).

Table 3 :

 3 Parameters

used to produce the presented results. All results were generated based on the same parameters. Some parameters might not be used or erased depending on the specific experimental protocol simulated.

Acknowledgments