P Moireau

C Bertoglio

N Xiao

C A Figueroa

C A Taylor

D Chapelle

J-F Gerbeau

Visco-elastic support has been previously established as a valuable modeling ingredient to represent the effect of surrounding tissues and organs in a fluidstructure vascular model. In this paper, we propose a complete methodological chain for the identification of the corresponding boundary support parameters, using patient image data. We consider distance maps of model to image contours as the discrepancy driving the data assimilation approach, which then relies on a combination of (1) state estimation based on the socalled SDF filtering method, designed within the realm of Luenberger observers and well-adapted to handling measurements provided by image sequences, and (2) parameter estimation based on a reduced-order UKF filtering method which has no need for tangent operator computations and features natural parallelism to a high degree. Implementation issues are discussed, and we show that the resulting computational effectiveness of the complete estimation chain is comparable to that of a direct simulation. Furthermore, we demonstrate the use of this framework in a realistic application case involving hemodynamics in the thoracic aorta. The estimation of the boundary support parameters proves successful, in particular in that direct modeling simulations based on the estimated parameters are more accurate than with a previous manual expert calibration. This paves the way for complete patient-specific

Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data 1 Introduction

In the past decade, numerous works have broached the topic of fluid-structure interaction (FSI) in arteries, with an increasing focus on the concept of patient-specific modeling. This polysemic term was initially mainly associated with the anatomical personalization of the modeled arteries, which justified major efforts to improve segmentation techniques in order to produce patientspecific geometries [START_REF] Cebral | Efficient pipeline for image-based patient-specific analysis of cerebral aneurysm hemodynamics: technique and sensitivity[END_REF][START_REF] Gerbeau | Fluid-structure interaction in blood flows on geometries based on medical imaging[END_REF][START_REF] Alastrué | Numerical framework for patientspecific computational modelling of vascular tissue[END_REF][START_REF] Xiong | Simulation of blood flow in deformable vessels using subject-specific geometry and assigned variable mechanical wall properties[END_REF]. In addition, many modeling specific components have been incorporated into the FSI formulation to improve the representation of the patient flow characteristics. In particular, a significant emphasis has been placed on fluid boundary conditions -see for example (Figueroa et al 2006;[START_REF] Formaggia | On the coupling of 3D and 1D Navier-Stokes equations for flow problems in compliant vessels[END_REF][START_REF] Kerckhoffs | Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation[END_REF][START_REF] Quarteroni | Coupling between lumped and distributed models for blood flow problems[END_REF] -to model the impact of the part of the arterial tree truncated from the region of interest. Given the complexity of the truncated arterial tree, it is common to substitute it with simplified lumped parameter models at the inlet and outlets, specifically calibrated with the available pressure and flow data available for the patient (see for example [START_REF] Brown | Importance of realistic LVAD profiles for assisted aortic simulations: evaluation of optimal outflow anastomosis locations Buerger C, Schaeffter T, King A[END_REF]Troianowski et al 2011)).

With recent advances in medical imaging technology such as cardiac-gated Computed Tomography (CT) and 4D Magnetic Resonance Imaging (MRI) that can now provide time-resolved data on the motion of the arterial wall over the cardiac cycle, a new wealth of in-formation is available to inform the simulation efforts still further. An approach -sometimes called "imagebased CFD", see e.g. [START_REF] Moreno | Nonlinear-transformation-field to build moving meshes for patient specific blood flow simulations[END_REF][START_REF] Piccinelli | 4D image-based CFD simulation of a compliant blood vessel[END_REF] -has been proposed whereby the flow is computed in a domain whose deformation is directly deduced from the medical images, thus avoiding the resolution of the coupled fluid-solid problem. However, this type of approach features important limitations in that inaccuracies inherent to such measurements then directly produce significant error sources in the modeling simulations, and some important physical phenomena intrinsic to FSI -such as pressure waves -are not adequately captured. By contrast, a different approach has been adopted in some cardiac and arterial FSI applications [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF]Bertoglio et al 2012b), using the time-resolved image data extracted from images and a cardiac mechanical model in [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF], and synthetic data with a fully coupled FSI model in (Bertoglio et al 2012b), to estimate some patient-specific characteristics -hence, enabling an additional degree of model personalization -based on a data assimilation strategy.

In this paper, we further demonstrate this strategy applied with FSI in arteries, with a particular focus on taking advantage of the available image data to estimate modeling parameters pertaining to vessel wall boundary conditions. In (Moireau et al 2011), we presented a framework to characterize the impact of external tissues and organs -in particular the spine for the aorta -on the arterial wall, by a simple viscoelastic term prescribed on the surface boundaries of the arterial wall. When applied to an inlet or outlet wall boundary, this viscoelastic term represents the truncated arterial tree extending from the geometric model. Comparably, when applied along the exterior surface of the arterial wall, the viscoelastic term models tethering to the external organs. This approach was shown to dramatically improve the quality of the simulations, providing accurate deformation patterns of the arterial tree. Of course, this requires a careful calibration of the parameters in this boundary support model. To that effect, time-resolved medical image data provides information on the apparent motion of the walls, which is valuable to determine how strong the external tissue support should be.

Whereas the calibration of the boundary support parameters was performed manually in (Moireau et al 2011), in the present article we propose a method to automate and refine this calibration process using the segmentations extracted from the image sequence, without compromising the simulation time which can already be quite long for this type of complex simulation. To accomplish this, we rely on an original estimation frame-work developed in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF][START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]Moireau and Chapelle 2011b) which proposes a strategy for using the available additional time-resolved data to identify various uncertain modeling quantities, while improving the quality of the numerical computations. The possibility to benefit from these additional data in the computation can be contemplated in the framework of data assimilation [START_REF] Blum | Data assimilation for geophysical fluids[END_REF][START_REF] Blum | Data assimilation for geophysical fluids[END_REF], which consists in merging model and data information in order to circumvent the initial model uncertainties and to improve the numerical prediction. In this context, two types of strategies are often considered:

-A variational approach, which estimates the uncertainties by minimizing a least square criterion involving the discrepancy between the data and the corresponding outputs from a model simulationsee [START_REF] Blum | Data assimilation for geophysical fluids[END_REF][START_REF] Chavent | Nonlinear Least Squares for Inverse Problems[END_REF]) and references therein, and [START_REF] Perego | A variational approach for estimating the compliance of the cardiovascular tissue: An inverse fluid-structure interaction problem[END_REF] in the context of hemodynamics; -A sequential approach, which filters the uncertainties over time to stabilize the computed numerical system onto the actual partially-observed systemsee [START_REF] Simon | Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF]) and references therein.

These two approaches are very different in their practical use, but rely on the same fundamental observability condition, which determines that observing the system even partially -through a time sequence of boundary shapes in our study -is sufficient for circumventing the lack of initial knowledge of the system [START_REF] Simon | Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF]. Furthermore, in various cases it can be proven that the two strategies are actually equivalent [START_REF] Kailath | Linear estimation[END_REF]. With the objective of a fast model identification "on the fly", we adopt herein a sequential approach whereby the uncertainties are evaluated during the simulations by considering the various available observations sequentially. In a very general framework, we consider two classes of uncertainties:

-Parameter uncertainties. These uncertainties are considered by assuming that some of the model constants, in some regions, remain imperfectly assessed after calibration. In particular, in this article we consider that the external tissue support boundary conditions -which constitute a rather simplified model of the complex attachment of the aorta on the spine -are difficult to calibrate manually. Since we define only a few regions with homogeneous coefficients, the total number of parameters to be identified is finite and reasonably small. Note that considering parameters with the same discretization level as the PDE field is possible in theory, but would certainly lead to a very sensitive or even ill-posed identification problem. Moreover, the small num-ber of parameters to estimate should allow the use of classical optimal methods in a sequential framework inspired from Kalman filtering, as presented in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF]Moireau and Chapelle 2011b) and applied to the identification of elastic properties of arterial walls in the fluid-structure interaction context in (Bertoglio et al 2012b). -State uncertainties. Other sources of uncertaintiessometimes neglected -can be categorized in what we call the state uncertainties. This concerns the lack of knowledge on the initial conditions, the approximation errors coming from the discretization procedure, and even any other uncharacterized model error and approximation. As discussed in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF], for a mechanical model these uncertainties always exist in practice and affect the solution fields of the PDEs. Therefore, they should be adequately controlled, i.e. filtered jointly with the identification of the parameters through the data assimilation procedure to avoid an unstable estimation behavior. Furthermore, as shown in [START_REF] Chapelle | Improving convergence in numerical analysis using observers. the wave-like equation case[END_REF] in the context of second-order hyperbolic equations, even with perfectly known parameters, filtering the state uncertainties with adequate -reasonably noisy -data always improves the numerical approximation of the system.

Concerning sequential approaches, the main difficulty in taking into account these uncertainties lies in the fact that most of these methods are based on optimization principles which make them not easily suited to PDE systems due to the "curse of dimensionality" coined by R.E. Bellman [START_REF] Bellman | Dynamic Programming[END_REF]. This is the reason why we adopt here instead an observer strategy which consists in adding to the original model a consistent stabilization term based on the measured discrepancy between the computed outputs of the model at a given time and the observation at hand. When the numerical system matches the actual observations, the stabilization term vanishes. Otherwise, the dynamics of the model are corrected to track the actual trajectory. This idea was initially introduced in the control community [START_REF] Luenberger | An introduction to observers[END_REF], but was only recently applied to PDEs as the so-called Luenberger estimation in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF][START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF] or as, in a more recent terminology, the nudging assimilation method in [START_REF] Auroux | A nudging-based data assimilation method: the Back and Forth Nudging (BFN) algorithm[END_REF]. In this article we present an observer strategy based on [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF], which can be used with the type of measurements available -i.e., the segmented surfaces of the fluid-structure interface -and has already been studied in its linearized form in the context of a solid coupled with a fluid in (Bertoglio et al 2012a).

The outline of this paper is as follows. First, we recall the model considered for the case of interest and describe the available data. Here, we focus on one of the cases studied in (Moireau et al 2011), a subject in whom the aortic wall motion is rather large, hence the fluid-structure interaction model involves an ALE formulation for the fluid with a non-linear large displacement hyperelastic solid model for the arterial wall. Then, we present our sequential estimator and carefully describe the handling of the specific measurements available here, namely, a time sequence of segmented surfaces. Finally we report on the results and discuss the accuracy, robustness and computational effectiveness of our method.

Models

We consider here a fluid-structure interaction (FSI) model of blood flow compatible with large displacements of the vessel walls. The model -described in detail in (Moireau et al 2011) -consists of incompressible Navier-Stokes equations coupled with a visco-hyperelastic formulation. The other salient modeling features concern the boundary condition specification, namely, -For the fluid, the downstream boundary conditions are modeled using a "coupled multi-domain" approach in which a so-called Windkessel lumped parameter model representing the downstream vasculature provides a dynamic relationship between outlet pressure and flow for the FSI solution (Vignon-Clementel et al 2010); -For the solid, a visco-elastic support model is used to represent the interaction of the vessel outer boundaries with the surrounding tissues, as well as the effect of the truncated arterial trees. A forcing term was used at the inlet boundary of the vessel wall to account for the heart motion; these simple boundary conditions were found to be very effective in producing simulations that compare well with the available clinical data (Moireau et al 2011).

Lastly, another difficulty arising from patient-specific FSI modeling is that the segmented medical image data correspond to deformed configurations of the vessel walls. Hence, a rather complex inverse problem must be solved in order to determine a corresponding stress-free configuration [START_REF] Gee | A computational strategy for prestressing patient-specific biomechanical problems under finite deformation[END_REF]Moireau et al 2011), which in turn implies that the initial condition in the solid -also taken as the reference configuration -is pre-stressed.

Model formulation

We now introduce some notation in order to summarize the model equations. Starting with geometrical and kinematical considerations, we denote by Ω f (t) the fluid domain at any time t, with u f the fluid velocity, whereas the solid domain is denoted by Ω s (t) with u s the solid velocity, and

Σ(t) = Ω f (t) ∩ Ω s (t) is the fluid-solid in- terface, see Figure 1. The deformed solid domain Ω s (t)
is the image of a reference configuration Ω s 0 through the solid deformation mapping given by

φ s : Ω s 0 × [0, T ] -→ Ω s (t), (ξ, t) -→ x = ξ + y s (ξ, t)
where y s denotes the solid displacement, hence u s = ∂ t y s . With a view to an Arbitrary Lagrangian-Eulerian fluid formulation, we also define a fixed fluid domain

Ω f 0 such that Ω f 0 ∩ Ω s 0 = Σ 0 = φ -1 s Σ(t) ,
φ A : Ω f 0 × [0, T ] -→ Ω f (t), (ξ, t) -→ x = ξ + Ext f y s | Σ0
which is compatible with the solid mapping in that φ A | Σ0 = φ s | Σ0 . Note that it would be sufficient to ensure the compatibility of the normal components, but it is simpler, and adequate for the present work, to ensure the continuity of the three components. We introduce the corresponding Jacobian

J f (ξ, t) = det ∇ ξ φ A (ξ, t)
and the fluid domain velocity

w = ∂ t φ A = Ext u s | Σ0 .
Regarding the remaining physical notation, for the fluid we denote by ρ f the constant mass per unit volume, p f the pressure, σ f the stress tensor given by

σ f = -p f 1 + 2µ ε(u f ),
with µ the dynamic viscosity and ε the symmetrized gradient operator. As for the solid we define ρ s as the (constant) mass per unit volume in the reference configuration, σ s the Cauchy stress tensor and Π s the associated first Piola-Kirchhoff stress tensor, both defined based on a hyperelastic potential W s , with additional viscous effects represented by a viscous modulus η s , see (Moireau et al 2011) for details. Finally, the strong form of the FSI problem can be written as follows

               w = Ext f (u s | Σ0 ), in Ω f 0 ρ f J f ∂J f u f ∂t ξ + ρ f (u f -w) • ∇ x u f -2µdiv x (ε(u f )) + ∇ x p f = 0, in Ω f (t) div x u f = 0, in Ω f (t) (1a)                ρ s ∂u s ∂t -div ξ (Π s ) = 0, in Ω s 0 ∂y s ∂t = u s , in Ω s 0 F -1 Π s = ∂W s ∂e + η s ė, in Ω s 0 (1b) u f = u s , on Σ(t) Π s • n s = J f σ f (u f , p f ) • (F f ) -T • n s , on Σ 0 (1c)
where ∂ ∂t ξ

represents the so-called ALE time derivative, see [START_REF] Fernández | A coupled momentum method for modeling blood flow in three-dimensional deformable arteries[END_REF] and references therein, and n s denotes the outward unit normal vector at the boundary of the reference solid domain Ω s 0 . Note that these modeling equations have been gathered into 3 groups corresponding to the fluid, solid and coupling equations, sequentially. Concerning the solid boundary conditions, we distinguish between Γ s in (see Figure 1), the vessel wall boundary at the inlet face of the model where a time-dependent forcing term y b representing the cardiac motion is precribed, and Γ s , the outer vessel wall and outlet face vessel wall boundaries where the viscoelatic tissue support boundary condition is prescribed:

σ s • n s = -k s (y s -y b ), on Γ s in , (2)
σ s • n s = -k s y s -c s u s , on Γ s . (3)
We do not include a viscous effect on Γ s in , primarily to avoid differentiating the displacement y b -extracted from the image sequence, hence containing substantial measurement errors -to compute the corresponding forcing velocity. In practice, we adopt a rather large value for the stiffness coefficient k s , and thus the elastic boundary condition can be seen as a penalized Dirichlet condition.

As for the boundary conditions for the fluid, the boundary of the fluid domain Ω f is decomposed into the inlet (connected to the aortic root) Γ f in , the outlets Γ f out , and the fluid-solid interface Σ. On Γ f in , a velocity field

u f = u in f is prescribed. The outlet boundary Γ f out (t) is divided into four parts, Γ f,(i)
out , i = 1, .., 4 denoting the brachiocephalic trunk, left common carotid artery, left subclavian artery and distal part of the thoracic aorta, respectively. On each Γ f,(i) out , an outlet pressure p i is applied by coupling the Navier-Stokes equations to a three-element Windkessel model. More precisely,

p i = P w i + R p,i Q i ,
where Q i is the volumetric flow rate,

Q i = Γ f,(i) out u f • n,
and P w i is the solution to the ordinary differential equation

C i dP w i dt + P w i R d,i = Q i .
The parameters R d,i , R p,i and C i represent the distal resistance, the proximal resistance and the capacitance of the distal vasculature, respectively. The values used in this work can be found in (Moireau et al 2011), Table 1. Lastly, initial conditions must be provided for the velocity and the displacement in the solid, for the velocity and the domain displacement in the fluid, and the four Windkessel pressures P w i .

Computational aspects

FSI Master

Fluid Solid

σ f f f • n f σ f • n f ( y s | Σ , u s | Σ )
Fig. 2 "FSI-master code" organization where fluid and solid are handled by separate processes which exchange physical quantities at the interface Σ

In practice we adopt a "multi-code strategy". Concerning the FSI direct simulation in itself, we solve the interaction using a partitioned method inspired from domain decomposition strategies [START_REF] Fernández | A coupled momentum method for modeling blood flow in three-dimensional deformable arteries[END_REF]) and a semi-implicit coupling scheme [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]: at every time step, the advectiondiffusion-ALE stage is solved only once for the sake of efficiency, while the pressure Poisson equation is strongly coupled to the solid, for the sake of stability. This non-linear problem is solved with a Newton algorithm. Besides its sheer effectiveness, this scheme also allows to manage both solid and fluid subsystems in their specific implementations with a so-called "FSI-master code"see Figure 2 -that exchanges the Dirichlet or Neuman information requested by each individual code. Therefore, at each iteration of the Newton algorithm the FSImaster:

receives a fluid "variational residual", corresponding to the load on the interface, and sends it to the solid; -receives solid displacements and velocities at the interface, and sends them to the fluid.

In our case, the fluid and structure codes were developed independently. The solid model is solved using a classical mid-point time scheme with a P 1 -Lagrange finite element discretization. The fluid system is also discretized using P 1 -Lagrange elements and is stabilized with a streamline diffusion term. The time scheme compatible with the partitioned strategy is a Chorin-Teman scheme. On a standard workstation (8 cores, 2.5 Ghz) the model solution for a complete heart-beat with a time step of 10 -3 s meshes of about 40, 000 vertices -in total for the fluid and the solid -takes about 8 hours, which gives a reference T ref timing to which we can compare estimation computational times.

Cardiovascular data

In this study, we focus on a subject presenting significant thoracic wall motion. The wall motion was acquired in a typical clinical data acquisition setup, as previously presented in (Moireau et al 2011). The data -Figure 3(a) -was obtained using a Computed Tomography (CT) angiographic scan -a 64-row Multi-Detector CT scanner (Somatom Sensation Cardiac 64, Siemens Medical Solutions). Ten phases of a 512x512x299 voxel image were reconstructed through the cardiac cycle with temporal resolution of 90 ms, in-plane resolution of 0.62x0.62 mm and a slice thickness of 1.00 mm. From these images, ten 3D segmentations were generated using a level-set method [START_REF] Caselles | Geodesic active contours[END_REF] illustrated in Figure 3(b). We then created the corresponding triangular meshes modeling the lumen boundaries [START_REF] Xiong | Simulation of blood flow in deformable vessels using subject-specific geometry and assigned variable mechanical wall properties[END_REF]Moireau et al 2011). These surface meshes of about 6,000 nodes and 12,000 triangles are presented in Figure 3(c). The computational mesh was obtained from the segmentation corresponding to the diastolic phase after trimming the inlet and outlet to produce flat surfaces suitable for boundary conditions specification in the flow domain -see Figure 3(d). Then, the fluid domain was discretized into a P 1 -Lagrange tetrahedral mesh. Next, we generated the tetrahedral solid mesh by extruding the lumen boundary triangular mesh with a varying thickness through the model. For example, the thickness of the three smaller arteries -brachiocephalic trunk, left carotid and left subclavian arteries -is set to 1 mm, while we use 2 mm in the rest of the aorta with a rapid smooth variation at the junctions. The final computational mesh consists of two conforming subdomains, namely, the fluid domain of approximatively 23,000 nodes and 110,000 tetrahedra, and the solid domain of approximatively 18,000 nodes and 24,000 tetrahedra as reproduced in Figure 3(e). We finally divided the exterior boundary of the solid mesh into several regions that have distinct external tissue support characteristics -see Figure 3(e). Each region requires the specification of different numerical values for the coefficients of the visco-elastic boundary conditions. In particular, we identified a "spine region" by locating the area between the origin of the intercostal arteries. Adjacent to this region, we defined a "spine vicinity" region containing the elements of the interface having their nodes within 1 cm of geodesic distance [START_REF] Peyré | Heuristically driven front propagation for fast geodesic extraction[END_REF] from the spine region. The smaller arteries also define a specific segment for the external tissue support prescriptions, and we further assign separate parameter values on the solid inlet and outlet to set up the artificial truncation boundaries. The rest of the exterior boundary will be considered as a region where the external support is assumed to be light at this stage. This region decomposition will be discussed from a physiological standpoint and possibly refined as described in Figure 7 and in the Discussion section.

As additional patient-specific data, we have recordings of diastolic and systolic pressures. In order to generate the prestress state, we need to associate a stressfree configuration with the solid reference mesh. This is achieved by solving the inverse problem consisting of finding a stress-free configuration that corresponds to a given target deformed geometry under a prescribed loading. In our case the target geometry is the solid reference mesh, and we take the measured diastolic pressure as the prescribed loading.

The above-described sequence of segmented meshes is also used to define the displacements y b enforcedweakly in our approach -by the heart beats. As described in (Moireau et al 2011), a motion is extracted from the segmentation by identifying the sinotubular junction -separating the sinus of Vasalva from the ascending aorta -and the left coronary artery attachment, see the red circle and red dot in the segmentations of Figure 3(c). In fact, by defining the barycenter of the sinotubular junction, its normal vector and the projection of the left coronary artery attachment point on this section, we can extract a rigid motion. More complex enforced displacement patterns can also be considered by further localizing the right coronary artery attach- ment, but we found that it was not segmented with sufficient accuracy in our ten phases.

Segmentations

The last type of data used in this framework consists of the inlet flows. From physiological considerations on the subject considered, we used an idealized flow waveform mapped to a flat velocity profile to obtain u in f on the aortic root section Γ f in . The time variation of the flow is representative of a typical aortic root flow, with a systolic phase spanning over one third of the cardiac cycle and no flow during diastole. This generated a typical cardiac output of 5 liters/minute and a heart rate of 67 beats per minute (Moireau et al 2011).

Sequential estimator

In this section, we describe our estimation approach by starting with the estimation of the state -meaning the system trajectory -then we extend the strategy to parameter identification. At the end of this section, our methodology is validated with a simple example created with synthetic data.

A sequential estimator -or observer -is a system that recursively uses the data throughout the simulation to retrieve the system trajectory. In the next section, we will see that this idea extends to parameter identification by allowing the estimated parameters to evolve during the simulation to pursue a final value starting from given a priori values. We commonly use a specific notation to distinguish the estimator system from the original system (1). Hence, if we denote by x(t) = (y s (t), u s (t), u f (t), p(t)) the system state variables -and in the next section also the additional variable θ representing the parameters to be identifiedwe now simulate (ŷ s (t), ûs (t), ûf (t), p(t)) and an additional θ following a modified version of the dynamics (1). The modification -classically called correction -uses the discrepancy between the computed system and the available observation in order to counteract the various sources of uncertainties.

State observer

Discrepancy measure and state correction

To introduce this observer, let us first specify how we can compare model and data through a discrepancy measure involving a computed state x = (y s , u s , u f , p) and a given observation. In fact, assuming sufficient regularity for the real solution and its resulting observation, we can compute at a time t k a discrepancy oper-ator as the signed distance field dist S k :

(L 2 (Σ)) 3 → (L 2 (Σ)) 3 x(ξ) → dist S k (x(ξ)) = dist S k (x(ξ))n S k (x(ξ))
between the simulated system and the available observed surface S k with n S k (x(ξ)) the vector normal to S k at the projection point. Note that the signed distance allows to determine on which side of the surface the model lies. At every other time, we can then interpolate these signed distance fields to compute appropriate distance fields consistent with the hypothetical position of the unavailable surface S(t). For example using a linear interpolation we define dist :

L 2 (Σ) × [0, T ] → L 2 (Σ) (x(ξ), t) → dist(x(ξ), t) = t k+1 -t ∆T dist S k (x(ξ)) + t-t k ∆T dist S k+1 (x(ξ)) × t k+1 -t ∆T n S k (x(ξ)) + t-t k ∆T n S k+1 (x(ξ)) , if t ∈ [t k , t k+1 ]
Based on these distance fields, we can build a discrepancy measure defined for example as the L 2 (Σ)norm of the distance on the surface

D(y s , t) = meas L 2 (Σ) {S k }(x, t) = Σ |dist (ξ + y s , t)| 2 dΓ, (4)
which is null for the pursued actual displacement field, up to measurement errors and a consistency term deriving from the linear interpolation between two consecutive surfaces, i.e.

∆T (t) = D(y ref s , t) = O(∆T ) + noise. Therefore, if this discrepancy vanishes -up to the error ∆T -then the given state is consistent with the observations. Note that it does not imply that the simulated system exactly represents the true state, but it indicates that the observation cannot provide any new quantitative information on the quality of our numerical approximation with respect to reality. By contrast, when "D = 0", the simulated trajectory should be corrected with a term that nudges the simulated system in the direction of the actual trajectory. This nudging direction is obtained by differentiating the measured distance with respect to the state. Hence, if we define the interpolated normal vector as

n {S k } (ξ + y s , t) = t k+1 -t ∆T n S k (x(ξ)) + t-t k ∆T n S k+1 (x(ξ)) ,
which is directly related to be the tangent field of the distance -see [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF] -we can define a discrepancy operator in a variational form on the space of admissible displacements V s by considering

P L 2 (Σ) {S k }(x, t)(v s ) = Σ dist(ξ + y s , t)n {S k } (ξ + y s , t) • v s dΓ, (5) for all v s ∈ V s .
The distance (4) considers an L 2 -norm on the interface Σ but other norms are possible. In particular, from [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]Bertoglio et al 2012a) we have seen that an H 1 2 -norm defined with a lifting operator Ext s on the solid domain is a more effective discrepancy measure for the state observer. We then introduce meas

H 1 2 (Σ) {S k }(x, t) = Ω s 0 Ext s (dist (ξ + y s , t)) 2 E s l dΩ, (6) with . 2 E s l
the norm associated with the linearized elastic energy in the solid, namely equivalent to the H 1norm. Therefore its corresponding discrepancy operator operator becomes

P H 1 2 (Σ) {S k }(x, t)(v s ) = Ω s 0 Ext s (dist(ξ+y s , t)), Ext s (v s •n {S k } )n {S k } E s l dΩ,
for all v s ∈ V s . This operator can be shown to be equivalent to

P H 1 2 (Σ) {S k }(x, t)(v s ) = Ω s 0 Ext s (dist(ξ + y s , t)), v s E s l dΩ, (7) using the extension characterization ∀v s such that v s Σ = 0, Ω s 0 Ext s (dist(ξ + y s , t)), v s E s l dΩ = 0.
The discrepancy is used as a correction term in the equation relating the solid displacement time derivative to the solid velocity to ensure the stabilization of the state estimator on the actual trajectory as fully analyzed -with linearization arguments -in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. This gives in a weak formulation

∀v s ∈ V s , Ω s 0 ∂ŷ s ∂t , v s E s l dΩ = Ω s 0 ûs , v s E s l dΩ + γ Ω s 0 Ext s (dist(ξ + ŷs , t)), v s E s l dΩ,
where γ represents a gain parameter to adjust the nudging. Hence, for the complete system (1) in strong formulation we substitute (ŷ s , ûs , ûf , p) for (y s , u s , u f , p) in the solid formulation (1b) with the additional identity

∂ŷ s ∂t = ûs + γExt s (dist(ξ + ŷs , t)), (8)
replacing the canonical identity ∂y s ∂t = u s . The constant γ must be carefully calibrated since it represents a trade-off between enhanced stability and the undesirable amplification of the measurement noise, see [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF][START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]Bertoglio et al 2012b,a).

4.1.2 Computational aspects FSI Master Fluid Stabilized Solid Solid Discrepancy Measure Fluid Discrepancy Measure Processed Observations S k , S k +1 dist( ξ + y s | Σ , t )

Observation

Managers:

Fig. 4 "FSI-master code" with additional "observation managers" to compute the state observer

Based on the FSI code presented in Figure 2 we handle the computations relative to the observations by interfacing the solid code with an "observation manager" that loads the observed surfaces and computes the discrepancy with a mesh-to-mesh distance algorithm (Baerentzen and Aanaes 2005), as was done in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. When using the state observer model, the modification of the velocity to time-differentiated displacement identity into (8) induces an additional computational cost because velocities cannot be directly eliminated any longer, which doubles the size of the system, hence roughly multiplies the inversion time in the solid by a factor 4, as explained in detail in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. Since the solid is not the most CPU-intensive part in the complete FSI model this finally only roughly doubles the solution time of the complete problem, meaning about 2T ref . Note that this modified architecture presented in Figure 4 can also be supplemented by an analogous "observation manager" on the fluid side if, for example, time resolved velocity or pressure profiles are available -see the end of the discussion for more details.

Parameter estimation

We present here a sequential identification strategy originally proposed in (Moireau and Chapelle 2011b,a), itself inspired from [START_REF] Pham | Filtres de Kalman singuliers évolutifs pour l'assimilation de données en océanographie[END_REF][START_REF] Pham | Stochastic methods for sequential data assimilation in strongly nonlinear systems[END_REF]) and already evaluated in a fluid-structure context in (Bertoglio et al 2012b) for wall stiffness identification where the reader can find the additional specific technical material.

In this work, we concentrate on the estimation of the elastic tissue support k s (ξ) representative of the regionalized external tissue tethering. We recall that the available observations are the lumen boundaryi.e. the internal boundary of the wall corresponding to the fluid-structure interface Σ -whereas we want to identify parameters defined on the external boundary of the walls Γ s . We must therefore assume that we are confident in the value of the wall stiffness, and we will tackle the question of a joint estimation of external support stiffness and wall stiffness in a future work using additional measurements to differentiate the contributions in the motion of these two types of parameters.

The parameter k s (ξ) is a spatialized field and we assume that we can discretize it with a vector θ ∈ R p of small dimension -typically p < 100 and in this work p 10. Complex spatial representations can be considered in this respect, as done for example in [START_REF] Xi | Myocardial transversely isotropic material parameter estimation from in-silico measurements based on reduced-order unscented Kalman filter[END_REF], but we choose here to simply divide Γ s into p regions within each of which k s (ξ) will be considered as constant. Clearly, the choice of the region division should be consistent physiologically, and we will discuss our estimation results in this light.

Sequential parameter estimator description

The identification procedure in a sequential framework consists in adjusting the parameters during the simulation to decrease the discrepancy measure. In other words, the parameters of interest are no longer constant in the estimator and instead evolve with a dynamics of the type

∂ θ ∂t = KP * {S k }(x, t)(∇ θ y s ),
where K denotes a filter gain operator, P * {S k } a discrepancy operator which remains to be specified with respect to its associated norm which in general can differ from that used in state estimation, and ∇ θ y s quantifies the sensitivity of the displacement field with respect to the parameters [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF]Moireau and Chapelle 2011b). Then the state dynamics should be accordingly corrected relatively to the original dynamics by using the parameter sensitivity, resulting into

∂ x ∂t = ∂x ∂t + (∇ θ x) • ∂ θ ∂t ,
where " ∂x ∂t " is an abuse of notation representing the original model equations.

The small dimension of θ allows us to use optimal filtering to define a filter gain K that is best adapted to a given choice of discrepancy measure. In our case, we choose for simplicity the meas L 2 definition. The computation of this gain can be approximated by a Reduced-Order Extended Kalman Filtering approach (ROEKF) [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] whereby the sensitivity quantities are computed by differentiating the model equations. Alternatively (Moireau and Chapelle 2011b,a) have proposed a Reduced-Order Unscented Kalman Filtering (ROUKF) strategy where, instead of resorting to tangent sensitivity computations, the non-linear impact of the parameters on the state is determined from sampling points representing the parameters uncertainties around the current estimate. To that purpose, we rely on an estimation of the parameter covariance whichstarting from an initial a priori covariance -will evolve in time, fed by the inverse of the observability Grammian -also called here identifiability Grammian-of the problem [START_REF] Simon | Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF]Moireau and Chapelle 2011b). In this case, we manipulate what we refer to as "uncertainty sensitivities" which map uncertainty directions in the parameter space to associated directions in the state space.

In order to fully summarize the ROUKF approach, it is convenient to write the final equations after space and time discretization of the model. Let us denote by X n the spatial discretization of the state at time t n . After discretization of the complete fluid-structure interaction problem we can introduce a forward operator A n which computes X n+1 starting from X n , i.e.

X n+1 = A n (X n , θ).
By considering the additional state correction we modify A n into a new operator A n γ such that the state observer follows the discrete dynamics

Xn+1 = A n γ ( Xn , θ, {S k }).
The purpose of the state estimator is to stabilize the simulated trajectory on the actual trajectory. Of course, this is partly artificial as long as the parameters themselves have not been adequately identified, since the parameter error then induces a persistent error source in the model equations. Nevertheless, for parameter identification purposes, remaining in the vicinity of the true trajectory is essential in order for parameter sensitivities to be representative. Then, the parameter identification problem is more likely to be well-posed -and we can prove this in a linearized configuration [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF]. This motivates replacing the classical fluidstructure model A n by the observer A n γ . Note that in A n γ the distance field to the observation surfaces is discretized into a quantity denoted by D( Xn , t n ). Discretizing the L 2 (Σ)-norm with a finite element matrix M Σ we have for example

P L 2 (Σ) {S k }(x n , t n )(v s ) = V s ∂D ∂X M Σ D( Xn , t n ),
for all test functions v s ∈ V s and V s the vector of the corresponding degrees of freedom. The same logic applies to compute P

H 1 2 (Σ)
{S k }(x, t)(v s ), and we refer to (Bertoglio et al 2012a) for technical details.

The ROUKF modifies this dynamics by allowing θ to evolve in time, and we therefore write θn . The algorithm follows a prediction-correction formalism and is based on sampling points for computing the sensitivities. Given adequate sampling rules, we precompute so-called unitary simplex sigma-points [START_REF] Julier | Reduced sigma point filters for the propagation of means and covariances through nonlinear transformations[END_REF]

I [i]
, 1 ≤ i ≤ p + 1 as a discretization of the unit ball in R p . We associate with each sigmapoint a corresponding weight α i = 1 p+1 , and collect the weights in a diagonal matrix D α . Then, assuming that the we expect a variation of the parameters given by a covariance ellipsoid around a computed value θn , we can sample particles around θn using I [i] and the standard deviation associated with the covariance ellipsoid. Hence, from computed Xn

+ ∈ R N , θn + ∈ R p with additional operators U n ∈ SL p (R) (symmetric), L n x ∈ M p,N (R), L n θ ∈ GL p (R) (invertible), we compute in three steps 1. Sampling:        C n = (U n ) -1 Xn + [i] = Xn + + L n x • C n • I [i] , 1 ≤ i ≤ p + 1 θn + [i] = θn + + L n θ • C n • I [i] , 1 ≤ i ≤ p + 1 2. Prediction:        Xn+1 - [i] = A n γ ( Xn + [i] , θn + [i] , {S k }), 1 ≤ i ≤ p + 1 Xn+1 -= E α ( Xn+1 - * ) θn+1 -= θn+1 + 3. Correction:                                          L n+1 x = [ Xn+1 - * ]D α [I * ] L n+1 θ = [ θn+1 - * ]D α [I * ] D n+1 [i] = D( Xn+1 - [i] , t n+1 ) L n+1 d = [D n+1 * ]D α [I * ] U n+1 = 1 + w -1 L n+1 d M Σ L n+1 d ∈ M p θn+1 + = θn+1 - -w -1 L n+1 θ (U n+1 ) -1 L n+1 d M Σ E α (D n+1 * ) Xn+1 + = Xn+1 - -L n x (L n+1 θ ) -1 ( θn+1 + -θn -)
where for every set of particles

X [i] , 1 ≤ i ≤ p + 1, [X * ] ∈ M N,p+1
denotes the concatenation of the (p+1) column vectors arranged side by side, and

E α (X * ) = 1≤i≤p+1 α i X [i] computes the empirical mean of the X [i]
. The scalar w -1 is a coefficient playing a similar role to γ in the state observer to moderate the parameter correction depending on the measurement noise. In other words, we take into account the observation noise amplitude by considering weighted norms in the choice of our discrepancy measure in the parameter identification process, meaning that we consider P w -1 L 2 {S k }(x, t), where w -1 L 2 as a subscript means that we use the L 2 -norm weighted by a factor w -1 .

As described before, we see that based on L n+1 X and L n+1 θ which are two operators gathering particles information, a sensitivity operator L n+1 X (L n+1 θ ) -1 is computed to correct the state dynamics after the parameter update. Furthermore, the operator L n θ U -1 n L n θ can be interpreted as a covariance representation of the parameter variation. Therefore we recognize in L n θ C n a standard deviation used for generating the sampling points starting from the unitary sigma points

I [i] .

Algorithm aspects

From a practical standpoint, the algorithm can be described as follows. At each time step t n , for a given estimated state, parameters, parameters covariance and uncertainty sensitivities, 1. Sampling:

-Generate p + 1 parameters θn + [i] sampling the covariance ellipsoid around the current estimated value; -Infer the corresponding state sampling points Xn +

[i] based on the uncertainty sensitivities; 2. Prediction:

-For each sampling point ( Xn

+ [i] , θn + [i]
), solve a onetime step model analysis starting from Xn +

[i] at

FSI Master

Discrepancy Measure Fluid Discrepancy Measure and infer the new uncertainty sensitivities; 3. Correction:

#1 #1 #1 #1 #1 #1 #1 #1 #1 #1 #1 #1 #1 #1 #
-Gather the discrepancies. Combined with the uncertainty sensitivities compute the optimal gain; -Modify the estimated parameters and parameters covariance in consequence; -Modify the estimated augmented state ( Xn+1 + , θn+1 + ) by using the sensitivity.

In order to build on the modular architecture of our FSI code, we implement a master code -on top of the FSI-master -in charge of the identification process as illustrated in Figure 5. This "Identification-master" supervises the complete estimation simulation time-step after time-step. It starts by launching in parallel as many FSI-master processes as the identification procedure requires -in our case p + 1 processes to represent the parametric covariance ellipsoid in a space of dimension p. Then each FSI-master receives state-parameter sampling points from the Identification-master to be dispatched to the solid and fluid in order to solve a single time-step evolution. The results are gathered by the FSI-master which forwards them to the Identificationmaster for processing. This strategy has a very limited impact on the original FSI code, and each sampling point time-step simulation can be independently run in parallel. Hence, a complete heart beat simulation with parameter identification is only marginally longer than the previous 2T ref with a CPU-overhead almost entirely due to communications.

Synthetic test case for algorithm evaluation

Before moving to the actual results section where we will use the clinical data presented in Section 3, we first validate our complex estimation framework with a simplified problem using synthetic data. We consider a straight tube of dimensions comparable to the descending aorta with a 2 cm radius and 10 cm length. On this simplified geometry, we define 4 regions. The first one represents the "spine" region surrounded by 2 "vicinity" regions, and the fourth region corresponds to the "opposite" to the spine as represented in Figure 6. We define in each region a specific parameter value for the stiffness coefficient of the external tissue support condition in the form

k s (ξ) = 2 θi , on Ω i 1 ≤ i ≤ 4.
The corresponding mesh of about 1,000 nodes and 4,000 tetrahedra for the solid, and 2,000 nodes and 10,000 tetrahedra for the fluid, was designed to be representative of the target system without undue com- putational complexity to efficiently evaluate the identification algorithm. In this regard, the input flow was defined from the flow extracted from the real aortic simulation in (Moireau et al 2011) at the junction between the aortic arch and the descending aorta, and the Windkessel element parameters are also taken from (Moireau et al 2011). By setting θ ref = ( log(10) log(2) , 0, 0, -log(10) log( 2) ) corresponding to 2 θ = (10, 1, 1, 0.1) we generate synthetic surfaces at every time step. Then the identification procedure is performed starting from the a priori θ(0) = θ 0 = (0, 0, 0, 0) . The parametric covariance is set to identity and the noise covariance is purposely low since the generated data are noise-free, in order to only evaluate the identification strategy. The parameter evolution during the sequential identification procedure is presented in Figure 6, and we conclude to the very fast convergence of the estimated parameters to their reference values.

Clinical data results

In this section we will present estimation results obtained with the above-described clinical data when applying our proposed methodology. In order to provide quantitative comparison references we will compute some numerical errors based on distances between contours in several sub-regions of the fluid-solid interface, namely, the lumen contour. These regions are depicted in Figure 7, which also illustrates the comparison and distinction of the previously introduced boundary support regions subdividing the outer wall. Here, we first compare in Figure 8 contours extracted from the image sequence to the initial configuration of the model -obtained after regularization of the extracted contour in the first image -in order to provide an estimate of the apparent motion of each sub-region, both in L 2 and L ∞ norms. This should give valuable comparison references for the next results of this paper

We also reproduce in Figure 9 some error plots already given in (Moireau et al 2011) and providing the distances between the segmented contours and the lumen boundaries of the model simulated with homogeneous Dirichlet conditions prescribed at the inlet and outlet as only essential boundary conditions for the wall, namely, without tissue support. This shows that the errors induced by this naive modeling setup are of the same order as the motion itself, and even much larger in some areas such as the spine and its vicinity.

State estimation results

In this section we illustrate the efficacy of the state estimator in robustly tracking the actual trajectory. In the absence of absolute ground truth, we consider some reference direct simulations produced with parameter values manually calibrated to represent the external tissue support as in (Moireau et al 2011), where these support conditions were shown to significantly improve the distance errors compared to the simplistic boundary conditions considered in Figure 9.

We now compare the standard simulation with boundary support derived from System (1) and the corresponding state estimator -where we modify (1b) using (8) -for an initial condition error that should be adequately tackled by the estimator. To that purpose, we generate the initial condition by running the standard simulation during one cardiac cycle and saving the last configuration. In the results presented in Figure 10, we compare the results obtained with increasing values of the gain γ characterizing the nudging effect, in order to obtain a satisfactory trade-off between the dynamics stability improvement and the data noise amplified by the gain. The first row shows the error results obtained with a zero gain, namely, a direct simulation. The second row displays the errors corresponding to a finite gain value γ 0 = 3(cgs) and we observe a decrease of the distances in all regions as expected, even if the parameter values are only manually calibrated. The third row corresponds to γ = 5γ 0 , but convergence in the simulations was only obtained after restricting the observation surfaces in the estimator to disregard an area of 1 cm extent along the vessel, starting from the inlet ring on which large displacements are prescribed. This was necessary to avoid difficulties in the ALE mesh motion computation induced by large prescribed displacements directly adjacent to regions where strong estimation corrections are applied. In fact, for γ = 10γ 0 such difficulties with the ALE mesh appear with the estimation corrections alone -i.e. away from any boundary condition -which draw the boundaries toward segmented surfaces featuring limited regularity. Note additionally that using a very large gain also reduces the system identifiability, see the Discussion section below. Finally, for all these reasons we select γ = γ 0 for all our subsequent estimation computations.

Parameter identification

We now consider a joint state and parameter estimation of 4 regions of interest -see Figure 7, 2 nd left. The first three regions are the "spine" region, one "spine vicinity" region surrounding the spine, and a "light support" zone aggregating the region opposite to the spine with the arch and the ascending aorta. In these three regions we identify the external tissue support stiffness representing the external tethering. The last region is made of the outlet surface ring at the end of the descending aorta. In this region we also want to estimate the stiffness parameter of the tissue support boundary condition which, here, artificially represents the rest of the arterial tree. We gather the region parameters corresponding to the 4 regions ("spine","spine vicinity","light support","outlet") in a vector θ = (θ i ) 1≤i≤4 , and define as we did for our validation test case

k s (ξ) = 2 θi , in Γ s i 1 ≤ i ≤ 4.
In this real case estimation, we are primarily interested in the external support of the aorta itself, hence we use prescribed (small) values of k s (ξ) as calibrated in (Moireau et al 2011) for the "small arteries". Furthermore, the artificial inlet support is large enough to assume that the heart motion is almost imposed as a Dirichlet boundary condition and does not necessitate an actual identification. We also refer to Section 6 for further discussions on whether or not to include these quantities in the estimation.

We then consider two starting a priori choices: 1. Firstly, we choose θ 0 = 0 for all regions and the covariance is diagonal with

P 0 (1, 1) = 1.25, P 0 (i, i) = 0.125, {i = 2, 3}, P 0 (4, 4) = 0.625,
which reflects our a priori on the expected variations of the parameters θ 1 and θ 4 for the "spine" and "outlet", respectively. 2. Secondly, we use θ 0 = θ manual , the value manually calibrated in (Moireau et al 2011) and

P 0 = α1, with α = 0.125,
to represent the fact that we assume the parameters should vary around the a priori value with the same expectation.

Considering the observation noise amplitude we choose a weighted L 2 -norm in the choice of our discrepancy measure in the parameter identification process and set w = 10 -2 /180 calibrated to represent the amplitude of the measurement noise covariance (about 1 mm 2 in cgs) divided by the surface of the regions measuredof about 180 cm 2 -to cancel the surface extent effect in the norm. Note that in practice after discretization in time of the model, this term is multiplied by δt ∆T the ratio of the time discretization δt time-step versus the data time sampling ∆T to take into account the difference of time-sampling between the simulation and the data. Indeed, it is logical to decrease the covariance when the data sampling increases, and when we increase the model accuracy for a given data time-sampling. The results in Figure 11 show the evolution of the parameters in the course of the sequential estimation process, for the above two choices. The solid line represents the estimated mean value, whereas the corresponding colored bands show the estimated standard deviation around the parameters used by the algorithm to evaluate the sensitivity matrices. The band shrinking around the mean value indicate a strong confidence in the identification. The parameters associated with the "spine" region and the "outlet" converge toward positive and rather large values, while the "'spine vicinity" coefficient remains close to 0 indicating a light spine tethering in this region. However, we see oscillations in the parameter evolution of the "light support" region reflecting a lack of convergence in the identification process for this particular parameter. The observations are consistent for two starting a priori choices and we further compare the respective estimation results in Section 6.

In the absence of absolute ground truth, following [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model[END_REF] we validate these first results with a second estimation procedure using more regions. In particular, the descending aorta is subdivided into 3 segments -see Figure 7 left -leading to 3 "spine" subregions, 3 "spine vicinity" subregions and 3 "opposite to spine" subregions, now distinct from the separate "arch" region. The initial covariance is set proportional to identity -i.e. P 0 = α1, with α -0.067 -since we start from the estimated values obtained with the previous 4 regions. Figure 12 shows the evolution of these 11 parameters with a stable convergence for all values except for the "arch" region which continues to oscillate.

We finally run a direct simulation -namely without the identification procedure -using the estimated values and compare the distances with the results obtained with very light external tissue support in these 4 regions -i.e. θ 0 = 0 -or with the calibrated values of (Moireau et al 2011) -i.e. θ = θ manual . The computed distances are presented in Figure 13, and we also display in Figure 14 a comparison of the simulated contours with CT images at peak systole, as well as computed fluid pressures and velocities at the inlet and outlets.

Discussion

Successfully identified parameters -The results demonstrate that we obtain stable and consistent identified parameters for the "spine" region and its vicinity, throughout the various estimation setups considered above, namely, with different initial parameter and covariance values (Figure 11), and when increasing the number of regions (Figure 12). Note that the final estimated parameter in the "spine" region is in fact smaller than the calibrated value of (Moireau et al 2011), which is found to be justified when closely examining in whereas we see in Figure 8 that this part is not completely fixed, see also Figure 14.

A slightly consistent -although still stablebehavior is observed for the outlet parameter. In this case, as we will see below with Figure 15, identifiability is not an issue since this parameter converges very fast to its synthetic value when using synthetic data. However, this parameter corresponds to an ad hoc modeling ingredient compensating the truncation of the arterial tree, hence it certainly does not have an absolute value independently of other modeling parameters, and in particular its adjustment can be expected to depend on the spine support parameter value.

This validates the double interest of our solid boundary terms in order to represent (1) physical boundary conditions accounting for external anatomical structures, and (2) adjustable artificial boundary conditions modeling the rest of the arterial tree. Moreover, this also indicates that no tissue support is really necessary in the "spine vicinity" and "light support" regions.

Identifiability considerations regarding heterogeneous regions -The parameter associated with the region containing the arch shows an oscillating behavior, typical in sequential estimation of non convergence due to identifiability issues. We claim that this is due to two main reasons. First, the computation of distances between two complex curved geometries like arches is not very informative. In fact, for a straight geometry with a main deformation normal to the surface, the distance to the target object is a good estimate of the displacement error between the current simulation and the target system. However, for the aorta, taking into account the large displacements induced by the heart, the distance computation in the aortic arch does not allow to follow actual material points. Then various configurations can lead to a rather small discrepancy measure and the total identifiability is reduced. Secondly and probably more importantly, the region called "light support region" in this paper in fact consists of heterogeneous components which are: the "opposite to the spine" part, the "arch" in itself and even the ascending aorta. Each region cor-responds to a different motion and behavior and also to a different tissue support -for example the arch is in reality sustained from below by the pulmonary trunk. Hence, by aggregating heterogeneous regions in a single region we decrease the identifiability of each region in itself.

These remarks are corroborated by our estimation presented in Figure 12 where we subsampled the 4 initial regions into 11 regions. In this estimation starting from the converged parameter values obtained with 4 regions, we see that the estimation is stable for the "spine" and "vicinity" subregions. Furthermore the "opposite to the spine" subregions do not exhibit parameter oscillations anymore and these estimated values indicate the expected light support behavior. However, in the last region corresponding to the heterogeneous aggregation of the arch and the ascending aorta the parameter continues to oscillate.

Data time sampling identifiability issues -The accuracy of the estimation should also be interpreted in the light of the limiting identifiably due to the time sampling of the real data which is quite low, i.e. ∆T = 10 2 δt = 0.1T . To illustrate this argument we present in Figure 15 an additional estimation procedure with synthetic data using the real geometry. Here, we generate target surfaces with two different sampling rates: a reasonably high resolution ∆T = 10δt = 0.01T , and a low resolution ∆T = 100δt = 0.1T analogous to the real data resolution. The observations correspond to a model where the external tissue support parameters are those manually calibrated in (Moireau et al 2011), θ ref = θ manual . Then the identification is performed with an initial state error as in Section 5.1, and with zero a priori values for all parameters. For these estimations, we consider a rather small observation noise covariance (w = 10 -2 w real with respect to the value chosen in the real configuration) identical for the two sampling rates, arguing that the observation surfaces are noise-free since they are originally synthetically produced by the model. However, there still exist some errors in the observations due to time sampling and the resulting linear interpolation used in the estimation sequence. The impact of this sampling error is therefore purposely amplified by the small observation noise covariance.

We see that the outlet solid boundary condition is well estimated in absolute value already with a low sampling rate, which confirms its good identifiability. Concerning the spine external tissue support we see that the sensitivity is correct with the two sampling rates even if we start from a wrong a priori. However, absolute value convergence is obtained only with the higher sampling rate. Finally concerning the "vicinity" region and the "light support" region, accurate trends can be globally observed since the "vicinity" region parameter value is in between the "spine" region parameter value and the "light support" region parameter value. However absolute values are incorrect. We argue that it is mainly due to an identifiability issue in the "light support" region -as developed above -which contaminates the "vicinity" region. Note that since the noise covariance is low, all effects are magnified in this illustrative example.

Effect of the state estimator -In fact, conceptually the essential purpose of the state estimator is to remain "sufficiently close" to the actual trajectory during the identification process while the parameters evolve in time. We emphasize that this idea is fundamentally distinct from the objective of making the simulated system simply "fit the data", namely driving the computed discrepancy to zero, which is the main purpose of data-based tracking. In the case of filtering estimation, by contrast, the objective is achieved by concurrently nudging the system "towards the data" and enhancing the stability of the simulated system. This induces some specific constraints in the choice of the state estimator gain, in addition to concerns related to measurement noise amplification as already seen in tracking. Namely, in many types of dynamical systems of interest -as e.g. in solid continuum mechanics and FSI -making the gain parameter arbitrarily large has a detrimental impact on the stability [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF][START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. Moreover, as already discussed and assessed in (Chabiniok et al 2011), using a very large gain also reduces the system identifiability, since the state estimator would then follow the observation regardless of the parameter values, which would compromise the concurrent parameter identification process.

Integration of multiple observations and material parameters -In the work presented here, the observations were given by a set of medical image data segmentations defining the position of the fluid-structure interface Σ through the cardiac cycle. This observation is therefore concerned with a single state variable: the (Eulerian) motion of the fluid-solid interface. However, in a more general setting, one may have to deal with multiple observations of different state variables, such as recordings of pressure and flow waveforms, in addition to the aforementioned motion of the fluid-solid boundary. These different observations of the state must therefore be considered simultaneously in the state and parameter estimation endeavor. Pressure can be acquired noninvasively via a sphygmomanometer, or invasively via a cannula needle placed through an arterial line. Flow waveforms can be acquired non-invasively using ultrasound or phase-contrast magnetic resonance imaging.

As already mentioned, these measurements can directly enter our estimation formalism on the fluid side -see Figure 5. Moreover, in addition to state variables and their directly derived quantities, one may have measurements of material parameters such as vessel wall thickness, stiffness or distensibility. In this regard, it should be noted that parameters are considered by extension as generalized state variables in the joint stateparameter estimation framework [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF], hence these measurements can directly be used in the filtering corrections. However, the integration of certain types of observations would require further extensions of the proposed strategy, especially when corresponding to quantities not instantaneously associated with the current state and parameter variables at each observation time. For instance, in the case of blood pressure, when acquired non-invasively only the bounds of the state variable (i.e., systolic and diastolic pressures) are registered. Another example is given by pressure pulse velocity observations easily measurable non-invasively for the assessment of hypertension. This situation may preclude the estimation "on the fly" since usually a complete cycle of the solution should be produced in order to quantify the pulse velocity.

Simultaneous estimation of stiffness and external tissue support -In this article, we illustrated the application of a sequential framework based on filtering methods to the estimation of the material parameters defining the support provided by the different tissues and organs around a subject-specific model of the thoracic aorta.

Here, a given value for the stiffness of the arterial wall was assumed to be known. In a previous investigation (Bertoglio et al 2012b) the same framework was applied to the identification of the elastic properties of the arterial wall in an idealized fluid-structure interaction model, neglecting the support of external tissues and organs. The simultaneous estimation of arterial wall stiffness and external tissue support has therefore not been attempted thus far. This task will ultimately provide the most valuable information since it would give estimates on parameters that have been known to be predictors of disease (i.e., arterial stiffness) together with the parameters that enable the representation of realistic motion patterns of the arterial structures (i.e., the external tissue support).

In order to simultaneously estimate arterial stiffness and tissue support, richer information on the arterial kinematics is needed. This information could potentially be obtained using newer modalities of timeresolved medical image data to discriminate between the various components of the motion of the arterial tissue: namely, the deformation due to the oscillating pressure field, the motion resulting from respiration, and the cardiac motion. This would, indeed, enable the identifiability of arterial stiffness -via measurements related to pressure-induced contribution -distinctly from the tissue support properties, while cardiac and respiratory contributions would provide directly relevant information on the corresponding external supports. Notable progress has been made in cardiac and respiratory compensation in Magnetic Resonance Imaging using cardiac gated and respiratory navigators (Buerger et al 2011). Further improvements in image data acquisition techniques, together with the sequential estimation framework presented here will undoubtedly bring us closer to the goal of extracting physiologically meaningful mate-rial data from images. Lastly, it is important to note that the final goal of the framework sketched here is to provide estimates not only for the material properties of the arterial tissue, but also on the input and output flow boundary conditions, via proper integration of image data, pressure, flow, derived state quantities and discrete measurements of material parameters. It is through the simultaneous consideration of these various observations that this framework will be capable of automatically tuning the material parameters and inflow and outflow conditions required in the simulation process.

Algorithm effectiveness and estimation benefits -In Figure 13 we observed (1) compared to light tissue support, an improvement of the distances, in particular near the end of the time window and (2) compared to the manually calibrated tissue support a slight improvement all along the simulation. This improvement is further illustrated in Figure 14 where we can see that the contours of the simulation with the estimated values are significantly closer to the segmented contours than for the simulation with the values taken from (Moireau et al 2011). This final validation of the results of the estimation procedure -in particular in the "spine region" and at the outlet -justifies in itself the interest of such a sequential estimation methodology, especially when considering that this estimation procedure can be performed in a time comparable to a classical simulation. In (Moireau et al 2011), the interest of the external tissue support was already clearly established by comparing the distances obtained with external tissue support against the configuration where no tissue support was considered with only homogeneous Dirichlet conditions imposed at the inlet and outlets -recall Figure 9. However, in (Moireau et al 2011) this manual calibration performed required great care and multiple simulations to quantitatively set up the parameters, whereas here the procedure is fully automatic and produces comparatively better results.

The computational complexity of our algorithm should be compared to other state-of-the-art identification procedures. For example, let us consider a variational identification strategy where we minimize a criterion based on the discrepancy measure integrated over time plus additional regularization terms, as was used in the context of blood flows in [START_REF] Martin | Parameter identification for a one-dimensional blood flow model[END_REF] with a very effective minimization algorithm. Even so, the complete procedure for one single parameter value in 1D requires 6 to 10 iterations, each requiring a complete heart beat forward simulation and a backward simulation of the adjoint problem. The adjoint system -being linearis less intensive to solve, notwithstanding the difficulty of storing all the state solutions and discrepancies over time. Overall, the computational time-cost of the complete identification is then at least 10T ref , meaning several days of computations on a standard workstation. Even worse, estimating multiple parameter values as in the present study would lead to many additional minimization iterations which increase the computational costs all the more. Here we recall for comparison that our estimation procedure CPU time is roughly 2T ref irrespective of the number of parameter values to be estimated.

Concluding remarks

We have proposed a complete methodological chain for the identification of boundary support parameters in a fluid-structure vascular model using patient image data, based on well-adapted sequential data assimilation procedures. We subsequently demonstrated the use of this framework in a real life application case in which we successfully performed the automatic characterization of the tethering to the spine and to the downstream vasculature, and most importantly identified the actual regions in which the boundary support modeling is necessary. The simulations performed with the estimated parameter values proved more accurate than with the previous expert manual calibration presented in (Moireau et al 2011). In this realistic example the estimation procedure itself was shown to be computationally effective -comparable to a direct simulation in terms of computation time, indeed -and robust, in particular with respect to parameter initial values. This methodological framework paves the way for complete patient-specific vascular modeling, in which ultimately all the available measurements could be used to estimate the various uncertain parameters needed for the patient-specific simulations. Moreover, some of these parameters are in itself of utmost clinical value for diagnosis purposes, which represents an important additional motivation for estimation. Finally, estimation can also be envisioned as an integrated tool for investigating intrinsic modeling validity, since a lack of convergence of some estimated modeling parametersprovided identifiability has been established by preliminary synthetic trials and measurement errors are adequately controlled -can only be accounted for by modeling shortcomings. This is especially valuable in the context of patient-specific modeling in which a priori testing cannot be used to assess the model validity.

Fig. 1

 1 Fig. 1 Schematic representation of aorta and its surrounding organs. Model domain and boundary conditions

  and a lifting operator Ext f (called extension map) from Σ 0 to Ω f 0 respecting the fluid Dirichlet boundary conditions on the fluid-structure interface. We can then define the fluid deformation mapping

Fig. 3

 3 Fig. 3 (a) Volume rendering of image containing one phase of the cardiac-gated CT data. (b) 3D level set is grown to detect the lumen boundaries. (c) Segmentations of the ten phases. The red-circle represents the sinotubular junction and the red dot the left coronary artery attachment. (d) Triming and smoothing of the first frame in order to generate the computational mesh (e)

Fig. 5

 5 Fig.5Complete identification code where the "ROUKF master" receives states, parameters and discrepancies and sends updated states and parameters

Fig. 6

 6 Fig. 6 Synthetic case corresponding to an idealized descending thoracic aorta. Estimation of 4 parameters (spine, vicinity left, vicinity right, opposite) with synthetic noise-free surfaces

Fig. 7

 7 Fig. 7 Anatomical schematic (second right) underlying the hierarchical tissue support division of the outer wall (to the left) and the division of the fluid-structure interface for error computations (to the right)

Fig. 8

 8 Fig. 8 Apparent motion analysis -L 2 and L ∞ (left and right resp.) norms of distances from the segmented lumen in each image with respect to the initial configuration of the model -Distances versus time plots in cgs units (for comparison, vessel radius is approximately 1 cm)

Fig. 9 L

 9 Fig. 9 L 2 and L ∞ (left and right resp.) norms of distances from the simulated fluid-structure interface of a model with fixed inlet and outlet as only Dirichlet boundary conditions to the segmented lumen surfaces in the image sequence -Distances versus time plots in cgs units (for comparison, vessel radius is approximately 1 cm)

Fig. 10 L

 10 Fig. 10 L 2 and L ∞ (left and right resp.) norms of distances from the segmented lumen boundary to the simulated model fluid-structure interface using tissue support with additive state observer and initial condition error -Distances versus time plots in cgs units (for comparison, vessel radius is approximately 1 cm)

Fig. 11

 11 Fig. 11 Time identification of the 4 parameters using real data -Standard deviations diagonal coefficient (L n θ C n ) i,i for i = 1..4 are represented by a colored band around the mean value

Fig. 12

 12 Fig. 12 Time identification of the multi-regions parameters (11 values) using real data L 2

Fig. 13 L

 13 Fig. 13 L 2 and L ∞ (left and right resp.) norms of distances from the segmented lumen to the simulated fluid-structure interface using the heart motion only (top) with additive manually calibrated tissue support (middle) and estimated tissue support (bottom)

Fig. 14

 14 Fig. 14 Left: comparison of segmented (dashed green) vs. simulated (blue with estimated values and pink with manually calibrated values) lumen contours at peak systole -Right: computed fluid pressures and velocities at inlet and outlets

Fig. 15

 15 Fig. 15 Time identification of 4 parameters with synthetic data, data sampling ∆T = 100δt = 0.1T (left) and ∆T = 10δt = 0.01T (right). Dashed lines: Target values; Solid lines: Estimated values

  1

								( y s , u s	,	, u y , p	)
							FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master
		FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master FSI Master #j				FSI Master #p+1	FSI Master #p
	Fluid Fluid	Stabilized Solid	Fluid Fluid	Fluid	Stabilized Solid	Fluid	Stabilized Stabilized Stabilized Solid
		Fluid Fluid Fluid Fluid Fluid		Solid Solid Solid	Discrepancy Discrepancy Discrepancy	Discrepancy Discrepancy	Solid Solid Solid Solid
	Fluid	Discrepancy Discrepancy Discrepancy Discrepancy Discrepancy Discrepancy Discrepancy	Measure Measure Discrepancy Discrepancy Discrepancy Measure Measure Discrepancy Discrepancy	Measure Measure Measure	Discrepancy Discrepancy Discrepancy Discrepancy Discrepancy Fluid Discrepancy Discrepancy	Discrepancy Discrepancy Discrepancy Discrepancy Discrepancy Solid Discrepancy Discrepancy
	Measure Measure Measure Measure Measure Discrepancy Discrepancy Discrepancy Discrepancy Discrepancy Measure Measure Measure Measure Measure Measure Measure Discrepancy Discrepancy Discrepancy Solid Fluid	Fluid	Measure Measure Measure Measure Measure Discrepancy Discrepancy Discrepancy Measure Measure Discrepancy Discrepancy Discrepancy Discrepancy Discrepancy Measure Measure Measure Measure Measure Measure Measure Measure Measure Solid Solid Fluid
	Measure Measure Measure Measure Discrepancy Discrepancy Measure		Discrepancy	Measure Measure Measure Discrepancy Discrepancy Measure Measure Measure Measure Measure Measure Discrepancy
		Measure	Measure		Measure	Measure	Measure	Measure
							Processed Observations

Identification Master based on Reduced Order Unscented Kalman Filter


Laurianne Georgeton
email: laurianne.georgeton@univ-paris3.fr

Nicolas Audibert

Cécile Fougeron
email: cecile.fougeron@univ-paris3.fr

ROUNDING AND HEIGHT CONTRASTS AT THE BEGINNING OF DIFFERENT PROSODIC CONSTITUENTS IN FRENCH

Keywords: segmental contrast, prosody, French

Variations in the acoustic manifestations of the rounding and height contrasts in French vowels are investigated in different prosodic positions. Four speakers produced sentences containing the vowels /i, e, a, y, ø/ at the beginning of different prosodic constituents: Intonational Phrase, Accentual Phrase and Word. Acoustic cues for the rounding contrast (F2, F3, F3-F2) and for the vowel height contrast (F1) are found to be enhanced in Intonational Phrase initial positions.

INTRODUCTION

This study investigates the influence of prosodic positions on the acoustic realization of speech sounds. Variation depending on prosodic position has been described in different languages and for different types of segments. So far, variation at the left edge of a constituent (so-called initial strengthening) has been observed mostly on consonants for lingual, glottal or jaw articulation [START_REF] Cho | Effects of prosody on articulation in English[END_REF], [START_REF] Fougeron | V ariations articulatoires en début de constituants prosodiques de différents niveaux en français[END_REF]. Variation in the realization of vowels located close to prosodic boundaries has mainly been studied for domain-final vowels, or for vowels in CV initial syllables of prosodic constituents to question the domain of initial strengthening [START_REF] Cho | Prosodic strengthening in transboundary V -to-V lingual movement in A merican English[END_REF], [START_REF] De | On the status of redundant features: the case of backing and rounding in A merican English[END_REF], [START_REF] Fougeron | V ariations articulatoires en début de constituants prosodiques de différents niveaux en français[END_REF].

In the present study, acoustic properties of vowels in absolute initial position (post-boundary) are compared in different prosodic domains. The question of the nature of the strengthening process is addressed by testing whether variation affects the acoustic cues of the phonological contrast between round and unrounded vowels and between vowels of different height in French. Fougeron [START_REF] Fougeron | A rticulatory properties of initial segments in several prosodic constituents in French[END_REF] showed that the vowel /i/ in the /p#ip/ sequence (where # stands for a prosodic boundary) has more linguopalatal contact (i.e. is more closed) following higher prosodic boundaries. This result challenges the view that prosodically-driven variations increase the sonority of vowels (sonority expansion hypothesis [START_REF] Beckman | Prosodic structure and tempo in a sonority model of articulatory dynamics[END_REF], [START_REF] Cho | Prosodic strengthening in transboundary V -to-V lingual movement in A merican English[END_REF], [START_REF] Fougeron | A rticulatory strengthening at edges of prosodic domains[END_REF]). It rather supports an enhancement of contrastive attributes in strong positions, as predicted by the local hyperarticulation hypothesis (e.g. [START_REF] De | On the status of redundant features: the case of backing and rounding in A merican English[END_REF], [START_REF] Fougeron | A rticulatory strengthening at edges of prosodic domains[END_REF]). Predictions of these two competing hypotheses will be further compared here by testing vowels contrasting in rounding (/i-y/, /e-ø/) and in vowel height (/i-e-a/).

MATERIAL AND METHODS

Speech material and speakers

Acoustic realization of the vowels /i, e, a, y, ø/ are studied in sentences cueing different prosodic phrasing. As illustrated in Table 1, the test vowels

(V 2 ) appear in a V 1 C 1 #V 2 C 2 sequence
where V 2 C 2 are the beginning of a (fake) first name, and # holds for either an Intonational Phrase (IP), an Accentual Phrase (AP), or a Word (Wd) boundary. Pre-boundary context (V 1 C 1 ) and following consonant (C 2 ) are kept constant (V 1 =/i/, C 1 & C 2 =/p/). Four native French female speakers (aged 25 to 40 years, with no identifiable regional accent) read the sentences in a random order. Each sentence was read twice consecutively at a comfortable speech rate and the entire list was read eight times, for a total of 1075 repetitions. The first author coded prosodic phrasing of each rendition, and for 180 utterances for which the phrasing was not obvious, the judgment of 3 expert phoneticians was used. 50 sentences were discarded due to lack of inter-judge agreement, leaving 1025 tokens to be considered in this study.

Rounding contrast 2.1. Acoustic cues of the rounding contrast

As expected and illustrated in Figure 1, acoustic differences between front rounded /y, ø/ and unrounded /i, e/ vowels are found in F3 and F2. All prosodic positions and all speakers pooled, the contrast between /i/ and /y/ is cued by a lower F2 (F(1,398)=1125, p<.0001) and a lower F3 (F(1,398)=1483, p<.0001) for the rounded vowel /y/. Similarly, the contrast between /e/ and /ø/ is also cued by a lower F2 (F(1,401)=1527, p<.0001), a lower F3 (F(1,401)=737, p<.0001) for /ø/. In terms of F3-F2 distance, /y/ show more compacity than /i/ (F(1,398)=1483, p<.0001). This compacity between F3-F2 is not found for the other rounded vowel /ø/, which has a rather central F2 and thus a larger F3-F2 distance than /e/ (F(1,401)=190, p<.0001). These acoustic differences are found to be all speaker-dependent (all interactions between vowel type and speakers at p<.0001). Subsequent analyses by speaker show that effect sizes vary, but that all speakers mark the distinctions found for F2 and F3. However, for the F3-F2 dimension, the distinction between /i/ and /y/ holds for 2 speakers only (S2, S3) and the distinction between /e/ and /ø/ holds for 3 speakers (all except S3).

Effect of prosodic position

The effect of 'Prosodic Position' is tested for each vowel (all speakers pooled) on the three acoustic dimensions found to be linked to the rounding contrast (F2, F3, F3-F2) with an ANOVA with the prosodic position as fixed factor. Results are presented in Table 2 and illustrated in Figure 2 for each speaker individually. Over all speakers, the effect of prosodic position on the realization of domain initial vowels varies according to the vowel, the acoustic dimensions affected and the number of position distinguished.

Distinctions are found to be more frequent and systematic between the IP initial position and the two lowest positions in the hierarchy, AP and Wd (except for F3-F2, for /ø/ where IP=AP). Vowels in AP initial position are generally similar to that in Wd initial position (except for F3 /e/ where the only 3 way distinction is found IP>AP>Wd).

Variations can be observed according to vowel category. Recall that the unrounded vowels /i/ and /e/ differ from their rounded counterparts by a higher F2 and a higher F3. Prosodic position is found to have an effect on these dimensions. Compared to AP and Wd-initial positions, IPinitial /i/s have a higher F3; and IP-initial /e/s have a higher F2 and F3. In terms of F3-F2, both unrounded vowels have a larger F3-F2 distance in IP initial position (the increase in F3 being greater than the change in F2).

These effects of prosodic positions on the realization of the two unrounded vowels vary according to speaker (all speaker*position interactions significant at p<.01). As illustrated on Figure 2, for /i/ the effects on F3 and F3-F2 hold for all speakers except S2. S3 also shows a higher F2 in IP-initial position than in AP/Wd. For /e/, higher F2 and F3 in IP-initial position is found for all speakers, while a larger F3-F2 distance holds for all speaker except S2.

As far as the rounded vowels are concerned, the effects of prosodic position are less consistent. All speakers pooled, /y/ is found to have a lower F3 and a smaller F3-F2 distance when initial in IP compared to AP/Wd. /ø/ has a higher F2 in IP initial position compared to AP/Wd and a smaller F3-F2 distances when initial in IP and AP compared to Wd. No effect is found on F3. Variation in F2 for this vowel could be linked to rounding but also to variation in tongue fronting.

For the rounded vowels, interactions between the factors position and speaker are also found (all at p<.01). Examination of the patterns by speakers (see Figure 2) shows that prosodic position effects on rounded vowels are less systematic across speakers, and are shown principally by speakers S1 and S4. For /y/, the effect on F3 holds only for S4, and the effect on F3-F2 holds for S1 and S4. For /ø/, the effects on F2 and F3-F2 hold for S1 and S4 (S3 showing also a distinction AP<*Wd for F3-F2)

Vowel height contrast

Acoustic differences between /i, e, a/

As expected, all speakers and all positions pooled, the vowel height contrast between /i, e, a/ is cued by variation in F1 (F(2,599)=1976, p<.0001). Posthoc comparisons (Fisher's test) show a three-way distinction with an increase of F1 from /i/ to /e/ to /a/ (p<.0001). While the size of the differences vary according to speaker (all interactions at p<.0001), all speakers show the same patterns, except for S4 who do not distinguish /i/ and /e/ in terms of F1.

Effect of prosodic position

The effect of the prosodic position was tested for each vowel on F1 values. All speakers pooled, an effect of position is found only for the open vowel /a/ (F(2,198)=82, p<.0001). Both all speakers confounded, and for all speaker tested individually, /a/ in IP-initial position have a higher F1 than in both AP and Wd initial positions.

An effect of position is also found for speaker S4 for /e/ (F(2,41=4), p=0.02), with IP-initial /e/ having a higher F1 (vs. AP/Wd). A surprising effect is observed for speaker S3 for /i/ with a higher F1 in IP and Wd initial positions compared to AP (F(2,42=9), p<.001). However, as illustrated on Figure 3, differences in Hz are rather small.

Conclusions

In this study, the effect of prosodic position on the realization of initial vowels is tested on acoustic features related to vowel contrasts in French. For rounded and unrounded vowels, variations in F2, F3, and F3-F2 are observed. For vowels contrasting in three degrees of aperture, variations in F1 are tested.

Results show that the prosodic position of the vowels can affect these acoustic dimensions. However, as observed previously in the literature (e.g. [2][7]), among the positions tested distinctions are more frequent and follow a more regular trend between the highest prosodic domain (here Intonational Phrase) on one side, and the lowest domains (here Accentual Phrase and Word) on the other side. Speaker differences are also frequently observed, some speakers showing more effects (e.g. S1 and S4 in 2.2) than others (e.g. S2 in 2.2).

Interestingly, our results show that the direction of the variations observed in IP initial positions compared to AP and Wd initial positions depend on vowel category. For the unrounded vowels /i/ and /e/, acoustic cues marking the contrast with their rounded counterparts are enhanced. These effects hold for at least three, if not all four speakers, depending on the comparisons. For the rounded vowels, there is also a global trend for an enhancement of the acoustic manifestation of rounding in IP initial position, but two speakers mainly support these effects. For the open vowel /a/, there is also an enhancement of the acoustic cue of vowel height (F1).

Globally, when an effect of prosodic position is found, /i, e/ are acoustically more 'unroundedlike', /y, ø/ are more 'rounded-like', and /a/ is more 'open-like' when they are in initial position of the the highest prosodic domain (IP). These results can be discussed in the light of the two hypotheses proposed in the literature to explain the nature of the influence of prosody on segmental realization. According to the sonority expansion hypothesis [START_REF] Beckman | Prosodic structure and tempo in a sonority model of articulatory dynamics[END_REF], segments in strong positions are expected to display an increase in their intrinsic sonority. The variations we observed in F2, F3, and F3-F2 can not be related to a sonority expansion. They rather support the predictions of local hyperarticulation hypothesis [START_REF] De | On the status of redundant features: the case of backing and rounding in A merican English[END_REF], according to which the contrastive characteristics of the unrounded and rounded vowels are enhanced by a local hyperarticulation of the gestures toward their presumed targets.

Concerning the effects found for /a/, the increase of F1 in IP initial position found for all speakers can be interpreted either as an increase of sonority, or as a hyperarticulation of the low target of the tongue. In order to distinguish these two interpretations, we were expected to find an effect on the F1 of the closed vowel /i/ or of the midclosed vowel /e/. Only one of the four speaker shows an effect for /i/, but the trend does not follow the prosodic hierarchy (IP,Wd>AP) and is hard to interpret. An effect is also found for one speaker for /e/, suggesting an increase in sonority. Further interpretation of these last results is not reasonable at this stage. This study needs to be completed by the observation of other French vowels (examination of the mid-open vowel // and back vowels is underway).

Figure 1 :

 1 Figure 1: Distribution of the vowels /i, y/ (top), /e, ø/ (bottom) in the F2/F3 dimensions (all speakers and all prosodic conditions confounded).

  182)=4 ns F(2,182)=32** F(2,182)=36** ns IP>*AP,Wd IP>*AP,Wd /e/ F(2,195)=13* F(2,195)=38** F(2,195)=6 * IP>*AP,W IP>*AP>*Wd IP>*AP,Wd /y/ F(2,200)=1 ns F(2,200)=11** F(2,200)=13** ns IP<*AP,Wd IP<*AP,Wd /ø/ F(2,190)=0 ns F(2,190)=3 ns F(2,190)=7 ** IP>*AP,W ns IP,AP<*Wd

Figure 2 :

 2 Figure 2: F2, F3 and F3-F2 values for /i/, /y/, /e/ and /ø/ for each speakers in IP, AP and Wd prosodic conditions (Hz).

Figure 3 :

 3 Figure 3: F1 values for /a/, /e/ and /i/for all speakers in IP, AP and Wd prosodic conditions (Hz).

Table 1 :

 1 Sample

		sentences containing the sequence
	V1C1#V2C2 with different prosodic boundaries (#:
	Intonational Phrase, Accentual Phrase, Word)
	IP	Pour le roi Philippe, Ipicaz et Larsen se battent comme
		des dieux. (For king P., Ipicaz and Larsen fight as
		gods).
	AP	Oupoulo, Philippe, Ipicaz et David sont comme des
		frères. (Oupoulo, Philippe, Ipicaz and David are like
		brothers).

Wd

Dalida et Lippe-Ipicaz sortent un tube d'enfer.

Dalida and Lippe-Ipicaz launch a hell of a hit.

Table 2 :

 2 Main effect of 'prosodic position' and Fisher post-hoc comparisons on F2, F3 and F2-F3 values for each vowel. Significance levels: **p<.0001, *p<.05, ns:non significant.